Soil degradation effect dynamic simulation method based on residual film-soil-crop causal chain

Abstract

The invention discloses a soil degradation effect dynamic simulation method based on a residual film-soil-crop causal chain, and belongs to the technical field of environmental system simulation and prediction. The problems of dynamic identification and process deduction of the residual agricultural film degradation effect are solved. According to the method, multi-dimensional fusion of agricultural film residue information data, soil physicochemical property information data, agricultural management and crop growth information data and climatic environment information data is realized through time dimension alignment, spatial dimension registration and submerged space feature extraction of multi-source data of an acquisition area, and multi-dimensional fusion of the agricultural film residue information data, the soil physicochemical property information data, the agricultural management and crop growth information data and the climatic environment information data is realized; forming a regional feature data set including a potential soil degradation index and a potential crop degradation index; constructing a residual film-soil-crop causal model based on a multilayer dynamic Bayesian network; dynamic updating of model parameters is realized through Bayesian flow; and performing residual agricultural film soil degradation effect dynamic simulation by using the trained residual film-soil-crop causal model based on the multilayer dynamic Bayesian network.

Description

Technical Field

[0001] This invention belongs to the field of environmental system simulation and prediction technology, specifically involving a dynamic simulation method for soil degradation effects based on the residual film-soil-crop causal chain. Background Technology

[0002] Residual agricultural film ages, breaks down, and accumulates in the soil, affecting soil physicochemical properties and crop growth, forming a typical causal chain effect of "residual film accumulation - soil degradation - crop yield reduction." There is a significant coupling relationship among residual film, soil, and crops. Residual film accumulation alters the soil physicochemical environment and weakens crop physiological functions, while crop growth status affects soil carbon and nitrogen cycling and the degradation rate of residual film. Identifying the causal relationship between residual agricultural film, soil degradation, and crop yield reduction, and dynamically analyzing degradation effects, identifying transmission pathways, and clarifying key driving factors, has become a critical bottleneck in the management of residual agricultural film.

[0003] Current mainstream models mostly employ structural equation regression, principal component regression, and simplified PLAM models, which generally suffer from two limitations: First, they tend to focus on local simulations of "residual agricultural film - soil degradation" or "residual agricultural film - crop yield reduction," failing to establish a unified modeling framework that integrates all three factors. Second, the model parameters are fixed, making it difficult to adaptively update them according to changes in regional characteristics, climate conditions, and management measures, resulting in decreased accuracy and insufficient applicability across different regions and scenarios. To overcome these limitations, there is an urgent need to establish a unified simulation method that unifies the causal chain of "residual film accumulation - soil degradation - crop yield reduction," possessing adaptive parameter updates and dynamic extrapolation capabilities. Summary of the Invention

[0004] The problem to be solved by this invention is to realize the dynamic identification and process deduction of the degradation effect of residual agricultural film, and to propose a dynamic simulation method of soil degradation effect based on the causal chain of residual film-soil-crop.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] A dynamic simulation method for soil degradation effects based on the residual film-soil-crop causal chain includes the following steps:

[0007] S1. Collect multi-source data from the region, including information on agricultural film residue, soil physicochemical properties, agricultural management and crop growth, and climate and environmental conditions.

[0008] S2. The regional multi-source data collected in step S1 is aligned in the time dimension, registered in the spatial dimension and extracted in the latent space to achieve multi-dimensional fusion of agricultural film residue information data, soil physicochemical property information data, agricultural management and crop growth information data and climate environment information data, forming a regional feature dataset including potential soil degradation index, potential crop degradation index and regulation function.

[0009] S3. Construct a causal model of residual film-soil-crop based on a multi-layer dynamic Bayesian network;

[0010] S4. Construct a dynamic parameter correction method for the residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network. Use the regional feature dataset obtained in step S2 to train the residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network to obtain a trained residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network.

[0011] S5. Dynamic simulation of soil degradation effects of residual agricultural film was performed using a trained causal model of residual film-soil-crop based on a multi-layer dynamic Bayesian network.

[0012] Furthermore, the specific implementation method of step S1 includes the following steps:

[0013] S1.1. Collect agricultural film residue information data based on multispectral, hyperspectral and shortwave infrared satellite imagery to obtain surface reflectance and raw data of each band, including agricultural film residue amount, aging degree, thickness and spatial distribution pattern;

[0014] S1.2. Collect soil physicochemical properties data by combining ground sampling, near-ground radar detection and laboratory testing, including soil porosity and bulk density, soil organic carbon content, pH, nitrogen, phosphorus and potassium nutrients, soil moisture and soil temperature;

[0015] S1.3. Collect agricultural management and crop growth information data, including vegetation index, crop type and mulch area, root length, rhizosphere biomass, crop yield, tillage method, irrigation system, fertilization intensity and crop rotation system;

[0016] S1.4. Collect climate and environmental information data, including regional temperature, precipitation, wind speed, sunshine, air humidity, freeze-thaw cycle frequency, and wind erosion index;

[0017] Furthermore, the specific implementation method of step S2 includes the following steps:

[0018] S2.1. Align the regional multi-source data collected in step S1 with the time dimension. Based on the typical phenological nodes of crops, dynamically interpolate the agricultural film residue information data, soil physicochemical property information data, agricultural management and crop growth information data, and climate environment information data in time so that the data from different years and observation sources are mapped to a unified phenological time axis.

[0019] S2.2. Perform spatial dimension registration on multi-source data in the region. Based on the boundaries of land parcels or regular grids, use ArcGIS zoning statistics and Kriging interpolation to achieve scale unification of data with different spatial resolutions and establish a unified spatial indexing system.

[0020] Then, the multi-source data from all regions were cleaned and standardized, and missing values ​​were filled with the dynamic mean to obtain the cleaned and standardized data matrix of the regional multi-source data. Where M is the amount of residual agricultural film after standardization treatment, A is the degree of aging after standardization treatment, D is the thickness after standardization treatment, P is the soil porosity after standardization treatment, ρ is the soil bulk density after standardization treatment, pH is the pH value of the soil after standardization treatment, C is the organic carbon content of the soil after standardization treatment, NPK is the nitrogen, phosphorus and potassium nutrients in the soil after standardization treatment, W is the soil moisture after standardization treatment, and NDVI is the vegetation index after standardization treatment. This refers to standardized crop growth information data other than NDVI.

[0021] S2.3. Perform latent space mapping and feature extraction on the data matrix after cleaning and standardization of regional multi-source data, including PCA preprocessing and autoencoder nonlinear mapping, to extract potential soil degradation index and potential crop degradation index;

[0022] S2.3.1. Calculate the covariance matrix and perform eigenvalue decomposition on the cleaned and standardized regional multi-source data matrix obtained in step S2.2. Select the top k eigenvectors whose cumulative variance contribution rate reaches a preset threshold to form a projection matrix, which is used to linearly transform the cleaned and standardized regional multi-source data matrix into linearly dimensionality-reduced eigenvalues. ;

[0023] S2.3.2. Features after linear dimensionality reduction Based on this, an autoencoder model is built. The autoencoder model consists of an encoder and a decoder. The encoder transmits data through several fully connected hidden layers. Mapping to latent space vector The hidden layer uses ReLU activation, and the output layer uses linear activation to form a stable representation; the decoder performs a reverse mapping of the latent space vectors to reconstruct the input, and the reconstruction error is used during training. Let be the objective function, and obtain the model parameters through iterative optimization, where For autoencoders of linear dimensionality reduction features The reconstructed output results;

