A time series data-driven agricultural ecosystem dynamic prediction and early warning method

Abstract

The application discloses a kind of based on timing data driving's agricultural ecosystem dynamic prediction and early warning method, belong to agricultural ecology monitoring technical field.The method first obtains remote sensing vegetation index, meteorology, soil physicochemical and crop yield multi-source long time series data, and generates multivariate equal length sequence matrix after preprocessing;Again, output reconstructed phase space state vector sequence by multivariate phase space reconstruction algorithm;Combined with sliding window entropy weight method and ridge regression, construct soil health degradation index SHDI;Then reconstruct sequence is constructed dynamic graph attention network with, predict future SHDI and quantile confidence interval;Finally, carry out double-path stability analysis, and the comprehensive result outputs three-level early warning.The application fuses multi-source time series information, accurately captures ecological evolution characteristics, effectively improves the accuracy of agricultural ecosystem prediction and early warning.

Description

Technical Field

[0001] This invention belongs to the field of agricultural ecosystem monitoring and intelligent early warning technology, specifically involving a method for dynamic prediction and early warning of agricultural ecosystems based on time-series data. Background Technology

[0002] In the rapeseed-maize rotation system, soil degradation is driven by a coupling of factors such as planting system, fertilization, and climate, exhibiting long-term nonlinear gradual-abrupt characteristics. Existing technologies, such as fuzzy comprehensive evaluation based on interval sampling, struggle to capture the time-delayed feedback between multiple variables. When using LSTM or Transformer for soil index prediction, the time-varying characteristics of the coupling strength between variables with the rotation process are not considered, and the inherent dynamic structure of the system cannot be revealed. Existing graph neural network methods mostly use static predefined graphs, which cannot adaptively reflect the changes in variable relationships during the rotation node transition period. Regarding stability analysis, existing early warning technologies mostly rely on single statistical thresholds, lacking a quantitative characterization of system dynamic stability loss, such as the Lyapunov exponent approaching zero, making it difficult to reliably capture precursors of degradation abrupt changes. There is no evidence of combining phase space reconstruction with time-varying graph attention mechanisms to quantitatively characterize the dynamic structure of agricultural ecosystems for prediction, nor is there evidence of using a dual-path complementary early warning mechanism combining Lyapunov exponent stability analysis and Bayesian online change point detection for advanced graded early warning of soil degradation. Summary of the Invention

[0003] The purpose of this invention is to provide a time-series data-driven method for dynamic prediction and early warning of agricultural ecosystems, addressing the limitations of existing methods in simultaneously extracting multivariable nonlinear dynamic structures, learning time-varying coupling relationships, and providing early warning of system stability degradation. By jointly reconstructing using PECUZAL to preserve the multivariable phase space dynamics, employing a time-varying graph attention network for adaptive learning of time-varying dependencies between variables, and combining Lyapunov stability analysis and Bayesian change point detection, a dual-path comprehensive early warning system is achieved. This enables reliable early warning before abrupt soil degradation in the rapeseed-maize rotation system.

[0004] To achieve the above objectives, the present invention adopts the following technical solution:

[0005] A time-series data-driven method for dynamic prediction and early warning of agricultural ecosystems includes the following steps:

[0006] Step S1: Obtain a multi-source long time series dataset, which includes remote sensing vegetation index sequence, meteorological sequence, soil physicochemical sequence and crop yield sequence. Preprocess each sequence in sequence to generate a multivariate equal-length sequence matrix.

[0007] Step S2: The multivariable phase space reconstruction joint embedding algorithm is used to perform nonlinear state space reconstruction on the multivariable sequences in the multivariable equal-length sequence matrix, mapping multiple one-dimensional observation sequences to the joint phase space, and outputting the reconstructed phase space state vector sequence.

[0008] Step S3: Based on the soil index data corresponding to the reconstructed phase space state vector sequence, the time-varying weight of each index is calculated using the sliding window entropy weight method, and the yield response coefficient of each index is obtained by ridge regression fitting of the soil index lagged by one season and the crop yield of the current season, and the Soil Health Degradation Index (SHDI) is constructed.

[0009] Step S4: Using the reconstructed phase space state vector sequence as input, a dynamic graph attention temporal prediction network is constructed to generate a time-varying adjacency matrix with adaptive variable interaction intensity. After aggregating node features through a multi-head graph attention mechanism, the SHDI prediction values ​​and quantile confidence intervals for multiple future time steps are output through a gated recurrent unit encoder and a direct multi-step output decoder.

[0010] Step S5: Conduct dual-path complementary stability analysis. Path A estimates the maximum Lyapunov exponent of the reconstructed phase space trajectory obtained from the reconstructed phase space state vector sequence to determine the stability of the temporal evolution. Path B performs Bayesian online change point detection on the SHDI prediction sequence composed of the SHDI prediction values ​​to determine abrupt changes. The combined dual-path analysis results output a three-level early warning.

[0011] Preferably, in step S2, the multivariate phase space reconstruction joint embedding algorithm sets the joint embedding dimension upper limit to 30 dimensions. In the iteration, the variable-delay combination that maximizes the reduction of the false nearest neighbor ratio is selected until the false nearest neighbor ratio is lower than 5% or the dimension upper limit is reached, and the output is a reconstructed phase space state vector sequence with a dimension not exceeding 30.

[0012] Preferably, in step S3, the sliding window entropy weight method applies a minimum weight constraint to the time-varying weights of the calculated indexes below 0.05 and renormalizes them; ridge regression is performed only on the training set, and the ridge parameters are determined through 10-fold cross-validation.

[0013] Preferably, the formula for calculating the Soil Health Degradation Index (SHDI) in step S3 is as follows: ;

[0014] In the formula, For the first Soil health degradation index at each time step; For the first The final entropy weight of each soil indicator; For the first Yield response coefficient of soil indicators; For the first Soil index in the first Normalized values ​​at each time step; This is a soil index with values ​​from 1 to 5, corresponding to organic matter, total nitrogen, available phosphorus, available potassium, and pH, respectively.

[0015] Preferably, in step S4, the time-varying adjacency matrix is ​​generated by scaling dot product attention calculation of node embedding, and after activation by scaling exponential linear units, entropy regularization is applied to prevent the adjacency matrix from degenerating and to retain variable interaction information.

[0016] Preferably, in step S4, the multi-head graph attention mechanism sets up 4 attention heads. Each attention head selects the 5 neighbors with the largest weights in the time-varying adjacency matrix to calculate the attention coefficients and aggregate them with weights. After the multi-head is spliced, the node features are activated by the exponential linear unit.

[0017] Preferably, in step S4, the gated cyclic unit encoder has a two-layer structure and a hidden dimension of 128. The direct multi-step output decoder outputs the median and 0.05 and 0.95 quantiles of the SHDI for the next 12 steps at once. The training loss includes mean square error, quantile loss and entropy regularization term.

[0018] Preferably, in step S5, path A uses a 96-step fixed-length sliding window and a 4-step sliding step size to estimate the maximum Lyapunov exponent. During the cold start period when the historical trajectory has less than 96 steps, path A is not activated and only path B performs the warning.

[0019] Preferably, in step S5, path B adopts a normal-inverse gamma conjugate prior. The initial hyperparameters are determined by the moment estimation of the first 48 steps of the training set. The risk rate is set to 1 / 250. A mutation alarm is triggered when the posterior probability of the change point is greater than 0.95 and the standard deviation of the residuals in 15 steps exceeds twice the baseline.

[0020] Preferably, in step S5, the three-level warning is a blue, orange, and red warning. The blue warning is when the SHDI is higher than the average of the same quarter in the past three years by one standard deviation for two consecutive periods. The orange warning is superimposed with the path stability judgment condition. The red warning is further superimposed with the change point and residual exceeding the limit condition.

[0021] The present invention provides a time-series data-driven method for dynamic prediction and early warning of agricultural ecosystems, which achieves several technical advantages:

[0022] This invention endows the system with dynamic structure through PECUZAL phase space reconstruction, enabling the attention of time-varying graphs to learn physically meaningful variable coupling relationships, thus solving the problem that existing black-box models cannot reflect the inherent evolutionary laws of agricultural ecosystems. The dual-path early warning mechanism of Lyapunov exponent and Bayesian change point detection collaboratively judges from two dimensions: deterministic stability degradation and probabilistic mutation, overcoming the shortcomings of single-path models in missing or false alarms during gradual-mutation processes. In the embodiments, the predicted RMSE is reduced by approximately 27% compared to static graph LSTM, and the early warning can be approximately 5 months ahead, with a significantly better synergistic effect than simple superposition of individual modules. Attached Figure Description

[0023] Figure 1 This is a flowchart of the method for dynamic prediction and early warning of agricultural ecosystems based on time-series data driven by the present invention.