[0024] S2.3.3. After the autoencoder model is trained, the output latent space vector... As a fusion feature, and through linear mapping, potential soil degradation index and potential crop degradation index are formed, expressed as follows:

[0025]

[0026]

[0027] in, As a potential soil degradation index, For use in latent space vectors Mapped to potential soil degradation index The weight vector, For the corresponding bias term, As a potential crop degradation index, For use in latent space vectors Mapped to potential crop degradation index The weight vector, For the corresponding bias term;

[0028] S2.4. Combine climate and environmental information data with agricultural management information to form an external control function. This is used to correct the transmission strength of degradation paths under different regional conditions, resulting in a regional feature dataset. t represents time.

[0029] Furthermore, the specific implementation method of step S3 includes the following steps:

[0030] S3.1. Design the model structure, using the potential soil degradation index and the potential crop degradation index as the core, and establish a causal main chain. , Let be the set of characteristics of residual agricultural film at time t, including the amount of residual agricultural film. aging degree and thickness ;

[0031] Extended subchains were constructed by combining soil porosity, bulk density, organic carbon content, pH, nitrogen, phosphorus and potassium nutrients, moisture, temperature, root length, rhizosphere biomass, NDVI and yield.

[0032] Each node in a multi-layer dynamic Bayesian network H-DBNs The conditional probability is expressed in parameterized form, as follows:

[0033] in, For each node The conditional probability, For the set of parent nodes, These are parameters for the local condition model;

[0034] The local conditional model employs a linear or lightweight nonlinear structure, comprising one or two fully connected hidden layers, and is expressed as follows:

[0035]

[0036] in, This represents the state value of the i-th node at time t. For its parent node set, For the parent node to The connection weight matrix, For bias terms, The activation function is ReLU or Sigmoid; the activation function of the hidden layer is ReLU or Sigmoid, and the output layer uses a parameterized representation of the conditional distribution formed by a linear mapping.

[0037] S3.2. Constructing a time extrapolation mechanism for a residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network: dividing the observation period into discrete time slices. Establish time dependencies based on the Markov assumption:

[0038]

[0039] in, The set of variables representing all nodes at time t. The set of variables represents time t-1; the model dynamically extrapolates along the time dimension through parameter updates, for each time step... Using Bayesian update, we get:

[0040]

[0041] in, Represents the set of model parameters. This represents the sequence of observed data from time 1 to time t. Let be the posterior distribution of the parameter. Let be the likelihood function of the data at the current time. Let be the posterior distribution of the parameters at the previous time step;

[0042] S3.3. Construct a spatial constraint mechanism for a causal model of residual agricultural film-soil-crop based on a multi-layer dynamic Bayesian network. In the spatial dimension, a spatial weight matrix is ​​introduced to capture the diffusion and transmission effects of residual agricultural film, soil conditions, and crop growth among neighboring units. The expression is as follows:

[0043]

[0044] Where W is the spatial weight matrix, w ij This indicates that units i and j are spatially adjacent. ;

[0045] right State introduces spatial propagation terms:

[0046]

[0047] Where j represents a neighboring unit that has a spatial proximity relationship with the i-th node, and n is the total number of spatial units. Spatial weighting coefficient, This represents the state value of adjacent unit j at time t-1, used to reflect the relationship between adjacent regions. Spatial propagation effects;

[0048] S3.4. Embedding the external regulation function into the conditional probability of the residual film-soil-crop causal model based on a multilayer dynamic Bayesian network, the expression is obtained as follows:

[0049]

[0050] in, This represents the external control function composed of climate and environmental information and agricultural management information, where W is the spatial weight matrix. for Local condition model parameters;

[0051] The joint distribution representation of the residual film-soil-crop causal model based on a multilayer dynamic Bayesian network is obtained as follows:

[0052] ;

[0053] in, This represents the sequence of all node states from time 1 to time T, where T is the time length of the deduction. The prior distribution of the initial state of the model;

[0054] S3.5. Set structural constraints for the residual film-soil-crop causal model based on multilayer dynamic Bayesian network, use the mechanism of agricultural film influence as a priori constraint to limit the causal direction, and score and screen candidate network structures through Bayesian Information Criterion (BIC) and Bayesian Factor (BF) to determine the final causal structure of the residual film-soil-crop causal model based on multilayer dynamic Bayesian network.

[0055] Furthermore, the specific implementation method of step S4 includes the following steps:

[0056] S4.1. Establish an adaptive prior update method, setting the model to update the posterior distribution of parameters after completing one causal inference at time step t. It is automatically used as the prior distribution input for the next time step t+1 to carry out continuous parameter evolution and information transmission;

[0057] The dataset D corresponding to the next time step t+1 t+1 After input, the model generates the likelihood function for the next time step t+1. The parameter distribution is dynamically corrected using a Bayesian update formula, expressed as:

[0058]

[0059] in, To provide a model with given parameters Spatial weight matrix and external control function Dataset under conditions The likelihood function;

[0060] S4.2. Establish a parameter correction and weight adjustment method, introduce an equivalent sample size correction mechanism, dynamically adjust the prior strength through the predicted residuals, and use Markov chain Monte Carlo (MCMC) sampling to generate parameter sample sequences. Calculate the posterior mean and confidence interval To quantitatively describe the trend of parameter changes and quantify the uncertainty;

[0061] S4.3. Establish a feedback explanation and dynamic attribution method. The SHAP value is calculated independently as the explanation layer. The SHAP explanation feedback mechanism is embedded in the Bayesian flow update process. The Model-Agnostic interpreter is selected, and the causal main chain is analyzed through the SHAP value. The node contributions are quantitatively explained to achieve a dynamic attribution of soil degradation effects through reverse extrapolation.

[0062] Furthermore, the specific implementation method of step S5 includes the following steps:

[0063] S5.1. Using a pre-trained causal model of residual film-soil-crop based on a multi-layer dynamic Bayesian network, through the causal main chain... Causal inversion was performed on crop phenotypic observation data, and the inversion expression is as follows:

[0064]

[0065]

[0066] in, For crop phenotypic observation data at time t, As a potential soil degradation index, This represents the set of residual agricultural film characteristics at the corresponding time point. For the set of model parameters, This is a Bayesian inversion operator based on posterior probability maximization, used to invert soil degradation status and upstream residual film driving factors under given crop phenotypic observation data;

[0067] S5.2. The model under given control parameters Under these conditions, forward extrapolation is performed on future time slices to generate a continuous evolutionary sequence of residual film-soil-crop, yielding:

[0068]

[0069] in, This represents the predicted potential crop degradation index at time t+k. For model parameters The determined forward inference function, W is the spatial weight matrix.

[0070] The beneficial effects of this invention are:

[0071] This invention presents a dynamic simulation method for soil degradation effects based on a causal chain of residual agricultural film, soil, and crop. It constructs a spatiotemporal causal structure of "residual film accumulation - soil degradation - crop yield reduction," enabling dynamic identification and process deduction of degradation effects. Combined with a Bayesian flow parameter update mechanism, the model parameters can be adaptively corrected with time and regional characteristics, enhancing the stability and generalization ability of causal inference. Furthermore, through a causal inversion algorithm, it quantitatively identifies key driving factors and transmission paths of the degradation process, revealing the dynamic response laws of soil degradation and crop yield reduction induced by residual agricultural film. This method overcomes the limitations of fragmented causal chains and static parameters in traditional models, achieving systematic dynamic modeling and precise analysis, providing methodological support for the analysis and scientific management of the correlation between residual agricultural film-induced soil degradation and crop yield reduction. Attached Figure Description