[0024] Figure 2 This is a schematic diagram of the PECUZAL joint embedding reconstruction of the present invention.

[0025] Figure 3 This is the dual-path integrated early warning decision tree of the present invention. Detailed Implementation

[0026] The following detailed implementation of a time-series data-driven method for dynamic prediction and early warning of agricultural ecosystems according to the present invention is provided in conjunction with specific embodiments. These embodiments are for illustrative purposes only and are not intended to limit the scope of protection of the present invention.

[0027] Example 1: Implementation of a time-series data-driven method for dynamic prediction and early warning of agricultural ecosystems.

[0028] Combined with appendix Figures 1-3 As shown, this invention provides a method for dynamic prediction and early warning of agricultural ecosystems based on time-series data.

[0029] like Figure 1 As shown, Figure 1 This is a flowchart of the method of the present invention. The present invention provides a method for dynamic prediction and early warning of agricultural ecosystems based on time-series data, the method comprising the following steps:

[0030] Step S1: Construction and preprocessing of multi-source long time-series datasets.

[0031] Step S1 involves acquiring a multi-source long-time-series dataset, which includes remote sensing vegetation index sequences, meteorological sequences, soil physicochemical sequences, and crop yield sequences. Each sequence is preprocessed sequentially to generate a multivariate, equal-length sequence matrix. The specific process is as follows:

[0032] Using rapeseed-maize rotation fields as units, with a time step of 15 days, at least 7 complete rotation cycles were collected. Remote sensing was performed using Landsat- Surface reflectance products were obtained using LaSRC correction and CFMAS cloud masking. Meteorological data were acquired from meteorological stations in neighboring countries. Collected variables included: NDVI, MSAVI; effective accumulated temperature (baseline 10℃), cumulative rainfall, and sunshine duration; soil organic matter, total nitrogen, available phosphorus, available potassium, and pH, measured after each harvest and interpolated to 15 days using PCHIP; and rapeseed and maize grain yields. Median filtering (window 3), outlier removal (IQR 3x), and Min-Max normalization were performed sequentially to obtain the matrix. , , .

[0033] Step S2: Multivariable nonlinear state-space reconstruction.

[0034] Step S2 employs a multivariable phase space reconstruction joint embedding algorithm to perform nonlinear state space reconstruction on the multivariable sequences in the multivariable equal-length sequence matrix, mapping multiple one-dimensional observation sequences to the joint phase space and outputting a reconstructed phase space state vector sequence.

[0035] In step S2, the multivariate phase space reconstruction joint embedding algorithm sets the upper limit of the joint embedding dimension to 30. During iteration, it selects the variable-delay combination that maximizes the reduction of the false nearest neighbor ratio until the false nearest neighbor ratio is below 5% or the upper limit of the dimension is reached. The output is a reconstructed phase space state vector sequence with a dimension not exceeding 30. The specific process is as follows:

[0036] like Figure 2 As shown, Figure 2 This is a schematic diagram of the PECUZAL joint embedding reconstruction of the present invention. The PECUZAL joint embedding algorithm (i.e., the multivariate phase space reconstruction joint embedding algorithm) is used to optimize the embedding parameters for all variables within a unified iterative framework, preserving mutual information between variables. The joint embedding dimension is set to an upper limit of 30 dimensions (sample size 240, ratio 8:1). At each step, the variable that maximizes the reduction in the false nearest neighbor ratio—the delayed combination—is added to the reconstruction vector until the false nearest neighbor ratio is below 5% or the upper limit is reached. The experiment actually yielded a 28-dimensional joint reconstruction state vector. Forming a phase space trajectory matrix This reconstruction not only provides low-noise characteristics for downstream predictions, but also makes the system state evolution correspond to the identifiable attractor geometry in phase space, providing the necessary dynamic basis for the Lyapunov stability analysis in step S5.

[0037] Step S3: Construction of the Soil Health Degradation Index.

[0038] Step S3 uses the soil index data corresponding to the reconstructed phase space state vector sequence as a basis, calculates the time-varying weight of each index using the sliding window entropy weight method, and obtains the yield response coefficient of each index by ridge regression fitting of the soil index lagged by one season and the crop yield of the current season, and constructs the Soil Health Degradation Index (SHDI).

[0039] In step S3, the sliding window entropy weighting method applies a minimum weight constraint to the time-varying weights of the calculated indices below 0.05 and then renormalizes them; ridge regression is performed only on the training set, and the ridge parameters are determined through 10-fold cross-validation. The specific process is as follows:

[0040] Using organic matter, total nitrogen, available phosphorus, available potassium, and pH as the five core indicators, an entropy weighting method was employed for dynamic weighting. For each time step, data from the current time step and the previous 47 steps (out of a total of 48 steps) were used to calculate the information entropy of each indicator. ,in Entropy weight To prevent the weight of slowly changing indicators from degenerating to zero, the following measures are taken: Apply minimum weight constraints to the indicators 0.05, then renormalize to get the final weights. .

[0041] in, For the first The information entropy of a soil index measures the degree of disorder in the values ​​of the index within a 48-step window. The value range is [0,1]. The closer it is to 1, the smaller the change and the lower the information content. This is the index of the time step within the sliding window, and its value is the current step plus the previous 47 steps, a total of 48 consecutive time steps. For the first The indicator in the first The proportion of the value at each time step to the total of the 48 steps of this indicator satisfies... ; For the first Soil index in the first The values ​​at each time step after normalization preprocessing; For the first The sum of the values ​​of the soil index within the 48-step window; For the first The initial entropy weight of the indicator is obtained by normalizing the information entropy redundancy, reflecting the objective weighting result that the more drastic the change of the indicator within the window period, the greater the weight. Information entropy redundancy, i.e., the first The degree of variation of an indicator within a window period; the greater the change, the larger the value. This is the sum of the information entropy redundancy of the five indicators, used to convert redundancy into weights that sum to 1; The minimum initial entropy weight among the five indicators; The i-th after minimum weight constraint and renormalization The final entropy weight of the five indicators is 1.

[0042] Determination of yield response coefficient: Soil indicators after the previous season's harvest were used as independent variables, and the normalized yield of the current season's crop was used as the dependent variable. Ridge regression was performed only on the training set. The lagged independent variable eliminates the endogeneity problem caused by simultaneous bidirectional causality. Ridge parameters. The comprehensive response coefficient was obtained by averaging and normalizing the absolute values ​​of the standardized regression coefficients using 10-fold cross-validation. Organic matter 0.32, total nitrogen 0.30, available phosphorus 0.17, available potassium , Test set The calculation uses the training set. . The regularization parameter for ridge regression is selected on the training set through 10-fold cross-validation to suppress multicollinearity.

[0043] Defined as:

[0044] ;

[0045] The normalized values ​​are as follows: Standard value for excellent organic matter. Degradation warning value Excellent standard value for total nitrogen Degradation warning value Excellent standard value for available phosphorus Degradation warning value Excellent standard value for readily available potassium Degradation warning value All the above metrics are linearly scaled and truncated to [0,1].

[0046] In the above formula, For the first The comprehensive yield response coefficient (absolute value) of a soil index represents the contribution intensity of that index per unit of degradation to SHDI. For the first The soil health degradation index at each time step, with a value range of [0,1]. The larger the value, the more severe the soil degradation. This is a soil index with values ​​from 1 to 5, corresponding to organic matter, total nitrogen, available phosphorus, available potassium, and pH, respectively. For the first The final entropy weight of the indicator; For the first The absolute value of the output response coefficient of the item indicator; : No. The indicator in the first The normalized value of each time step, with a range of [0,1], where 0 represents the optimal state and 1 represents the degeneration warning value. For the first The indicator in the first The degree of normalization deviation at each time step represents the proportion of degradation relative to the optimal state; The numerator is the weighted sum of the contributions of the five indicators to the degradation. The more severe the degradation of each indicator, the larger its weight and response coefficient, the larger the numerator. The sum of the products of the five indicator weights and the response coefficients (denominator) is used as a normalization factor to ensure SHDI .

[0047] pH is normalized using segmented normalization:

[0048] ;

[0049] In the above formula, for pH value after time standardization; for The original measured soil pH value at any given time; where 4.5 is the warning value for severe acidification, 6.5 is the target / neutral suitable pH value, and 8.5 is the warning value for severe alkalization. SHDI The larger the value, the more severe the degradation.

[0050] Step S4, Dynamic Graph Attention Temporal Prediction Network.

[0051] Step S4 uses the reconstructed phase space state vector sequence as input to construct a dynamic graph attention temporal prediction network, generates a time-varying adjacency matrix of adaptive variable interaction intensity, aggregates node features through a multi-head graph attention mechanism, and outputs SHDI prediction values ​​and quantile confidence intervals for multiple future time steps through a gated recurrent unit encoder and a direct multi-step output decoder.