[0072] Figure 1This is a flowchart of a dynamic simulation method for soil degradation effects based on the residual film-soil-crop causal chain, as described in this invention.

[0073] Figure 2 This is an analytical diagram of the driving factors of soil degradation and crop yield reduction in regions A and B of the present invention, where (a) represents region A and (b) represents region B;

[0074] Figure 3 The diagram shows the dynamic correction of key parameters of the model in different regions over time, where (a) represents region A and (b) represents region B.

[0075] Figure 4 The following diagrams show the simulation results of the residual film-soil-crop degradation effect in different regions according to the present invention. (a) shows the distribution of the potential soil degradation index in region A, (b) shows the distribution of the potential crop degradation index in region A, (c) shows the distribution of the potential soil degradation index in region B, and (d) shows the distribution of the potential crop degradation index in region B.

[0076] Figure 5 This is an analytical diagram showing the contribution of environmental factors to the residual film-induced soil degradation effect of the present invention.

[0077] Figure 6 This is a structural block diagram of a dynamic simulation method for soil degradation effects based on the residual film-soil-crop causal chain described in this invention. Detailed Implementation

[0078] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are only for explaining the invention and are not intended to limit the invention; that is, the described specific embodiments are merely a part of the embodiments of the invention, and not all of them. The components of the specific embodiments of the invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations, and the invention may also have other embodiments.

[0079] Therefore, the following detailed description of specific embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected specific embodiments of the invention. All other specific embodiments obtained by those skilled in the art based on these specific embodiments without inventive effort are within the scope of protection of this invention.

[0080] To further understand the invention's content, features, and effects, the following specific embodiments are provided, along with accompanying drawings. Figure 1 - Appendix Figure 6 Detailed explanation is as follows:

[0081] Example 1:

[0082] A dynamic simulation method for soil degradation effects based on the residual film-soil-crop causal chain includes the following steps:

[0083] S1. Collect multi-source data from the region, including information on agricultural film residue, soil physicochemical properties, agricultural management and crop growth, and climate and environmental conditions.

[0084] Furthermore, the specific implementation method of step S1 includes the following steps:

[0085] S1.1. Collect agricultural film residue information data based on multispectral, hyperspectral and shortwave infrared satellite imagery to obtain surface reflectance and raw data of each band, including agricultural film residue amount, aging degree, thickness and spatial distribution pattern;

[0086] Furthermore, the amount of residual agricultural film is used as input with the reflectance of short-wave infrared, red edge and visible light bands to construct an existing agricultural film identification index. Combined with the measured residual film amount in ground quadrats, the index is calculated using a multiple linear regression or machine learning inversion model and statistically analyzed at the plot scale.

[0087] The aging degree A is constructed by using indicators such as peak shift and reflectance decrease of agricultural film spectral curves in hyperspectral and shortwave infrared bands. The aging index of agricultural film is obtained by combining the measured spectral data of film samples at different aging stages through empirical formulas or principal component regression inversion.

[0088] Thickness D is obtained by fitting the existing film thickness inversion formula with the measured film thickness and image reflectance or transmittance as independent variables in the short-wave infrared sensitive band, and then using this inversion model to calculate the remote sensing pixels.

[0089] The spatial distribution pattern is initially distinguished using the brightness features of remote sensing images, and the film-covered areas are identified by combining existing supervised classification methods. This yields a spatial distribution map of agricultural film and its coverage intensity classification results, which are used to characterize the coverage intensity and degradation characteristics of agricultural film in different regions and seasons.

[0090] S1.2. Collect soil physicochemical properties data by combining ground sampling, near-ground radar detection and laboratory testing, including soil porosity and bulk density, soil organic carbon content, pH, nitrogen, phosphorus and potassium nutrients, soil moisture and soil temperature;

[0091] Furthermore, soil porosity and bulk density can be selected based on the availability of regional data, or data fusion can be performed in typical areas. Based on Sentinel-1 / SAR multipolar scattering characteristics or ground-based radar (GPR) inversion dielectric constant, bulk density and porosity are estimated using existing soil physical models to form a temporally varying soil structure estimate. Using long-term stable bulk density, texture, and organic matter data from the national soil database, regional rasterized bulk density is obtained through geostatistical interpolation, according to P=1-ρ / ρ s Calculate porosity to construct the regional baseline background field for soil structure.

[0092] Organic carbon (C), pH, and nitrogen, phosphorus, and potassium (NPK) nutrients were analyzed using the continuous spectral characteristics of hyperspectral imagery in the visible-near-infrared band. The spectral sensitivity of organic carbon and major soil nutrients was analyzed, and the spatial distribution of organic carbon (C) and NPK nutrients was estimated using existing spectral inversion models. Soil pH was retrieved based on the color index and texture features of UAV or satellite RGB imagery, enabling the extraction of soil chemical information in large-scale areas lacking sampling points.

[0093] Soil moisture and temperature were obtained through inversion methods using satellite microwave and radar remote sensing. Surface soil moisture was retrieved using passive microwave brightness temperature, active microwave scattering, and multipolar radar information provided by satellites such as SMAP, Fengyun (FY), and Sentinel-1. The inversion results were then input into a land surface process model, and regional-scale soil moisture and soil temperature were calculated based on energy budget, water balance, and soil heat transfer mechanisms, forming continuous time-series monitoring data.

[0094] S1.3. Collect agricultural management and crop growth information data, including vegetation index, crop type and mulch area, root length, rhizosphere biomass, crop yield, tillage method, irrigation system, fertilization intensity and crop rotation system;

[0095] Furthermore, crop type and mulch area can be directly obtained through classification and identification using UAV or high-resolution satellite imagery.

[0096] The vegetation index NDVI is calculated from red and near-infrared reflectance using existing formulas. Root length and rhizosphere biomass are directly measured at typical sampling points; in large-scale areas, biomass can be inverted based on LAI or vegetation growth models. Crop yield is inverted based on actual measurements at typical sampling points, and at the regional scale, it is inverted using remote sensing yield estimation models (light energy utilization efficiency models, vegetation index models, etc.). Management information such as tillage methods, irrigation systems, fertilization intensity, and crop rotation systems are verified temporally by combining agricultural sector databases and farmer records, as well as irrigation events detected by remote sensing.

[0097] S1.4. Collect climate and environmental information data, including regional temperature, precipitation, wind speed, sunshine, air humidity, freeze-thaw cycle frequency, and wind erosion index;

[0098] Furthermore, regional temperature Precipitation (R), wind speed (V), solar radiation (Ls), and air humidity (H) are obtained directly or calculated from radiation products. Freeze-thaw cycle frequency is based on temperature sequences and statistically analyzed according to existing freeze-thaw criteria, counting the number of times the temperature crosses 0°C. Wind erosion index. The calculations were performed using existing empirical wind erosion models based on wind speed, wind direction, and surface roughness.

[0099] S1.5. Collect monitoring data on the surface and subsurface environment, including changes in surface elevation, fluctuations in groundwater level, and surface runoff information.

[0100] Furthermore, surface elevation changes are directly calculated using GNSS control points, multi-temporal DEM, or InSAR data. Groundwater level fluctuations are obtained through remote sensing gravity inversion and hydrological model estimation. Groundwater storage changes are acquired based on satellite-based time-series water storage change data, and this result is input into a regional hydrological model to estimate groundwater depth and its time-varying sequence based on the water balance process. Surface runoff information is simulated at a regional scale using a hydrological model, combined with rainfall data, to obtain the intensity of non-point source runoff.