[0052] In step S4, the time-varying adjacency matrix is ​​generated by scaling dot product attention of node embedding. After activation by scaling exponential linear units, entropy regularization is applied to prevent the adjacency matrix from degenerating and to preserve variable interaction information.

[0053] In step S4, the multi-head graph attention mechanism sets up 4 attention heads. Each attention head selects the 5 neighbors with the largest weights in the time-varying adjacency matrix to calculate the attention coefficients and aggregate them with weights. After the multi-head is spliced, the node features are activated by the exponential linear unit.

[0054] In step S4, the gated recurrent unit encoder has a two-layer structure with a hidden dimension of 128. The direct multi-step output decoder outputs the median and 0.05 and 0.95 quantiles of the SHDI for the next 12 steps at once. The training loss includes mean squared error, quantile loss, and entropy regularization. The specific process is as follows:

[0055] (1) Time-varying adjacency matrix learner: Aggregate the current and previous 4 reconstructed vectors by variables and reduce the dimension to node embedding. (d=32, Xavier initialization). Adjacency matrix elements are calculated using scaled dot product attention:

[0056] ;

[0057] Applying entropy regularization , .

[0058] in, For the first The variable in the first... The node embedding vectors at each time step are obtained by dimensionality reduction through aggregation of the reconstructed vectors by variables. ; Dimension of the node embedding vector, with a value of 32; Xavier: parameter initialization method, which samples from a uniform distribution to make the input and output variances consistent, and is used to initialize the node embedding and the subsequent learnable weight matrix; This is the query weight matrix, used to map node embeddings to query vectors. ; This is the key weight matrix, used to map node embeddings to key vectors. ; and These are the numbers in the animated diagram. The node, the first Each node Feature embedding vector at time step; For the first Each time step node With nodes The original attention score between the query vector and the key vector is calculated by the scaled dot product of the query vector and the key vector; The scaling factor is the square root of the node embedding dimension, which scales the dot product result to prevent the softmax gradient from disappearing due to an excessively large dot product in high dimensions. For the first Each time step node With nodes The adjacency matrix element values ​​between variables represent the standardized interaction strength between the variables; The scaling exponential linear unit activation function performs a nonlinear transformation on the input to make its output approach zero mean and unit variance; To perform an exponential operation on the original attention score, any real value is mapped to a positive number; For nodes For all neighboring nodes The sum of the exponential attention scores is used for softmax normalization; This is the adjacency matrix entropy regularization loss term, which encourages the adjacency matrix to maintain a moderate level of sparsity and information content. This is the entropy regularization coefficient, with a value of 0.01, which controls the regularization strength. This is the information entropy of adjacent matrix elements. When the matrix tends to be uniformly distributed, this value is negative and has the largest absolute value. Adding a negative sign serves as a penalty term.

[0059] (2) Attention propagation in multi-head graphs: K=4 heads, the first head... Head attention coefficient:

[0060] ;

[0061] in , , Take the top-5 neighbors from the adjacency matrix. The node output is a concatenated multi-head ELU activation.

[0062] in, The number of attention heads is 4, and each head independently learns the feature representations of different subspaces. This is the attention head index, with values ​​ranging from 1 to 4; For the first Nodes in each attention head For nodes The normalized attention coefficient represents the neighbor's attention coefficient during information aggregation. Contribution weight; For the first The learnable parameter vector of each attention head This is used to map the concatenated node features to scalar attention scores; The hidden dimension inside the attention head, with a value of 32; For the first The linear transformation matrix of each attention head. Mapping input features to dimensional subspace; The dimension of the input features for the graph attention layer; For nodes The input feature vector of the current graph attention layer; For nodes The input feature vector of the current graph attention layer; This is a vector concatenation operator that joins two vectors end-to-end along a dimensional direction. To make nodes and nodes After the first The feature vectors after the head linear transformation are concatenated and used as input for the attention coefficient calculation; For the leakage-corrected linear unit activation function, a small positive slope is given to the negative input instead of setting it to zero directly; typically, the negative slope is 0.2. The scalar attention original score is obtained by taking the vector dot product of the transposed attention parameter vector and the concatenated features. For nodes The neighbor set, take the adjacency matrix. Zhongyu The top 5 connected nodes with the highest weight; Gathering for neighbors Traversal index in; For nodes The sum of the exponential attention scores of all neighbors is used for softmax normalization; ELU is the exponential linear unit activation function, which smoothly approaches the negative saturation value exponentially for negative inputs and is an identity mapping for positive inputs, used for the final activation after multi-head splicing. for The vectors output by the head attention are concatenated, connecting the node features calculated independently by each head end-to-end along the dimensional direction.

[0063] (3) Temporal Encoding and Multi-Step Output: The graph aggregation features are flattened and input into a two-layer GRU encoder (128 hiddens, Dropout 0.2). The decoder adopts a direct multi-step output (MIMO) structure, outputting the SHDI prediction values ​​and quantiles for the next H=12 steps at once. Quantile regression is jointly optimized by Pinball Loss, with quantiles q=0.05, 0.5, and 0.95. The total loss is the sum of MSE loss, the two quantile losses, and the entropy regularization term. Adam is trained for 200 rounds with a learning rate of 0.001 cosine annealing.

[0064] Step S5: Hierarchical early warning based on dynamic stability analysis.

[0065] Step S5 involves conducting a dual-path complementary stability analysis. Path A estimates the maximum Lyapunov exponent of the reconstructed phase space trajectory obtained from the reconstructed phase space state vector sequence to determine the stability of the temporal evolution. Path B performs Bayesian online change point detection on the SHDI prediction sequence composed of the SHDI prediction values ​​to determine abrupt changes. The combined dual-path analysis results output a three-level early warning.

[0066] In step S5, path A uses a 96-step fixed-length sliding window and a 4-step sliding step size to estimate the maximum Lyapunov exponent. During the cold start period when the historical trajectory has less than 96 steps, path A is not activated and only path B performs the warning.

[0067] In step S5 path B, a normal-inverse gamma conjugate prior is used. The initial hyperparameters are determined by the moment estimation of the first 48 steps of the training set. The risk rate is set to 1 / 250. A mutation alarm is triggered when the posterior probability of the change point is greater than 0.95 and the standard deviation of the residuals in 15 steps exceeds twice the baseline.

[0068] In step S5, the three-level early warning system consists of blue, orange, and red alerts. A blue alert is triggered when the SHDI (Shanghai Stock Exchange Analytical Index) exceeds the average of the same quarter over the past three years by one standard deviation for two consecutive periods. An orange alert is overlaid with path stability criteria, and a red alert is further overlaid with change point and residual exceedance criteria. The specific process is as follows:

[0069] like Figure 3 As shown, Figure 3 This is the dual-path integrated early warning decision tree of the present invention. Path A: Phase space Lyapunov exponent estimation. For the reconstructed trajectory... The Rosenstein small data method was used to estimate the maximum Lyapunov exponent. A fixed-length sliding window (96 steps, approximately 4 crop seasons) with a 4-step sliding step and a Theiler window with a 4-step step are used. The average slope of the logarithmic divergence rate curve of the reference point-nearest neighbor pair within each window is calculated as the maximum Lyapunov exponent. and standard deviation .when When three consecutive adjacent windows are valid, the system stability is determined to be significantly reduced. When the system cold start trajectory is less than 96 steps, path A is temporarily inactive, and path B drives the system.

[0070] in, The maximum Lyapunov exponent measures the average logarithmic divergence of neighboring point pairs in the phase space trajectory over time. A positive value indicates a chaotic system, while a negative value indicates that the system converges to a stable attractor. for The standard deviation of the slope obtained by linear regression fitting within the sliding window is used to quantify the local uncertainty of the Lyapunov exponent estimate. To determine the statistical criteria for a significant decrease in system stability, estimation uncertainty is taken into account. It is required that the divergence threshold condition is still met even when considering a twofold deviation, so as to reduce misjudgments caused by data fluctuations.

[0071] Path B: Bayesian Online Change Point Detection. The BOCPD (Bayesian Online Change Point Detection) framework of Adams & MacKay (2007) is adopted, using the student prediction distribution under the normal-inverse gamma conjugate prior, with hyperparameters updated online. Initial hyperparameters are estimated from the first 48 moments of the training set. Hazard rate When the posterior probability changes Furthermore, an alarm is triggered when the predicted residual 15-step moving standard deviation exceeds twice the standard deviation of the training period baseline.