[0101] S2. The regional multi-source data collected in step S1 is aligned in the time dimension, registered in the spatial dimension and extracted in the latent space to achieve multi-dimensional fusion of agricultural film residue information data, soil physicochemical property information data, agricultural management and crop growth information data and climate environment information data, forming a regional feature dataset including potential soil degradation index, potential crop degradation index and regulation function.

[0102] To address the challenge of uniformly representing multi-source heterogeneous data across time and space, this paper proposes a method for constructing regional feature datasets based on latent space mapping. This method achieves multi-dimensional fusion of climate, soil, crop, and agricultural film information through temporal alignment, spatial registration, and latent space feature extraction. This results in a high-quality input dataset capable of characterizing regional differences and degradation dynamics, providing a structurally consistent feature foundation for causal inference models.

[0103] Furthermore, the specific implementation method of step S2 includes the following steps:

[0104] S2.1. Align the regional multi-source data collected in step S1 with the time dimension. Based on the typical phenological nodes of crops, dynamically interpolate the agricultural film residue information data, soil physicochemical property information data, agricultural management and crop growth information data, and climate environment information data in time so that the data from different years and observation sources are mapped to a unified phenological time axis.

[0105] S2.2. Perform spatial dimension registration on multi-source data in the region. Based on the boundaries of land parcels or regular grids, use ArcGIS zoning statistics and Kriging interpolation to achieve scale unification of data with different spatial resolutions and establish a unified spatial indexing system.

[0106] Then, the multi-source data from all regions were cleaned and standardized, and missing values ​​were filled with the dynamic mean to obtain the cleaned and standardized data matrix of the regional multi-source data. Where M is the amount of residual agricultural film after standardization treatment, A is the degree of aging after standardization treatment, D is the thickness after standardization treatment, P is the soil porosity after standardization treatment, ρ is the soil bulk density after standardization treatment, pH is the pH value of the soil after standardization treatment, C is the organic carbon content of the soil after standardization treatment, NPK is the nitrogen, phosphorus and potassium nutrients in the soil after standardization treatment, W is the soil moisture after standardization treatment, and NDVI is the vegetation index after standardization treatment. This refers to standardized crop growth information data other than NDVI.

[0107] Furthermore, outliers are removed, missing values ​​are filled with dynamic means, and standardization processes include Z-score standardization for continuous variables, Min-Max normalization for ratio indicators, and Box-Cox transformation for skewed distribution variables to improve data accuracy and comparability.

[0108] S2.3. Perform latent space mapping and feature extraction on the data matrix after cleaning and standardization of regional multi-source data, including PCA preprocessing and autoencoder nonlinear mapping, to extract potential soil degradation index and potential crop degradation index;

[0109] S2.3.1. Calculate the covariance matrix and perform eigenvalue decomposition on the cleaned and standardized regional multi-source data matrix obtained in step S2.2. Select the top k eigenvectors whose cumulative variance contribution rate reaches a preset threshold to form a projection matrix, which is used to linearly transform the cleaned and standardized regional multi-source data matrix into linearly dimensionality-reduced eigenvalues. ;

[0110] S2.3.2. Features after linear dimensionality reduction Based on this, an autoencoder model is built. The autoencoder model consists of an encoder and a decoder. The encoder transmits data through several fully connected hidden layers. Mapping to latent space vector The hidden layer uses ReLU activation, and the output layer uses linear activation to form a stable representation; the decoder performs a reverse mapping of the latent space vectors to reconstruct the input, and the reconstruction error is used during training. Let be the objective function, and obtain the model parameters through iterative optimization, where For autoencoders of linear dimensionality reduction features The reconstructed output results;

[0111] S2.3.3. After the autoencoder model is trained, the output latent space vector... As a fusion feature, and through linear mapping, potential soil degradation index and potential crop degradation index are formed, expressed as follows:

[0112]

[0113]

[0114] in, As a potential soil degradation index, For use in latent space vectors Mapped to potential soil degradation index The weight vector, For the corresponding bias term, As a potential crop degradation index, For use in latent space vectors Mapped to potential crop degradation index The weight vector, For the corresponding bias term;

[0115] S2.4. Combine climate and environmental information data with agricultural management information to form an external control function. This is used to correct the transmission strength of degradation paths under different regional conditions, resulting in a regional feature dataset. t represents time.

[0116] S3. Constructing a causal model of residual film-soil-crop based on multilayer dynamic Bayesian networks; To characterize the soil degradation and crop growth restriction processes driven by residual film accumulation, this invention constructs a spatiotemporal causal model based on multilayer dynamic Bayesian networks (H-DBNs). The model consists of a three-layer structure: an agricultural film layer, a soil degradation layer, and a crop growth layer, with each layer's nodes corresponding to the residual film characteristics. Potential soil degradation index and potential crop degradation index .

[0117] Furthermore, the specific implementation method of step S3 includes the following steps:

[0118] S3.1. Design the model structure, using the potential soil degradation index and the potential crop degradation index as the core, and establish a causal main chain. , Let be the set of characteristics of residual agricultural film at time t, including the amount of residual agricultural film. aging degree and thickness ;

[0119] Extended subchains were constructed by combining soil porosity, bulk density, organic carbon content, pH, nitrogen, phosphorus and potassium nutrients, moisture, temperature, root length, rhizosphere biomass, NDVI and yield.

[0120] Each node in a multi-layer dynamic Bayesian network H-DBNs The conditional probability is expressed in parameterized form, as follows:

[0121] in, For each node The conditional probability, For the set of parent nodes, These are parameters for the local condition model;

[0122] The local conditional model employs a linear or lightweight nonlinear structure, comprising one or two fully connected hidden layers, and is expressed as follows:

[0123]

[0124] in, This represents the state value of the i-th node at time t. For its parent node set, For the parent node to The connection weight matrix, For bias terms, The activation function is ReLU or Sigmoid; the activation function of the hidden layer is ReLU or Sigmoid, and the output layer uses a parameterized representation of the conditional distribution formed by a linear mapping.

[0125] S3.2. Constructing a time extrapolation mechanism for a residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network: dividing the observation period into discrete time slices. Establish time dependencies based on the Markov assumption:

[0126]

[0127] in, The set of variables representing all nodes at time t. The set of variables represents time t-1; the model dynamically extrapolates along the time dimension through parameter updates, for each time step... Using Bayesian update, we get:

[0128]

[0129] in, Represents the set of model parameters. This represents the sequence of observed data from time 1 to time t. Let be the posterior distribution of the parameter. Let be the likelihood function of the data at the current time. Let be the posterior distribution of the parameters at the previous time step;

[0130] S3.3. Construct a spatial constraint mechanism for a causal model of residual agricultural film-soil-crop based on a multi-layer dynamic Bayesian network. In the spatial dimension, a spatial weight matrix is ​​introduced to capture the diffusion and transmission effects of residual agricultural film, soil conditions, and crop growth among neighboring units. The expression is as follows:

[0131]

[0132] Where W is the spatial weight matrix, w ij This indicates that units i and j are spatially adjacent. ;

[0133] right State introduces spatial propagation terms:

[0134]

[0135] Where j represents a neighboring unit that has a spatial proximity relationship with the i-th node, and n is the total number of spatial units. Spatial weighting coefficient, This represents the state value of adjacent unit j at time t-1, used to reflect the relationship between adjacent regions. Spatial propagation effects;