[0072] in, The initial value of the prior sample size represents the estimate of the initial mean. The confidence level is set to 48. The initial parameter for the prior degrees of freedom, together with the scaling parameter of the inverse gamma distribution, controls the uncertainty of the variance, and is set to 24. The initial values ​​for the prior scale parameters are given by... Calculate the scaling parameters for the inverse gamma distribution; The variance of the first 48 steps of the SHDI sequence in the training set is calculated from the initial data segment using the method of moments. The Hazard Rate represents the prior probability of a change point occurring at any given time step. It has a value of 1 / 250, corresponding to an expected structural change occurring approximately every 250 steps (about 10 years).

[0073] Comprehensive early warning rules:

[0074] Blue alert: The SHDI is predicted to be higher than the average of the same quarter in the past three years by one standard deviation for two consecutive periods (30 days);

[0075] Orange alert: The blue condition is met, and path A is activated. Three consecutive windows are valid (when not active, the determination is based solely on criterion B).

[0076] Red alert: Orange conditions are met, and the probability of path B changing points is high. residual sliding standard deviation Double the baseline.

[0077] Example 2 is a simulation example based on the time-series data-driven method for dynamic prediction and early warning of agricultural ecosystems in Example 1.

[0078] Experimental data were generated by simulating typical paddy soil conditions for rapeseed-maize rotation in the middle and lower reaches of the Yangtze River. Soil parameters were referenced from the typical value range of Grade II paddy soil in the middle and lower reaches of the Yangtze River in the "Grade of Cultivated Land Quality" (GB / T 33469-2016), and meteorological parameters were referenced from historical meteorological statistics of Wuhan / Nanchang stations from 2014 to 2024. 240 time steps were used, with an 8:2 ratio for training / testing. Table 1 shows the results of the ablation experiment (mean ± standard deviation of 5 random seed trials):

[0079] Table 1

[0080]

[0081] Paired t-tests showed that each improvement was significant (p<0.05). PECUZAL reconstruction made the time-varying adjacency matrix of the dynamic graph attention learning correspond to the attractor topology, giving the model dynamics interpretability; it complements the dynamic graph attention, further improving prediction accuracy. In the fall of the second year of testing, SHDI increased from 0.42 to 0.67, and the system... Four consecutive windows >0, with a change point probability of 0.96, issued a red alert, about 5 months earlier than the usual monitoring alarm. The topsoil thickness was confirmed to have dropped to 12.5cm and the pH to 5.3.

[0082] Among them, DGATN: Dynamic Graph Attention Temporal Network.

[0083] DGATN_full: Refers to the complete method of this patent, namely the DGATN model that combines PECUZAL phase space reconstruction input and dynamic graph attention mechanism. It can be understood as: DGATN with Full PECUZAL-based input and dynamic graph attention.

[0084] DGATN_NoPEC: Refers to the control model in the ablation experiment that did not use PECUZAL reconstruction, but directly used the original 16-dimensional sequence as input. It can be understood as: DGATN without PECUZAL embedding.

[0085] PECUZAL: This is the name of a multivariate phase space reconstruction algorithm. Its core function is to optimize the embedding delay and dimension of all variables within a unified iterative framework, thereby achieving joint phase space reconstruction.

[0086] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0087] The above embodiments are merely illustrative examples and are not intended to limit the implementation. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is neither necessary nor possible to exhaustively list all possible implementations. However, obvious variations or modifications derived therefrom are still within the scope of protection of this invention.

Claims

1. A method for dynamic prediction and early warning of agricultural ecosystems based on time-series data, characterized in that, Includes the following steps: Step S1: Obtain a multi-source long time series dataset, which includes remote sensing vegetation index sequence, meteorological sequence, soil physicochemical sequence and crop yield sequence. Preprocess each sequence in sequence to generate a multivariate equal-length sequence matrix. Step S2: The multivariable phase space reconstruction joint embedding algorithm is used to perform nonlinear state space reconstruction on the multivariable sequences in the multivariable equal-length sequence matrix, mapping multiple one-dimensional observation sequences to the joint phase space, and outputting the reconstructed phase space state vector sequence. Step S3: Based on the soil index data corresponding to the reconstructed phase space state vector sequence, the time-varying weight of each index is calculated using the sliding window entropy weight method, and the yield response coefficient of each index is obtained by ridge regression fitting of the soil index lagged by one season and the crop yield of the current season, and the Soil Health Degradation Index (SHDI) is constructed. Step S4: Using the reconstructed phase space state vector sequence as input, a dynamic graph attention temporal prediction network is constructed to generate a time-varying adjacency matrix with adaptive variable interaction intensity. After aggregating node features through a multi-head graph attention mechanism, the SHDI prediction values ​​and quantile confidence intervals for multiple future time steps are output through a gated recurrent unit encoder and a direct multi-step output decoder. Step S5: Conduct dual-path complementary stability analysis. Path A estimates the maximum Lyapunov exponent of the reconstructed phase space trajectory obtained from the reconstructed phase space state vector sequence to determine the stability of the temporal evolution. Path B performs Bayesian online change point detection on the SHDI prediction sequence composed of the SHDI prediction values ​​to determine abrupt changes. The combined dual-path analysis results output a three-level early warning.