[0136] S3.4. Embedding the external regulation function into the conditional probability of the residual film-soil-crop causal model based on a multilayer dynamic Bayesian network, the expression is obtained as follows:

[0137]

[0138] in, This represents the external control function composed of climate and environmental information and agricultural management information, where W is the spatial weight matrix. for Local condition model parameters;

[0139] The joint distribution representation of the residual film-soil-crop causal model based on a multilayer dynamic Bayesian network is obtained as follows:

[0140] ;

[0141] in, This represents the sequence of all node states from time 1 to time T, where T is the time length of the deduction. The prior distribution of the initial state of the model;

[0142] S3.5. Set structural constraints for the residual film-soil-crop causal model based on multilayer dynamic Bayesian network, use the mechanism of agricultural film influence as a priori constraint to limit the causal direction, and score and screen candidate network structures through Bayesian Information Criterion (BIC) and Bayesian Factor (BF) to determine the final causal structure of the residual film-soil-crop causal model based on multilayer dynamic Bayesian network.

[0143] S4. Construct a dynamic parameter correction method for the residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network. Use the regional feature dataset obtained in step S2 to train the residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network to obtain the trained residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network.

[0144] To address the shortcomings of traditional Bayesian models, which suffer from "static priors, fixed parameters, and delayed updates" leading to insufficient long-term learning capabilities, a dynamic parameter update method based on Bayesian flow is proposed. This method constructs a closed-loop evolutionary system of "adaptive prior update - dynamic parameter correction - posterior validation feedback," enabling continuous self-learning and temporal convergence of model parameters, and ensuring the model's stability and generalization ability under varying regional, climatic, and management conditions.

[0145] Furthermore, the specific implementation method of step S4 includes the following steps:

[0146] S4.1. Establish an adaptive prior update method, setting the model to update the posterior distribution of parameters after completing one causal inference at time step t. It is automatically used as the prior distribution input for the next time step t+1 to carry out continuous parameter evolution and information transmission;

[0147] The dataset D corresponding to the next time step t+1 t+1 After input, the model generates the likelihood function for the next time step t+1. The parameter distribution is dynamically corrected using a Bayesian update formula, expressed as:

[0148]

[0149] in, To provide a model with given parameters Spatial weight matrix and external control function Dataset under conditions The likelihood function;

[0150] S4.2. Establish a parameter correction and weight adjustment method, introduce an equivalent sample size correction mechanism, dynamically adjust the prior strength through the predicted residuals, and use Markov chain Monte Carlo (MCMC) sampling to generate parameter sample sequences. Calculate the posterior mean and confidence interval To quantitatively describe the trend of parameter changes and quantify the uncertainty;

[0151] S4.3. Establish a feedback explanation and dynamic attribution method. The SHAP value is calculated independently as the explanation layer. The SHAP explanation feedback mechanism is embedded in the Bayesian flow update process. The Model-Agnostic interpreter is selected, and the causal main chain is analyzed through the SHAP value. The node contributions are quantitatively explained to achieve a dynamic attribution of soil degradation effects through reverse extrapolation.

[0152] S5. Dynamic simulation of soil degradation effects of residual agricultural film was performed using a trained multi-layer dynamic Bayesian network-based residual film-soil-crop causal model.

[0153] Furthermore, the specific implementation method of step S5 includes the following steps:

[0154] S5.1. Utilize the trained residual film-soil-crop causal model based on a multilayer dynamic Bayesian network through the causal main chain. The causal inversion is performed, and the inversion expression is as follows:

[0155]

[0156]

[0157] in, For crop phenotypic observation data at time t, As a potential soil degradation index, This represents the set of residual agricultural film characteristics at the corresponding time point. For the set of model parameters, This is a Bayesian inversion operator based on posterior probability maximization, used to invert soil degradation status and upstream residual film driving factors under given crop phenotypic observation data;

[0158] Furthermore, input crop phenotypic observation data Through Bayesian inversion operator Obtain the potential soil degradation index and upstream driving factors The contribution interval, combined with the spatial weight matrix By analyzing the spatial diffusion characteristics of degradation sources, a causal inversion can be achieved from "crop restriction" to "soil degradation" and "residual film accumulation".

[0159] S5.2. The model under given control parameters Under these conditions, forward extrapolation is performed on future time slices to generate a continuous evolutionary sequence of residual film-soil-crop, yielding:

[0160]

[0161] in, This represents the predicted potential crop degradation index at time t+k. For model parameters The determined forward inference function, W is the spatial weight matrix.

[0162] Furthermore, the final output will consist of three core types:

[0163] Dynamic analysis of degradation effect: through degradation effect index and The temporal changes describe the degree and rate of soil and crop degradation and their response to external disturbances, forming a characterization of the degradation effects in different regions.

[0164] Key driver identification: combining The contribution interval and SHAP interpretation help identify the impact. , The main climate or management factors are identified, and their relative contributions and directions of action are quantified to clarify the core driving forces of degradation effects in different regions.

[0165] Degradation pathway identification: based on and The posterior distribution and conditional probability of the degradation effect were used to assess the actual transmission intensity of the "residual film-soil-crop" main chain in different regions and time periods, revealing the transmission path of degradation effect in the causal chain.

[0166] The following is a practical application example of this embodiment:

[0167] The study areas were selected as two typical loess regions in western my country (A, 10km × 10km) and black soil regions in northern China (B, 10km × 10km), and each was divided into 100 1km × 1km plots (numbered A1-A100, B1-B100). The specific implementation process is illustrated using these two regions as examples:

[0168] 1. Dynamic acquisition of multi-source data in the region

[0169] A unified spatial indexing system was established in regions A and B, dividing the study area into plot-level grid units. Using each plot as the basic observation unit, agricultural production cycle data covering the entire process of spring sowing, summer management, and autumn harvest were collected for 2022-2023.

[0170] Land surface illumination, temperature, vegetation index (NDVI), and spectral characteristics of agricultural film were obtained using Sentinel-2, GF-6, Landsat-8, and MODIS remote sensing images. The residual amount of agricultural film (M), aging index (A), and thickness (D) were retrieved by combining shortwave infrared and hyperspectral data to characterize the spatial distribution and seasonal variation of agricultural film.

[0171] Temperature data was obtained using data from the China Meteorological Administration, regional automatic weather stations, and the ERA5-Land reanalysis product. ,precipitation Wind speed ,illumination ,humidity Freeze-thaw frequency With wind erosion index Climate factors, such as those mentioned above, are used to characterize the driving effect of external environmental disturbances on degradation processes.

[0172] Soil porosity was obtained through ground sampling, near-ground radar (GPR), and SoilGrids data. , bulk density , Organic carbon , , moisture Physicochemical parameters were used, and information on groundwater level, redox properties, and microbial abundance was supplemented by the National Soil Environmental Monitoring and Research Database.

[0173] By leveraging drone aerial photography and agricultural remote sensing monitoring, combined with local agricultural databases, data on crop types and mulch area were extracted. Irrigation intensity Fertilization intensity Farming methods Including agricultural management and crop growth indicators such as yield.

[0174] 2. Construction of Regional Feature Dataset

[0175] In this embodiment, to ensure the consistency and fusionability of multi-source data in both temporal and spatial dimensions, a latent space mapping method is used to construct a regional feature dataset. This process includes four steps: temporal alignment, spatial registration, data standardization, and latent space feature extraction.

[0176] Time dimension alignment: Climate, soil, crop, and agricultural film observation indicators are segmented according to the typical phenological stages of major crops (sowing, jointing, grain filling, and maturity). Missing time series data are filled using dynamic interpolation to ensure comparability of data from different years and observation sources on a unified time axis.

[0177] Spatial Dimension Registration: A unified spatial indexing system was established based on plot boundaries or a regular grid (100 m × 100 m). Regional statistics, kriging interpolation, and spatial resampling were performed using the ArcGIS Pro platform to unify the spatial resolution of multi-source data. To improve accuracy, all variables were cleaned and standardized: outliers were removed, missing data were filled with the time-moving mean, continuous variables were standardized using Z-score, ratio indicators were normalized using Min–Max, and skewed distribution variables were corrected using Box-Cox transformation.

[0178] Latent Space Feature Extraction and Fusion: Based on the Python environment, PCA (Principal Component Analysis) and an autoencoder were used to perform nonlinear dimensionality reduction and feature fusion on the cleaned multi-source data to extract the potential soil degradation index Z, which reflects the comprehensive degradation degree of soil physicochemical properties, and the potential crop degradation index Y, which characterizes the degree of crop growth restriction. The two types of latent variables output by the model serve as the core inputs for subsequent causal modeling.

[0179] Regional regulation parameter embedding: mapping regional climate, surface conditions, and agricultural management factors into regulation functions. The data is then input into subsequent causal models to dynamically correct the differences in the intensity of the "residual film-soil-crop" transmission pathway in different regions.

[0180] 3. Construction of a Causal Model of Residual Film-Soil-Crop Based on H-DBNs

[0181] In this embodiment, to achieve dynamic causal analysis between residual agricultural film, soil physicochemical properties, and crop growth status, a regionalized dynamic causal inference model is constructed using multilayer dynamic Bayesian networks (H-DBNs). The model comprehensively considers temporal evolution and spatial diffusion mechanisms to perform spatiotemporal causal inference on the three-layer system of "residual film-soil-crop".

[0182] First, the observational data for regions A and B from 2022-2023 were divided chronologically into a training set (70%) and a test set (30%) to ensure the temporal continuity of the samples and the model's generalization performance. The model uses the potential soil degradation index Z... t With potential crop degradation index Y t Construct the main causal chain using the core variable. Through conditional dependencies The temporal transition of the variables is characterized, and a spatial weight matrix is ​​introduced. Describe the spatial diffusion effects between land parcels. Regional climate and management parameters. Embedded in conditional probabilities, the transmission intensity is dynamically adjusted under different ecological zones to achieve regional adaptation. To verify the scientific validity of the model structure, Bayesian factor (BF) structure scores were performed on 84 causal chains in regions A and B. The results show that the BF values ​​of main chains I–III in both regions are greater than 2.5, significantly better than the other sub-chains, indicating that they are reasonable in terms of statistical confidence and agronomic logic. Table 1 shows the key causal chains of soil degradation-crop yield reduction in region A, and Table 2 shows the key causal chains of soil degradation-crop yield reduction in region B.

[0183] Table 1

[0184]

[0185] Table 2

[0186]

[0187] Based on the causal chain structure score, the Pearson correlation coefficient matrix of the main variables of the regional transmission path is calculated, and a correlation heatmap is plotted as follows: Figure 2 As shown, the rationality of the model path was verified. The results showed that in region A, agricultural film residue was stably negatively correlated with soil porosity, water content, NDVI, and yield, reflecting a continuous degradation chain of "porosity obstruction - weakened water supply - restricted crop growth - yield decline". In region B, freeze-thaw frequency and wind erosion intensity were strongly correlated with soil bulk density, organic carbon, and crop activity, showing a degradation pattern of "structural compaction - fertility decline - growth inhibition - yield loss".

[0188] 4. Dynamic updating of model parameters based on Bayesian flow

[0189] This embodiment constructs a Bayesian flow model based on the PyTorch framework in a Python environment, achieving continuous parameter evolution and adaptive updates. During model runtime, iterative updates of the posterior parameters are performed in continuous time series with a time step Δt = 1 month. The main path parameters include β. MZ (Agricultural film residue → soil degradation conduction intensity), β ZY (Soil degradation → Crop growth inhibition intensity) and β MY (Agricultural film residue → intensity of direct crop stress), the three correspond to the core causal chains I-III respectively.

[0190] First, a prior distribution of parameters was established using data from the spring sowing season of 2022 as the initial sample. And generate an initial parameter set through Bayesian sampling. Subsequently, an explicit time-series advancement strategy with a time step of Δt = 1 month was adopted. Observational data was progressively input into the continuous time series, the likelihood function was calculated, and the posterior distribution was iteratively corrected according to the Bayesian update formula, ensuring that the prior-likelihood-posterior distribution remained continuously flowing in the time dimension. During the optimization process, the torch.optim.AdamW optimizer was used to minimize the path prediction error and the negative log-likelihood of the posterior path parameters. Dynamic adjustments are made, and an equivalent sample size weighting mechanism is used to suppress the interference of outliers on the update direction. Gradient calculation and backpropagation are performed using the `torch.autograd` automatic differentiation mechanism, ensuring that the model balances data fitting accuracy and physical consistency during training. Finally, the model automatically calculates the posterior mean of the parameters at each time step. and confidence interval It adaptively adjusts the parameter intensity according to environmental disturbances (sunlight, freeze-thaw, wind erosion, etc.) to achieve self-learning generalization and dynamic convergence under different regional and climatic conditions.

[0191] Figure 3 This demonstrates the dynamic evolution of causal chain path parameters over time. In region A ( Figure 3 a) The rapid increase during the spring sowing season and the maintenance of high values ​​in summer indicate that residual film aging and high temperatures together accelerate the soil degradation process. The significant increase from July to September indicates that soil degradation has a stronger impact on yield formation during the critical growth period of crops. Overall, the values ​​were low, with only a slight increase during the autumn harvest season, reflecting the lagged effects of cumulative stress. Conditional probability and The synchronous increase over time indicates that the "residual film accumulation-soil degradation" pathway shows a synergistic strengthening trend over time.

[0192] In region B ( Figure 3 b), β MZ It rises significantly during the spring freeze-thaw period, followed by periodic fluctuations, demonstrating the surface fragmentation effect under the combined effects of freeze-thaw disturbance and wind erosion. It peaks during the summer heat and water concentration period and then declines, reflecting the seasonal superposition of structural damage and rhizosphere stress. It remains at a low level, consistent with the typical pattern of "structural failure-driven degradation." (Conditional probability) The soil condition in this region initially rose and then fluctuated, indicating that repeated freeze-thaw cycles and wind erosion drove the periodic degradation of the soil.

[0193] Table 3 shows the performance comparison between the dynamic Bayesian causal model of this invention and the traditional static Bayesian model. The results show that the model of this invention significantly outperforms the traditional model in terms of accuracy, stability, and convergence efficiency. Dynamic Model The MSE was increased to 0.993 (approximately 4% improvement over the static model), while the MAPE decreased to 0.027 and from approximately 9.1% to 6.5%. It decreased from 0.0093 to 0.0079. Both Cs and Cs showed significant improvement.

[0194] Table 3

[0195]

[0196] 4. Causal inversion and path identification of the soil degradation effect of residual agricultural film