2. The method for dynamic prediction and early warning of agricultural ecosystems based on time-series data as described in claim 1, characterized in that, In step S2, the multivariate phase space reconstruction joint embedding algorithm sets the upper limit of the joint embedding dimension to 30. During the iteration, the variable-delay combination that maximizes the reduction of the false nearest neighbor ratio is selected until the false nearest neighbor ratio is lower than 5% or the upper limit of the dimension is reached. The output is a reconstructed phase space state vector sequence with a dimension not exceeding 30.

3. The method for dynamic prediction and early warning of agricultural ecosystems based on time-series data as described in claim 1, characterized in that, In step S3, the sliding window entropy weight method applies minimum weight constraints to the time-varying weights of the calculated indicators below 0.05 and renormalizes them; ridge regression is performed only on the training set, and the ridge parameters are determined through 10-fold cross-validation.

4. The method for dynamic prediction and early warning of agricultural ecosystems based on time-series data as described in claim 1, characterized in that, The formula for calculating the Soil Health Degradation Index (SHDI) in step S3 is as follows: ； In the formula, For the first Soil health degradation index at each time step; For the first The final entropy weight of each soil indicator; For the first Yield response coefficient of soil indicators; For the first Soil index in the first Normalized values ​​at each time step; This is a soil index with values ​​from 1 to 5, corresponding to organic matter, total nitrogen, available phosphorus, available potassium, and pH, respectively.

5. The method for dynamic prediction and early warning of agricultural ecosystems based on time-series data as described in claim 1, characterized in that, In step S4, the time-varying adjacency matrix is ​​generated by scaling dot product attention of node embedding. After activation by scaling exponential linear units, entropy regularization is applied to prevent the adjacency matrix from degenerating and to preserve variable interaction information.

6. The method for dynamic prediction and early warning of agricultural ecosystems based on time-series data as described in claim 1, characterized in that, In step S4, the multi-head graph attention mechanism sets up 4 attention heads. Each attention head selects the 5 neighbors with the largest weights in the time-varying adjacency matrix to calculate the attention coefficients and aggregate them with weights. After the multi-head is spliced, it is activated by an exponential linear unit to output node features.

7. The method for dynamic prediction and early warning of agricultural ecosystems based on time-series data as described in claim 1, characterized in that, In step S4, the gated recurrent unit encoder has a two-layer structure and a hidden dimension of 128. The direct multi-step output decoder outputs the median and 0.05 and 0.95 quantiles of the SHDI for the next 12 steps at once. The training loss includes mean square error, quantile loss and entropy regularization term.

8. The method for dynamic prediction and early warning of agricultural ecosystems based on time-series data as described in claim 1, characterized in that, Step S5: Path A uses a 96-step fixed-length sliding window and a 4-step sliding step size to estimate the maximum Lyapunov exponent. Path A is not activated during the cold start period when the historical trajectory has less than 96 steps, and only Path B performs the warning.

9. The method for dynamic prediction and early warning of agricultural ecosystems based on time-series data as described in claim 1, characterized in that, Step S5 path B adopts normal-inverse gamma conjugate prior. The initial hyperparameters are determined by the moment estimation of the first 48 steps of the training set. The risk rate is set to 1 / 250. A mutation alarm is triggered when the posterior probability of the change point is greater than 0.95 and the standard deviation of the residuals in 15 steps exceeds twice the baseline.

10. The method for dynamic prediction and early warning of agricultural ecosystems based on time-series data as described in claim 1, characterized in that, In step S5, the three-level warning is divided into blue, orange, and red. The blue warning is when the SHDI is higher than the average of the same quarter of the past three years by one standard deviation for two consecutive periods. The orange warning is superimposed with the path stability judgment condition. The red warning is further superimposed with the change point and residual exceeding the limit conditions.

Images