[0197] To achieve causal attribution and critical path identification in the degradation process, an inversion operator is constructed based on a trained dynamic Bayesian causal model. Crop phenotypic observation data Perform reverse engineering to obtain the potential soil degradation index. and upstream driving factor group The posterior contribution range is determined. The inversion process is implemented in a Python environment, using PyTorch and PyMC3 frameworks for joint modeling. By generating a posterior sequence of parameters through Bayesian inversion and MCMC sampling, the transmission chain of "crop yield reduction - soil degradation - residual film accumulation" is recovered in time, realizing the reverse identification and dynamic analysis of degradation paths.

[0198] Figure 4 The spatial distribution characteristics of the potential soil degradation index Z and potential crop degradation index Y in regions A and B from 2022 to 2023 are presented. Both Z and Y are derived from the model based on multi-source observation data and directly reflect the model's core output on degradation effects: Z comprehensively reflects the physical and chemical degradation state, such as decreased soil porosity, increased bulk density, decreased organic carbon, and weakened water supply; Y reflects the degree of crop growth restriction, such as decreased root activity, weakened vegetation vigor (NDVI), and hindered yield formation.

[0199] In region A (Fig. 4a, b), Z significantly increased in the southwest and central areas, indicating soil degradation phenomena such as pore blockage, insufficient water supply, and accelerated decomposition of organic matter driven by high sunlight and large mulch area. Correspondingly, Y increased synchronously in the same area, manifesting as decreased NDVI and crop yield reduction, verifying the causal transmission chain of "residual film-soil-crop" and its spatiotemporal consistency revealed by the model. In region B (Fig. 4c, d), Z exhibited a "band-clod" composite distribution dominated by freeze-thaw disturbance and wind erosion intensity, reflecting a typical degradation pattern of structural fragmentation caused by freeze-thaw cycles, surface loss due to wind erosion, and increased bulk density. Y was synchronously restricted around this high-value zone, reflecting crop degradation effects such as root zone stress, decreased vegetation vitality, and increased yield sensitivity, perfectly consistent with the dominant environmental factors (freeze-thaw and wind erosion) identified by the model.

[0200] To enhance the interpretability of the inversion results, this invention introduces the SHAP explanation mechanism in the model's posterior inference stage to quantitatively decompose the conditional contributions of each node in the causal main chain. The model calls shap.Explainer to process the input feature matrix. Calculate feature contribution values ​​and implement them based on the Model-Agnostic framework. The dynamic interpretation allows for a direct and intuitive demonstration of the relative contributions of different environmental factors at the plot scale and their direction of influence.

[0201] like Figure 5 As shown, the SHAP contribution results reveal the dominant roles of different environmental factors in the regional degradation process. In region A, the duration of sunlight... With the area of ​​the film The positive contribution was the highest (+0.32 and +0.12), indicating that long-term sunlight accelerates the aging and embrittlement of agricultural film, and the accumulation of residual film in the cultivated layer leads to decreased porosity and enhanced decomposition of organic matter, thereby limiting root growth; in region B, the wind erosion index... With freeze-thaw frequency The presence of +0.08 and +0.06 as the main positive driving factors indicates that wind erosion and freeze-thaw disturbances exacerbate surface structure fragmentation and reduce aggregate stability, leading to increased bulk density and decreased permeability, ultimately affecting the vertical transport of water and roots. These results demonstrate that this invention can accurately identify degradation mechanisms based on the characteristics of climate disturbances in different regions, validating the model's dynamic source tracing and causal discrimination capabilities at cross-regional scales.

[0202] It should be noted that relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0203] Although this application has been described above with reference to specific embodiments, various modifications can be made and components can be replaced with equivalents without departing from the scope of this application. In particular, as long as there is no structural conflict, the features in the specific embodiments disclosed in this application can be combined with each other in any way. The lack of an exhaustive description of these combinations in this specification is merely for the sake of brevity and resource conservation. Therefore, this application is not limited to the specific embodiments disclosed herein, but includes all technical solutions falling within the scope of the claims.

Claims

1. A method for dynamic simulation of soil degradation effects based on the residual film-soil-crop causal chain, characterized in that, Includes the following steps: S1. Collect multi-source data from the region, including information on agricultural film residue, soil physicochemical properties, agricultural management and crop growth, and climate and environmental conditions. S2. Using the regional multi-source data collected in step S1, through temporal alignment, spatial registration, and latent spatial feature extraction, multidimensional fusion of agricultural film residue information, soil physicochemical property information, agricultural management and crop growth information, and climate environment information is achieved, forming a data set including a potential soil degradation index, a potential crop degradation index, and a regulatory function. Regional feature dataset; S3. Construct a causal model of residual film-soil-crop based on a multi-layer dynamic Bayesian network; S4. Construct a dynamic parameter correction method for the residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network. Use the regional feature dataset obtained in step S2 to train the residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network to obtain a trained residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network. S5. Dynamic simulation of soil degradation effects of residual agricultural film was performed using a trained causal model of residual film-soil-crop based on a multi-layer dynamic Bayesian network.

2. The method for dynamic simulation of soil degradation effects based on the residual film-soil-crop causal chain according to claim 1, characterized in that, The specific implementation method of step S1 includes the following steps: S1.

1. Collect agricultural film residue information data based on multispectral, hyperspectral and shortwave infrared satellite imagery to obtain surface reflectance and raw data of each band, including agricultural film residue amount, aging degree, thickness and spatial distribution pattern; S1.

2. Collect soil physicochemical properties data by combining ground sampling, near-ground radar detection and laboratory testing, including soil porosity and bulk density, soil organic carbon content, pH, nitrogen, phosphorus and potassium nutrients, soil moisture and soil temperature; S1.

3. Collect agricultural management and crop growth information data, including vegetation index, crop type and mulch area, root length, rhizosphere biomass, crop yield, tillage method, irrigation system, fertilization intensity and crop rotation system; S1.

4. Collect climate and environmental information data, including regional temperature, precipitation, wind speed, sunshine, air humidity, freeze-thaw cycle frequency, and wind erosion index.

3. The method for dynamic simulation of soil degradation effects based on the residual film-soil-crop causal chain according to claim 2, characterized in that, The specific implementation method of step S2 includes the following steps: S2.

1. Align the regional multi-source data collected in step S1 with the time dimension. Based on the typical phenological nodes of crops, dynamically interpolate the agricultural film residue information data, soil physicochemical property information data, agricultural management and crop growth information data, and climate environment information data in time so that the data from different years and observation sources are mapped to a unified phenological time axis. S2.

2. Perform spatial dimension registration on multi-source data in the region. Based on the boundaries of land parcels or regular grids, use ArcGIS zoning statistics and Kriging interpolation to achieve scale unification of data with different spatial resolutions and establish a unified spatial indexing system. Then, the multi-source data from all regions were cleaned and standardized, and missing values ​​were filled with the dynamic mean to obtain the cleaned and standardized data matrix of the regional multi-source data. Where M is the amount of residual agricultural film after standardization treatment, A is the degree of aging after standardization treatment, D is the thickness after standardization treatment, P is the soil porosity after standardization treatment, ρ is the soil bulk density after standardization treatment, pH is the pH value of the soil after standardization treatment, C is the organic carbon content of the soil after standardization treatment, NPK is the nitrogen, phosphorus and potassium nutrients in the soil after standardization treatment, W is the soil moisture after standardization treatment, and NDVI is the vegetation index after standardization treatment. This refers to standardized crop growth information data, excluding NDVI. S2.

3. Perform latent space mapping and feature extraction on the data matrix after cleaning and standardization of regional multi-source data, including PCA preprocessing and autoencoder nonlinear mapping, to extract potential soil degradation index and potential crop degradation index; S2.3.

1. Calculate the covariance matrix and perform eigenvalue decomposition on the cleaned and standardized regional multi-source data matrix obtained in step S2.

2. Select the top k eigenvectors whose cumulative variance contribution rate reaches a preset threshold to form a projection matrix, which is used to linearly transform the cleaned and standardized regional multi-source data matrix into linearly dimensionality-reduced eigenvalues. ; S2.3.

2. Features after linear dimensionality reduction Based on this, an autoencoder model is built. The autoencoder model consists of an encoder and a decoder. The encoder transmits data through several fully connected hidden layers. Mapping to latent space vector The hidden layer uses ReLU activation, and the output layer uses linear activation to form a stable representation; The decoder performs a reverse mapping of the latent space vectors to reconstruct the input, and during training, it uses the reconstruction error as a basis for calculation. Let be the objective function, and obtain the model parameters through iterative optimization, where For autoencoders of linear dimensionality reduction features The reconstructed output results; S2.3.

3. After the autoencoder model is trained, the output latent space vector... As a fusion feature, and through linear mapping, potential soil degradation index and potential crop degradation index are formed, expressed as follows: ； ； in, As a potential soil degradation index, For use in latent space vectors Mapped to potential soil degradation index The weight vector, For the corresponding bias term, As a potential crop degradation index, For use in latent space vectors Mapped to potential crop degradation index The weight vector, For the corresponding bias term; S2.

4. Combine climate and environmental information data with agricultural management information to form an external control function. This is used to correct the transmission strength of degradation paths under different regional conditions, resulting in a regional feature dataset. t represents time.

4. The method for dynamic simulation of soil degradation effects based on the residual film-soil-crop causal chain according to claim 3, characterized in that, The specific implementation method of step S3 includes the following steps: S3.

1. Design the model structure, using the potential soil degradation index and the potential crop degradation index as the core, and establish a causal main chain. , This is the set of characteristics of residual agricultural film at time t, including the amount of residual agricultural film. aging degree and thickness ; Extended subchains were constructed by combining soil porosity, bulk density, organic carbon content, pH, nitrogen, phosphorus and potassium nutrients, moisture, temperature, root length, rhizosphere biomass, NDVI and yield. Each node in a multi-layer dynamic Bayesian network H-DBNs The conditional probability is expressed in parameterized form, as follows: in, For each node The conditional probability, For the set of parent nodes, These are parameters for the local condition model; The local conditional model employs a linear or lightweight nonlinear structure, comprising one or two fully connected hidden layers, and is expressed as follows: ； in, This represents the state value of the i-th node at time t. For its parent node set, For the parent node to The connection weight matrix, For bias terms, The activation function is ReLU or Sigmoid; the activation function of the hidden layer is ReLU or Sigmoid, and the output layer uses a parameterized representation of the conditional distribution formed by a linear mapping. S3.

2. Constructing a time extrapolation mechanism for a residual film-soil-crop causal model based on a multi-layer dynamic Bayesian network: dividing the observation period into discrete time slices. Establish time dependencies based on the Markov assumption: ； in, The set of variables representing all nodes at time t. The set of variables represents time t-1; the model dynamically extrapolates along the time dimension through parameter updates, for each time step... Using Bayesian update, we get: ； in, Represents the set of model parameters. This represents the sequence of observed data from time 1 to time t. Let be the posterior distribution of the parameter. Let be the likelihood function of the data at the current time. Let be the posterior distribution of the parameters at the previous time step; S3.

3. Construct a spatial constraint mechanism for a causal model of residual agricultural film-soil-crop based on a multi-layer dynamic Bayesian network. In the spatial dimension, a spatial weight matrix is ​​introduced to capture the diffusion and transmission effects of residual agricultural film, soil conditions, and crop growth among neighboring units. The expression is as follows: ； Where W is the spatial weight matrix, w ij This indicates that units i and j are spatially adjacent. ; right State introduces spatial propagation terms: ； Where j represents a neighboring unit that has a spatial proximity relationship with the i-th node, and n is the total number of spatial units. Spatial weighting coefficient, This represents the state value of adjacent unit j at time t-1, used to reflect the relationship between adjacent regions. Spatial propagation effects; S3.

4. Embedding the external regulation function into the conditional probability of the residual film-soil-crop causal model based on a multilayer dynamic Bayesian network, the expression is obtained as follows: ； in, This represents the external control function composed of climate and environmental information and agricultural management information, where W is the spatial weight matrix. for Local condition model parameters; The joint distribution representation of the residual film-soil-crop causal model based on a multilayer dynamic Bayesian network is obtained as follows: ； in, This represents the sequence of all node states from time 1 to time T, where T is the time length of the deduction. The prior distribution of the initial state of the model; S3.

5. Set structural constraints for the residual film-soil-crop causal model based on multilayer dynamic Bayesian network, use the mechanism of agricultural film influence as a priori constraint to limit the causal direction, and score and screen candidate network structures through Bayesian Information Criterion (BIC) and Bayesian Factor (BF) to determine the final causal structure of the residual film-soil-crop causal model based on multilayer dynamic Bayesian network.

5. The method for dynamic simulation of soil degradation effects based on the residual film-soil-crop causal chain according to claim 4, characterized in that, The specific implementation method of step S4 includes the following steps: S4.

1. Establish an adaptive prior update method, setting the model to update the posterior distribution of parameters after completing one causal inference at time step t. It is automatically used as the prior distribution input for the next time step t+1 to carry out continuous parameter evolution and information transmission; The dataset D corresponding to the next time step t+1 t+1 After input, the model generates the likelihood function for the next time step t+1. The parameter distribution is dynamically corrected using a Bayesian update formula, expressed as: ； in, To provide a model with given parameters Spatial weight matrix and external control function Dataset under conditions The likelihood function; S4.

2. Establish a parameter correction and weight adjustment method, introduce an equivalent sample size correction mechanism, dynamically adjust the prior strength through the predicted residuals, and use Markov chain Monte Carlo (MCMC) sampling to generate parameter sample sequences. Calculate the posterior mean and confidence interval To quantitatively describe the trend of parameter changes and quantify the uncertainty; S4.

3. Establish a feedback explanation and dynamic attribution method. The SHAP value is calculated independently as the explanation layer. The SHAP explanation feedback mechanism is embedded in the Bayesian flow update process. The Model-Agnostic interpreter is selected, and the causal main chain is analyzed through the SHAP value. The node contributions are quantitatively explained to achieve a dynamic attribution of soil degradation effects through reverse extrapolation.

6. The method for dynamic simulation of soil degradation effects based on the residual film-soil-crop causal chain according to claim 5, characterized in that, The specific implementation method of step S5 includes the following steps: S5.

1. Using a pre-trained causal model of residual film-soil-crop based on a multi-layer dynamic Bayesian network, through the causal main chain... Causal inversion was performed on crop phenotypic observation data, and the inversion expression is as follows: ； ； in, For crop phenotypic observation data at time t, As a potential soil degradation index, This represents the set of residual agricultural film characteristics at the corresponding time point. For the set of model parameters, This is a Bayesian inversion operator based on posterior probability maximization, used to invert soil degradation status and upstream residual film driving factors under given crop phenotypic observation data; S5.

2. The model under given control parameters Under these conditions, forward extrapolation is performed on future time slices to generate a continuous evolutionary sequence of residual film-soil-crop, yielding: ； in, This represents the predicted potential crop degradation index at time t+k. For model parameters The determined forward inference function, W is the spatial weight matrix.

Images