NDVI time sequence prediction method combining physical equation constraint and long short-term memory network

Abstract

The invention discloses a physical equation constraint and long short-term memory network combined NDVI (normalized difference vegetation index) time sequence prediction method, which comprises the following steps of: acquiring meteorological variables, vegetation indexes and soil data related to NDVI change in a historical time period of a target area from a plurality of data sources, and preprocessing to obtain a spatio-temporal data set; extracting statistical features and spatial geometric features of an NDVI time sequence from the spatio-temporal data set, and constructing an NDVI spatio-temporal cube of which the dimension is a time-space-environment factor; the space-time cube is flattened after being subjected to seasonal trend decomposition, a trained mixed type multivariable LSTM-PINN model is input for time sequence prediction, and an NDVI time sequence prediction result of the target area in the target time period is obtained. The invention relates to the technical field of NDVI (normalized difference vegetation index) time sequence prediction, and realizes high-precision prediction of NDVI time sequence change of a target area by fusing multiple technical means such as multi-source three-dimensional space-time combination construction, physical-data dual-drive LSTM-PINN network modeling, constraint condition construction based on a vegetation ecological mechanism and the like.

Description

Technical Field

[0001] This invention relates to the field of NDVI timing prediction technology, and more specifically to an NDVI timing prediction method that combines physical equation constraints with long short-term memory networks. Background Technology

[0002] In fields such as vegetation dynamics monitoring, ecosystem assessment, and climate change research, the Normalized Difference Vegetation Index (NDVI) is a key indicator reflecting vegetation cover and growth status, and is widely used in vegetation growth analysis, land use classification, and ecological quality assessment. NDVI effectively reflects vegetation greenness information, providing crucial data support for large-scale studies of vegetation spatiotemporal changes, and plays an irreplaceable role in understanding the evolution of regional ecological environments.

[0003] However, there are still many challenges in accurately predicting NDVI time series over large-scale areas: on the one hand, different regions have different vegetation types and asynchronous phenological periods, and are affected by a combination of factors such as topography and climate, resulting in strong heterogeneity and complexity in the spatiotemporal changes of NDVI; on the other hand, NDVI time series data are easily affected by cloud cover, sensor errors, etc., leading to data loss or high noise, which hinders time series modeling; in addition, existing methods mostly rely on data-driven statistical or machine learning models, and rarely incorporate the physiological and ecological mechanisms of vegetation growth, making it difficult to accurately depict the intrinsic laws of NDVI changes, and easily leading to problems of decreased accuracy and insufficient generalization ability in long-term time series prediction.

[0004] To address these issues, researchers have attempted to combine multi-source data with deep learning techniques (such as Long Short-Term Memory (LSTM) networks and Convolutional Neural Networks (CNNs)) to conduct NDVI prediction research. However, existing work either focuses on temporal modeling of single pixels, neglecting the spatial dimension's correlation features, or lacks effective integration of physical mechanisms in model design, resulting in deficiencies in the prediction results' ecological interpretability and accuracy stability.

[0005] Therefore, improving the accuracy and reliability of NDVI time series prediction to better serve ecological monitoring and resource management is an urgent problem to be solved by those skilled in the art. Summary of the Invention

[0006] In view of the above problems, this invention is proposed to provide an NDVI time-series prediction method that combines physical equation constraints with long short-term memory networks to overcome or at least partially solve the above problems. It integrates multiple technical means such as multi-source three-dimensional spatiotemporal composite construction, physical-data dual-driven LSTM-PINN network modeling, and constraint condition construction based on vegetation ecological mechanisms. It aims to achieve high-precision prediction of NDVI time-series changes in target areas by deeply integrating physiological and ecological mechanisms such as vegetation photosynthesis and water stress with the inherent patterns of NDVI long-term data.

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

[0008] This invention provides an NDVI time series prediction method that combines physical equation constraints with long short-term memory networks, comprising the following steps: S1. Obtain meteorological variables, vegetation indices and soil data related to NDVI changes in the target area over historical time periods from multiple data sources, perform preprocessing, and obtain a spatiotemporal dataset. S2. Extract the statistical and spatial geometric features of the NDVI time series from the spatiotemporal dataset, and construct an NDVI spatiotemporal cube with dimensions of time-space-environmental factors. S3. Flatten the spatiotemporal cube after seasonal trend decomposition, input it into the trained hybrid multivariate LSTM-PINN model for time series prediction, and obtain the NDVI time series prediction result of the target area for the target time period. The hybrid multivariable LSTM-PINN model consists of cascaded multi-layer bidirectional LSTM units and multi-head self-attention layers, and is trained based on a physical constraint loss function.

[0009] Furthermore, in step S1, The meteorological variables include surface solar radiation (SSR), monthly average temperature (MAT), and monthly cumulative precipitation (MCP). The vegetation indices include Enhanced Vegetation Index (EVI), Leaf Area Index (LAI), and Normalized Difference Water Index (NDWI). The soil data includes soil moisture.

[0010] Furthermore, in step S1, The preprocessing includes geometric correction, projection stitching, resampling, maximum value synthesis, Savitzky-Golay filtering, and normalization.

[0011] Furthermore, step S2 specifically includes: The sliding window method was used to extract the statistical features of NDVI time series in the spatiotemporal dataset, and the spatial geometric features were extracted using the spatial analysis tools of the geographic information system. The statistical features included the trend of NDVI change and the vegetation growth level. The spatial geometric features included the topographic relief, slope aspect and distance from the water body. The interaction detection values ​​between each variable in the spatiotemporal dataset and NDVI are calculated using a geographic detector model to identify key driving variables; among which, the key driving variables include: meteorological variables, vegetation index and soil data; Based on the aforementioned key driving variables, statistical features, and spatial geometric features, an NDVI multidimensional spatiotemporal cube with dimensions of time-space-environmental factors is constructed.

[0012] Furthermore, in step S3, the spacetime cube is decomposed and then flattened, specifically including: The spatiotemporal cube is subjected to seasonal trend decomposition to separate the trend component, seasonal component and residual component, and the decomposed components are flattened into one-dimensional feature vectors.

[0013] Furthermore, in step S3, the data processing procedure for time series prediction using the hybrid multivariate LSTM-PINN model includes: The input one-dimensional feature vector is standardized through a data normalization layer; The standardized sequence data is fed into a multi-layer bidirectional LSTM unit. By utilizing its forward and backward memory mechanisms, the influence of historical states on the current moment and the implicit constraints of future context on the current state in the time series are modeled simultaneously, capturing the dynamic dependency of NDVI and its driving factors over a long time span. The high-dimensional time series representation output by the multi-layer bidirectional LSTM unit is passed to the multi-head self-attention layer. By calculating the correlation weights between each time step and between different environmental variable dimensions, the key time periods and key driving factors that play a dominant role in NDVI prediction are dynamically identified, and selective enhancement and redundancy suppression of the information channel are performed. The feature sequences weighted and fused by the self-attention mechanism are sequentially passed through a Dropout layer and a fully connected layer. The Dropout layer randomly masks some neuron connections to prevent overfitting. The fully connected layer maps high-dimensional features into single-step or multi-step NDVI prediction values, completing the end-to-end transformation from multi-source heterogeneous temporal input to target vegetation index output.

[0014] Furthermore, in step S3, the multi-layer bidirectional LSTM unit includes multiple stacked bidirectional LSTM layers; each bidirectional LSTM layer is composed of a forward LSTM sub-network and a backward LSTM sub-network in parallel, which are used to process forward and backward time series information respectively; each sub-network updates its hidden state and cell state through the LSTM unit at each time step.

[0015] Furthermore, each of the LSTM units includes a cell state, a forget gate, an input gate, and an output gate; The process of updating the hidden state and cell state through LSTM units at each time step specifically includes four interrelated computational steps: 1) Calculation of the forget gate; Based on the current time step t Input feature vector The hidden state of the previous time step Calculate the output of the forget gate :

[0016] in, It is the Sigmoid activation function. and These are the learnable weight matrix and bias term of the forget gate, respectively; 2) Input gate and candidate cell state generation; Update the input gates that are written to the cell state in the current input. :

[0017] Based on the current time step t New memory content added, candidate cell status updated. :

[0018] in, and These are the learnable weight matrix and bias term of the input gate, respectively; and These are the learnable weight matrix and bias term for the candidate cell state, respectively; tanh For activation functions; 3) Cell state renewal; Update the current time step by combining the forget gate, input gate, and candidate cell states. t Cellular state :

[0019] in, This represents the cell state at the previous time step. 4) Generation of output gates and hidden states; Calculate the output gate of the current cell state. :

[0020] calculate t Hidden state of time :

[0021] in, and These are the learnable weight matrix and bias term of the output gate, respectively.

[0022] Furthermore, in step S3, the physical constraint loss function includes: partial differential equation constraint loss L. PDE Initial condition loss L IN and boundary condition loss L BC ; Wherein, the L IN This includes physical range constraints for vegetation NDVI, minimum greenness constraints, and prohibition of negative photosynthesis; the L BC This includes constraints on vegetation growth rate, water stress response, temperature response curves, photosynthetic efficiency, vegetation index consistency, and seasonality of vegetation growth; the L PDE This includes cumulative biomass constraints.

[0023] As can be seen from the above technical solution, compared with the prior art, the present invention discloses an NDVI time series prediction method that combines physical equation constraints with long short-term memory networks, which has the following beneficial effects: This invention achieves in-depth integration and mechanistic modeling of multiple factors affecting NDVI changes, including meteorology, soil, and vegetation physiology, through the construction of a multi-source three-dimensional spatiotemporal composite, the construction of constraints based on vegetation ecological mechanisms, and a physical-data dual-driven LSTM-PINN network.

[0024] On the one hand, the construction of a three-dimensional spatiotemporal cube based on time, space, and environmental factors fully preserves the spatiotemporal dynamic correlation of NDVI. Combined with key covariate screening, it effectively characterizes the complex coupling relationship between NDVI and driving factors, laying the foundation for high-precision prediction. On the other hand, the bidirectional LSTM unit's ability to capture long-term time-series dependencies, combined with the constraints of physical equations embedded in the loss function, avoids the risk of overfitting in data-driven models while ensuring that the prediction results conform to vegetation physiological and ecological mechanisms, significantly improving the accuracy and reliability of NDVI time-series prediction. This invention helps to more accurately monitor regional vegetation growth dynamics, assess ecological environment quality, and provide scientific basis and decision support for agricultural production regulation, ecological restoration planning, and climate change response. Attached Figure Description

[0025] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0026] Figure 1 This is a flowchart of the NDVI time series prediction method combining physical equation constraints and long short-term memory networks provided in this embodiment of the invention. Figure 2 This is a schematic diagram illustrating the construction principle of the NDVI and multi-source driving variable relationship three-dimensional spatiotemporal combination provided in the embodiments of the present invention; Figure 3 This is a schematic diagram of the LSTM-PINN model processing provided in an embodiment of the present invention. Detailed Implementation

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

[0028] This invention discloses an NDVI time series prediction method that combines physical equation constraints with long short-term memory networks, referring to... Figure 1 As shown, it includes the following steps: S1. Obtain meteorological variables, vegetation indices and soil data related to NDVI changes in the target area over historical time periods from multiple data sources, perform preprocessing, and obtain a spatiotemporal dataset. S2. Extract the statistical and spatial geometric features of NDVI time series from the spatiotemporal dataset, and construct an NDVI spatiotemporal cube with dimensions of time-space-environment factors. S3. Flatten the spatiotemporal cube after seasonal trend decomposition, input it into the trained hybrid multivariate LSTM-PINN model for time series prediction, and obtain the NDVI time series prediction results for the target area and target time period. The hybrid multivariate LSTM-PINN model consists of cascaded multi-layer bidirectional LSTM units and multi-head self-attention layers, and is trained based on a physical constraint loss function.

[0029] Grassland ecosystems are highly sensitive to climate change and human activities, urgently requiring high-precision, interpretable vegetation dynamics predictions to support ecological early warning and grassland-livestock balance management. Traditional pure data-driven models have poor generalization ability in extreme drought years, while pure mechanistic models struggle to characterize complex nonlinear responses. This embodiment employs the method of this invention, integrating multi-source remote sensing and meteorological data to achieve highly robust NDVI time-series predictions for a specific grassland region under the constraints of physical laws.

[0030] The implementation process of this embodiment is described in detail below: Step S1: Acquisition and preprocessing of multi-source data.

[0031] This embodiment obtains NDVI changes and related variables such as meteorological, vegetation index, and soil data from multiple data sources, performs preprocessing, and generates a spatiotemporal dataset in a unified format.

[0032] In this embodiment, multiple data sources include, but are not limited to, MODIS (Moderate Resolution Imaging Spectroradiometer), CHIRPS (Climate Hazards Center InfraRed Precipitation with Station data), ERA5-Land (European Centre for Medium-Range Weather Forecasts Reanalysis), and GRIDMET (University of Idaho Gridded Surface Meteorological Dataset).

[0033] The data in this embodiment covers the period from 2003 to 2022, with spatial resolution varying depending on the data source. Vegetation indices include Enhanced Vegetation Index (EVI), Leaf Area Index (LAI), and Normalized Difference Water Index (NDWI); meteorological variables include surface solar radiation (SSR), monthly average temperature (MAT), and monthly cumulative precipitation (MCP); and soil data includes soil moisture.

[0034] The preprocessing operations in this embodiment include geometric correction, projection stitching, resampling, maximum value synthesis, and Savitzky-Golay filtering to reduce the impact of cloud contamination and observation noise. Through preprocessing, the spatial and temporal resolution differences between different data sources are unified into a form suitable for model input.

[0035] Step S2: Construction of a three-dimensional spacetime cube.

[0036] This embodiment extracts features of NDVI and key covariates, including temporal statistical features and spatial geometric features, to construct a spatiotemporal cube.

[0037] The key covariates are meteorological, vegetation index, and soil factors related to NDVI; the feature information includes statistical features in the time series data (such as mean, standard deviation, trend slope, etc.) and geometric features in the spatial distribution (such as adjacency similarity, slope variance, etc.). These feature information are constructed into a spatiotemporal cube.

[0038] This embodiment uses the sliding window method to extract the temporal features of NDVI, and uses GIS spatial analysis tools to extract spatial geometric features such as topographic relief, slope aspect, and distance from water bodies; and uses a geographic detector model to calculate the interaction detection values ​​between each driving variable and NDVI.

[0039] Here, the key covariates are those that have a high impact on NDVI, including statistical characteristics, trend characteristics, periodic characteristics, correlation characteristics, etc. in multi-source time series data.

[0040] This embodiment introduces methods such as geographic detectors to screen key covariates, ensuring that the selected variables can fully explain the spatiotemporal variations of NDVI, thereby reducing redundant information and lowering the complexity and uncertainty of the prediction results.

[0041] Reference Figure 2 As shown, a schematic diagram illustrating the construction principle of the NDVI and multi-source driving variable relationship three-dimensional spatiotemporal composite is presented. Figure 2Multi-source temporal imagery represents remote sensing data sequences from different times, sensors, or inversion products, including vegetation index products such as NDVI and EVI; influencing factors refer to the key driving variables after screening by the geospatial detector, with bidirectional arrows indicating correlation and coupling between different variables; feature alignment refers to spatial registration, temporal synchronization, and variable correspondence of multi-source data sampled at different sources, resolutions, and times at the pixel scale; the three-dimensional cube refers to combining the aligned multi-source variables in the spatial dimension (x, y), temporal dimension (t), and feature dimension to form a three-dimensional spatiotemporal data cube. Where t-1, t, and t+1 represent adjacent moments in a continuous time series, with each time point corresponding to a set of multi-source remote sensing images; in the 3D cube, the horizontal and vertical axes represent the spatial dimension (x, y), and the depth direction represents the time dimension (t); different color layers inside the cube represent different variables or feature channels; in the time series curve, the horizontal axis represents the time step, and the vertical axis represents the normalized or standardized values ​​of NDVI and related driving variables; different colored curves represent the time series variation curves of NDVI and the synchronous variation curves of different driving variables; the blue double-headed arrows represent the correlation and temporal relationship between variables in the time dimension.

[0042] Step S3: Construct a hybrid multivariate LSTM-PINN model for time series prediction.

[0043] The hybrid multivariate LSTM-PINN model in this embodiment cascades multiple bidirectional LSTM units and multi-head self-attention layers, and combines them with a physical constraint loss function. It can characterize long-term temporal dependence and vegetation physiological and ecological mechanisms, and achieve spatiotemporal prediction of NDVI by using multiple environmental factors as joint inputs.

[0044] In this embodiment, a hybrid LSTM-PINN model is constructed using the Keras neural network API in Python: First, two LSTM units are stacked to form a bidirectional LSTM layer; then, an encoder with positional encoding, scaled dot product multi-head self-attention (including residual connections and layer normalization) is inserted, followed by a fully connected layer; finally, the original loss function is replaced by a physical constraint equation. The entire process is accelerated using GPU training. This model, with its cascaded structure of "LSTM temporal memory + PINN physical modeling," gradually integrates the influence of multiple factors on the target variable, providing more accurate and comprehensive prediction results, and can automatically extract key features from the input data during training. Specific information for each layer is as follows: The input layer receives NDVI and multi-factor feature information extracted from time-series datasets from multiple sources. This feature information is constructed into a spatiotemporal cube, flattened, and used as input to the fusion neural network.

[0045] Each LSTM unit contains a cell state, a forget gate, an input gate, and an output gate. This gating mechanism controls information flow, preventing gradient explosion and vanishing problems, and effectively capturing long-term dependencies in the data. Compared to single-layer LSTMs, multi-layer structures can better learn abstract representations of data and capture feature information at different levels. Each LSTM layer not only processes temporal data of NDVI but also simultaneously processes multiple variables affecting NDVI, enhancing the model's multi-factor feature extraction capabilities.

[0046] The multi-head self-attention layer following LSTM embeds spatiotemporal location information through position encoding and uses multiple sets of query-key-value (QKV) mappings to calculate scaled dot product attention, thereby capturing long-distance dependencies between different time steps and spatial units globally. Residual connections and layer normalization ensure gradient stability and accelerate convergence. The subsequent feedforward network further refines the feature representation, which helps to integrate the interactive effects of multiple variables on NDVI.

[0047] The output layer constructs a loss function using physical constraint equations, replacing the traditional loss function. This loss function incorporates the partial differential equation constraint loss L. PDE Initial condition loss L IN and boundary condition loss L BC The model's prediction results are constrained from three dimensions: dynamic evolution patterns, initial state, and boundary behavior, ensuring that the predictions conform to both the temporal patterns of the data and the physical mechanisms of vegetation physiology and ecology. Then, the hidden states, iteratively optimized by the hybrid network, are received and transformed into predicted values ​​through a fully connected layer. The predicted values ​​include not only the NDVI values ​​but also their associated spatial distribution information.

[0048] In this embodiment, each LSTM unit contains a cell state. The three neural network layers are the forget gate. Input gate and output gate The calculation process of the LSTM unit at time t is as follows: 1. Forget Gate Calculation Based on the current time step t Input feature vector (Normalized multivariate fusion features) and the hidden state of the previous time step Calculate the output of the forget gate :

[0049] in, It is the Sigmoid activation function. and These are the learnable weight matrix and bias term of the forget gate, respectively; the gate signal. The decision is based on the cell state at the previous time step. How much historical information is retained, such as early vegetation growth trends or water accumulation effects.

[0050] 2. Input gate and candidate cell state generation The following two computations are performed in parallel: Input gate :

[0051] Controls how much new information from the current input will be written to the cell state; Candidate cell status:

[0052] Generate new memory content that may be added at the current time step, reflecting the immediate impact of factors such as monthly precipitation and radiation on vegetation. and These are the learnable weight matrix and bias term of the input gate, respectively; and , , represent the learnable weight matrix and bias term for the candidate cell state, respectively; tanh is the activation function.

[0053] 3. Cell state renewal By combining the forget gate, input gate, and candidate states, the cell state at the current time step is updated. :

[0054] in, This represents the cell state at the previous time step. 4. Generation of output gates and hidden states Output gate :

[0055] Determine which parts of the current cell state will be exported; Hidden state :

[0056] in, and These are the learnable weight matrix and bias term of the output gate, respectively. This represents the dot product operation.

[0057] In the feedforward LSTM subnetwork, time step t traverses the input sequence in natural order (t=1,2,…,T) to capture the causal dependencies from the past to the present. In the backward LSTM subnetwork, time steps t are processed in reverse order (t=T,T-1,…,1) to capture the implicit constraints of future context on the current state (such as the symmetry of seasonal cycles). The hidden states of the two sub-networks are output after each step of concatenation, forming bidirectional context-aware temporal features, which provide a high-dimensional semantic foundation for subsequent physical constraint embedding and attention mechanisms.

[0058] Reference Figure 3 As shown, Figure 3 The input layer contains the multi-dimensional spatiotemporal cube, where the horizontal and vertical axes represent the spatial dimension (x, y), and the depth direction represents the time dimension (t). After STL decomposition, the trend component represents long-term variation characteristics, the seasonal component represents periodic variation patterns, and the residual component represents short-term disturbances and abnormal fluctuations. Normalization processing involves normalizing the decomposed components to ensure that the data distribution of different variables and components is within a uniform numerical scale. The LSTM time series modeling layer includes multiple time-step LSTM units, each of which sequentially includes: a forget gate, an input gate, an output gate, and a memory state update structure. Represents the hidden state vector. Indicates the state of the memory cell. and Representing time steps +1 hidden state and memory state Here, is the Sigmoid activation function, and tanh is the activation function. + represents element-wise multiplication, and + represents element-wise addition; each of the eight attention heads in the self-attention layer performs a linear mapping on the input features, generating a query vector (Q), a key vector (K), and a value vector (V); PINN is the physical constraint loss function, where the hidden layer uses a neural network structure to process the physical state variables. , , , Mapping and fitting are performed, and automatic micro-leveling is used to automatically solve for the derivative of the network output with respect to time or physical variables, where... For the automatic differentiating operator, the total loss function includes the L... PDE L IN L BC L PDE This refers to the constraint loss of partial differential equations, L IN This refers to the initial condition loss, L BC This refers to boundary condition loss, L PDEIncluding biomass accumulation constraints, L IN Including vegetation NDVI range (0,1), minimum greenness constraint, prohibition of negative photosynthesis, L BC Including vegetation growth rate limitation, water stress response, temperature response curve, photosynthetic efficiency constraint, multi-index consistency, and vegetation growth seasonality constraint; total loss value. The weighted sum of the mean squared error loss function and the preset threshold Compare when less than or equal to The output predicted sequence is the final result. When it is greater than 100, the output predicted sequence is the final result. The network parameters are adjusted based on feedback, and training and optimization continue.

[0059] This embodiment first extracts the feature information of NDVI and its highly important covariates. Then, an LSTM-PINN model is constructed, and the extracted feature information is used as input for spatiotemporal prediction. Next, a physical constraint loss function is constructed based on vegetation growth mechanisms such as photosynthesis and physical constraints under water stress, and the prediction results are iteratively optimized. Finally, a temporal prediction method for NDVI coupled with physical equation constraints and a long short-term memory network is realized.

[0060] Regarding the physical constraint loss function, this embodiment relies on physiological and ecological processes such as vegetation photosynthesis and water stress to transform the basic equation into a partial differential equation, constructing partial differential equation constraint loss, initial condition loss, and boundary condition loss; and integrates the three to obtain a comprehensive physical constraint loss function.

[0061] In this embodiment, the physical constraint equation (PINN) is transformed into a partial differential equation form based on fundamental equations of vegetation physiological and ecological processes, such as the mechanism of vegetation photosynthesis and physical constraints of water stress. This is used to construct the partial differential equation constraint loss L. PDE The initial conditional loss L is constructed based on the reasonable NDVI range of vegetation and the NDVI value at the initial time. IN ;Analyze the characteristics of vegetation growth rate limitation and seasonal constraints to construct boundary condition loss L BC ; Fusion L PDE L IN L BC We construct a comprehensive physical constraint loss function to replace the traditional loss function based solely on data fitting.

[0062] In this embodiment, L IN Includes vegetation NDVI range (0,1), minimum greenness constraint, and prohibition of negative photosynthesis; L BC This includes vegetation growth rate constraints, water stress response, temperature response curves, photosynthetic efficiency constraints, multi-index consistency, and seasonal constraints on vegetation growth; L PDE This includes cumulative biomass constraints. The calculation is based on the following formula: Vegetation NDVI physical range constraints It can be expressed by the formula:

[0063] in, The total number of samples for NDVI sequences. Indicates the sample index. Indicates the first The normalized difference vegetation index (NDVI) corresponds to each sample, where 1.0 represents the physically reasonable upper limit threshold of NDVI. This represents the corrected linear unit function. Vegetation growth rate limitation It can be expressed by the formula:

[0064] in, The length of the time series. Indicates a time index. Represents the normalized vegetation index at time t. This represents the normalized vegetation index at time t+1. This represents the absolute change in NDVI between adjacent time steps. This represents the maximum reasonable threshold for NDVI variation between adjacent time steps.

[0065] Water stress response It can be expressed by the formula:

[0066]

[0067]

[0068] in, The length of the time series. Indicates a time index. express Comprehensive moisture status index at any given time. express t Rainfall at any given moment express t Soil moisture at any time, This represents the nonlinear response function of water conditions to vegetation growth. Indicates time The overall moisture status at any given time. For time lag terms, This represents the sensitivity adjustment parameter of the moisture response function.

[0069] Temperature response curve It can be expressed by the formula:

[0070]

[0071] in, The length of the time series. Indicates a time index. Indicates the vegetation's response to ambient temperature The response function value, Represents the ambient temperature variable. Indicates the optimal temperature for vegetation growth. The standard deviation of the temperature response function is represented.

[0072] Photosynthetic efficiency constraints It can be expressed by the formula:

[0073]

[0074] in, This represents the light energy utilization efficiency function of vegetation under available light conditions. Indicates variables related to light intensity. This represents the half-saturation constant of the light response function.

[0075] Consistency of vegetation index This includes: EVI constraints and LAI constraints.

[0076] EVI constraints are expressed by the following formula:

[0077] The LAI constraint is expressed by the formula:

[0078]

[0079]

[0080] in, The length of the time series. This represents the time index, and 0.95 represents the empirical saturation limit of NDVI. This represents the expected value of NDVI derived from LAI. This represents the NDVI and EVI consistency term, i.e., the EVI constraint; 2.0 represents the deviation loss term between NDVI and its expected value derived based on LAI, i.e., the LAI constraint; 2.0 represents the weighting coefficients for the EVI constraint and the LAI constraint.

[0081] Seasonal constraints on vegetation growth It can be expressed by the formula:

[0082]

[0083]

[0084] in, Indicates a time index. This represents the annual average NDVI value for the first year. This represents the arithmetic mean of 12 time steps over a year. This represents the annual average NDVI value for the second year. This represents the maximum reasonable threshold for interannual variation in NDVI.

[0085] Biomass accumulation constraint It can be expressed by the formula:

[0086]

[0087]

[0088] in, The length of the time series. Indicates a time index. This represents the state variable of vegetation biomass, corresponding to NDVI or its nonlinear mapping. This represents the intrinsic growth rate of vegetation. Indicates the environmental carrying capacity for vegetation growth. It represents the rate of change of vegetation biomass per unit time. Indicates at time step arrive Biomass increment between This represents the vegetation biomass state at time t. This represents the vegetation biomass state at time t+1. This represents the theoretically expected increase in biomass calculated based on the logistic growth model.

[0089] Prohibition of negative photosynthesis It can be expressed by the formula:

[0090] in, The length of the time series. This represents the time index, and 0.25 represents the maximum reasonable drop threshold allowed for NDVI within a single time step.

[0091] Minimum greenness constraint It can be expressed by the formula:

[0092] in, The length of the time series. Indicates a time index.

[0093] above This represents the physical constraint loss term constructed in this instance.

[0094] In this embodiment, during the construction of the hybrid LSTM-PINN model, L... PDE L IN L BC Construct a comprehensive physical constraint loss function. The constraints of the comprehensive physical constraint function are: vegetation NDVI range constraint, minimum greenness constraint, prohibition of negative photosynthesis constraint, vegetation growth rate limitation constraint, water stress response constraint, temperature response curve constraint, photosynthetic efficiency constraint, vegetation index consistency constraint, vegetation growth seasonality constraint, and biomass accumulation constraint.

[0095] Comprehensive physical constraint loss function Together with the mean squared error loss function, they constitute the overall model loss function. L It can be expressed by the formula:

[0096]

[0097]

[0098] In the formula, This represents the weights corresponding to different physical constraint loss terms. ; This represents the mean squared error loss function. This represents the weight corresponding to the mean squared error loss. Represents the total number of samples. Represents the actual observed value. This represents the predicted value output by the model.

[0099] The ten physical constraints in this embodiment are: vegetation NDVI range constraint, minimum greenness constraint, prohibition of negative photosynthesis constraint, vegetation growth rate limitation constraint, water stress response constraint, temperature response curve constraint, photosynthetic efficiency constraint, multi-index consistency constraint, vegetation growth seasonality constraint, and biomass accumulation constraint. These ten physical constraints are then fused with the mean squared error loss function using a weighted fusion method. Detailed weighting is shown in Table 1 below. Table 1. Weight combination table of physical equation constraints in the hybrid LSTM-PINN neural network model used in this invention.

[0100] This embodiment uses preprocessed data to train a hybrid LSTM-PINN model, and uses GridSearchCV in Scikit-learn to perform parameter tuning to improve prediction performance. It specifies the range of candidate hyperparameter values ​​and iterates through the combinations of these parameters to achieve the parameter tuning process.

[0101] During model training, the input sequence is divided into training and test sets. GridSearchCV is used to traverse candidate hyperparameter combinations to fine-tune parameters such as the number of hidden units, the number of layers, the learning rate, and the Dropout ratio to improve prediction performance.

[0102] This embodiment uses the following evaluation metrics to assess the performance of the hybrid LSTM-PINN model: coefficient of determination R², mean squared error (MSE), root mean square error (RMSE), mean absolute percentage error (MAPE), mean absolute error (MAE), and mean directional accuracy (MDA) to evaluate the accuracy and reliability of the model's prediction results. The model's stability and generalization ability are verified by calculating metrics such as the Population Stability Index.

[0103] In this embodiment, the model training process includes the following steps in sequence: data preparation, temporal feature extraction, model training, and prediction.

[0104] 1. Data Preparation 1.1 Data Sources: To perform NDVI time-series forecasting, various relevant data need to be collected, including NDVI data, meteorological data (such as surface solar radiation SSR, monthly mean temperature MAT, monthly cumulative precipitation MCP, etc.), vegetation index data (such as leaf area index LAI, enhanced vegetation index EVI, normalized difference water index NDWI, etc.), and soil data (such as topsoil moisture content SSM). Data sources mainly include MODIS, CHIRPS, ERA5-Land, and GRIDMET, covering the period from 2003 to 2022, with spatial resolution varying depending on the data source.

[0105] 1.2 Data Preprocessing: In the data preprocessing stage, interpolation and resampling operations are required to ensure data integrity and consistency. In particular, it is necessary to address the differences in spatial and temporal resolution between different data sources and unify them into a format suitable for model input. For missing values, methods such as nearest-neighbor interpolation or time-series interpolation are used for imputation. Savitzky-Golay filtering is then applied to the images to eliminate the influence of noise.

[0106] 2. Temporal Feature Extraction On the constructed time-series dataset, the system systematically extracts the feature information of NDVI (Normalized Difference Vegetation Index) and its highly important covariates (including meteorological data, vegetation indices, soil properties, etc.). Among them, the time-series dimension features cover statistical features such as the mean (reflecting the long-term average level), standard deviation (reflecting the degree of data fluctuation), maximum and minimum values ​​(characterizing the extreme range of values), and trend slope (revealing the time-series change trend); the spatial dimension features include geometric features such as the adjacency similarity of spatial distribution (reflecting the correlation between adjacent units) and slope variance (reflecting the differences in topographic relief), realizing a comprehensive characterization of the data in both the temporal and spatial dimensions.

[0107] By integrating spatiotemporal dual-dimensional feature information, a corresponding spatiotemporal cube is constructed, which fully preserves the spatiotemporal correlation attributes of NDVI and covariates. Subsequently, seasonal trend decomposition is performed on the spatiotemporal cube sequence to separate the trend term, seasonal term, and residual term, further exploring the temporal patterns. Finally, the decomposed spatiotemporal cube is flattened into a fixed-dimensional feature vector to ensure that the input data format meets the training requirements of the hybrid LSTM-PINN network, providing structurally complete and information-rich input data support for subsequent prediction tasks.

[0108] 3. Model Training 3.1 Physical constraint equations The physical constraint loss function in this embodiment includes: partial differential equation constraint loss L PDE Initial condition loss L IN and boundary condition loss L BC Among them, L IN Includes vegetation NDVI range (0,1), minimum greenness constraint, and prohibition of negative photosynthesis; L BC This includes vegetation growth rate constraints, water stress response, temperature response curves, photosynthetic efficiency constraints, multi-index consistency, and seasonal constraints on vegetation growth; L PDE This includes cumulative biomass constraints.

[0109] This embodiment focuses on the spatiotemporal dynamics of vegetation biomass, treating vegetation growth as a constrained dynamic process that follows the law of conservation of matter to construct an LPDE. Its significance lies in transforming the differential equation into residual loss at discrete time steps, thereby forcing the biomass accumulation simulated by the model to conform to the law of conservation of matter and avoiding the phenomenon of sudden reduction in biomass without reason and violating physical logic.

[0110] The NDVI range constraint term, based on the physical definition of NDVI for surface vegetation, forces the NDVI value output by the model to be limited to the [0,1] interval, excluding outliers exceeding the physical definition. The minimum greenness constraint term sets a lower threshold for NDVI in vegetated areas to characterize the basic greenness level that vegetation needs to maintain under stress conditions, avoiding the contradictory state of having vegetation but no greenness in the simulation. The prohibition of negative photosynthesis constraint term stipulates that the instantaneous decrease in NDVI must not exceed a physiologically permissible threshold, which is set to 0.25 in this study. Biomass accumulation is a slow process, and photosynthesis is the cornerstone of energy metabolism. The net photosynthetic rate is usually non-negative or changes slowly. This loss term, by penalizing excessively rapid NDVI decreases, indirectly forces the model's state changes between time steps t and t+1 to conform to the physiological inertia of vegetation material accumulation and dissipation, avoiding non-physical instantaneous collapse of biomass.

[0111] The vegetation growth rate constraint aims to reflect the physiological inertia and maximum potential limit of vegetation growth. Even under optimal environmental conditions, these fundamental cellular and tissue-level processes have a physical limit. This constraint prevents the model from simulating instantaneous collapse or overgrowth of biomass or greenness within short time steps. The water stress response constraint forces NDVI predictions to be consistent with the expected greenness under water constraints. The temperature response curve constraint forces NDVI predictions to be consistent with the limiting effect of temperature on vegetation growth. The photosynthetic efficiency constraint sets a theoretical upper limit on light energy use efficiency (LUE), preventing the model output efficiency from exceeding biological limits. This loss term indirectly constrains the logical relationship between NDVI and LUE, preventing the model from simulating the non-physical phenomenon of extremely high NDVI accumulation under low light energy use efficiency. The multi-index consistency constraint mainly forces NDVI to maintain synergistic consistency with other remote sensing indices such as EVI and LAI in terms of trends and values, thereby enhancing the physical credibility of the simulation results across multiple remote sensing dimensions. The vegetation growth seasonality constraint forces the model to adhere to the macroscopic inertia of vegetation under climatic background. It is crucial for enhancing the generalization ability of models and the reliability of long-term predictions, effectively preventing non-physical interannual drift or extreme instability in multi-decadal simulations.

[0112] 3.1 Hybrid LSTM-PINN Model This embodiment introduces an LSTM-PINN hybrid network. Building upon the strength of LSTM in capturing long-term temporal dependencies and its multi-head self-attention mechanism, a physical equation loss constraint function is superimposed to learn the ecological mechanisms of vegetation growth. Compared to a single LSTM or traditional machine learning models, this structure avoids both gradient explosion and overfitting caused by purely data-driven approaches, thus exhibiting higher accuracy and generalization ability in NDVI prediction.

[0113] This embodiment uses Keras, a neural network API in Python, to build a hybrid LSTM-PINN model. First, two layers of bidirectional LSTM units are stacked to form a deep temporal encoder. Input and output layers are then constructed, and GPU acceleration is used for model training. By cascading multiple layers of bidirectional LSTM units and a multi-head self-attention (Transformer Encoder) layer, the influence of multiple factors on the target NDVI is gradually fused, resulting in more accurate and comprehensive predictions. Dropout and Dense layers are added to suppress overfitting and achieve non-linear feature mapping. A fully connected layer is added before the output layer, using a weighted fusion loss function based on the constructed NDVI-related physical equations instead of the traditional deep learning MSE loss. This hybrid model leverages the gated memory of LSTM and the physical mechanism of PINN to automatically extract temporal features from the input data. Simultaneously, its parallel structure can process multiple input sequences and share weights among the physical mechanisms corresponding to different variables, learning both the patterns of multiple factors affecting NDVI and preserving the unique physiological and ecological processes of each variable.

[0114] 3.2 Model Training and Validation: The input-output flow of the hybrid LSTM-PINN model provided in this embodiment is as follows: Figure 3 As shown, before training begins, each input sequence is divided into two segments chronologically: the first 80% is used as the training set, and the latter 20% is used as the test set. Then, GridSearchCV is used to perform a grid search on the model's hyperparameters for each feature vector to improve prediction performance. This study pre-sets candidate values ​​for learning rate, number of hidden units, Dropout ratio, batch size, and physical loss weights, traversing all combinations and recording validation errors to determine the optimal parameters. All variables share the same physical constraint loss function, simultaneously optimizing the prediction results for all variables. Furthermore, 3x cross-validation (CV) is introduced during training to verify the reliability of the optimal parameters.

[0115] 4. Prediction: Input the test set data into the trained hybrid LSTM-PINN model, and iteratively validate and optimize to obtain the NDVI prediction value for each pixel. These prediction values ​​reflect the temporal variation characteristics of NDVI under given multi-factor constraints.

[0116] 5. Performance Evaluation: To verify the performance of the LSTM-PINN method, various evaluation metrics were used to assess the prediction results, including the coefficient of determination R², mean squared error (MSE), root mean square error (RMSE), mean absolute percentage error (MAPE), mean absolute error (MAE), and mean directional accuracy (MDA). Furthermore, this embodiment compares the proposed method with other mainstream prediction methods (such as GWconvLSTM, Transformer, and CNN) to demonstrate the superiority of the invention. The model execution results are compared in Table 2 below: Table 2. Evaluation Comparison of Different Models on the Validation Set

[0117] The proposed NDVI time-series prediction method, which couples physical equation constraints with a long short-term memory network, significantly improves the accuracy and generalization ability of NDVI time-series prediction by integrating multiple techniques, including the construction of a multi-source three-dimensional spatiotemporal composite, LSTM-PINN network modeling driven by both physical and data, and the construction of constraints based on vegetation ecological mechanisms. This method demonstrates excellent performance in validation applications, verifying its broad application prospects and promotional value in vegetation dynamic monitoring and ecosystem stability assessment. Future work will continue to optimize the network structure, incorporate arable land factors, and expand to national and even global scale scenarios, providing more reliable technical support for national arable land vegetation dynamic monitoring and disaster emergency response.

[0118] Compared to existing technologies, the embodiments of the present invention have at least the following advantages: Improved prediction accuracy: By integrating multiple covariates of meteorology, vegetation index, and soil moisture, and fully coupling the potential spatiotemporal relationships between NDVI and environmental factors, prediction errors are significantly reduced. Prediction results conform to ecological laws: By constructing L... PDE L IN L BC The three physical losses are combined with the LSTN-PINN model for iterative verification and optimization of prediction results, making the predictions more consistent with vegetation growth mechanisms and greatly improving the model's generalization ability. Key technical support is provided: high-precision, automated spatiotemporal prediction solutions are offered for regional vegetation dynamic monitoring, ecosystem stability assessment, and climate change response analysis.

[0119] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0120] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for predicting NDVI time series by combining physical equation constraints with long short-term memory networks, characterized in that, Includes the following steps: S1. Obtain meteorological variables, vegetation indices and soil data related to NDVI changes in the target area over historical time periods from multiple data sources, perform preprocessing, and obtain a spatiotemporal dataset. S2. Extract the statistical and spatial geometric features of the NDVI time series from the spatiotemporal dataset, and construct an NDVI spatiotemporal cube with the dimensions of time-space-environmental factors. S3. Flatten the spatiotemporal cube after seasonal trend decomposition, input it into the trained hybrid multivariate LSTM-PINN model for time series prediction, and obtain the NDVI time series prediction result of the target area for the target time period. The hybrid multivariable LSTM-PINN model consists of cascaded multi-layer bidirectional LSTM units and multi-head self-attention layers, and is trained based on a physical constraint loss function.

2. The method as described in claim 1, characterized in that, In step S1, The meteorological variables include surface solar radiation (SSR), monthly average temperature (MAT), and monthly cumulative precipitation (MCP). The vegetation indices include Enhanced Vegetation Index (EVI), Leaf Area Index (LAI), and Normalized Difference Water Index (NDWI). The soil data includes soil moisture.

3. The method as described in claim 1, characterized in that, In step S1, The preprocessing includes geometric correction, projection stitching, resampling, maximum value synthesis, Savitzky-Golay filtering, and normalization.

4. The method as described in claim 1, characterized in that, Step S2 specifically includes: The sliding window method was used to extract the statistical features of NDVI time series in the spatiotemporal dataset, and the spatial geometric features were extracted using the spatial analysis tools of the geographic information system. The statistical features included the trend of NDVI change and the vegetation growth level. The spatial geometric features included the topographic relief, slope aspect and distance from the water body. The interaction detection values ​​between each variable in the spatiotemporal dataset and NDVI are calculated using a geographic detector model to identify key driving variables; among which, the key driving variables include: meteorological variables, vegetation index and soil data; Based on the aforementioned key driving variables, statistical features, and spatial geometric features, an NDVI multidimensional spatiotemporal cube with dimensions of time-space-environmental factors is constructed.

5. The method as described in claim 1, characterized in that, In step S3, the spacetime cube is decomposed and then flattened, specifically including: The spatiotemporal cube is subjected to seasonal trend decomposition to separate the trend component, seasonal component and residual component, and the decomposed components are flattened into one-dimensional feature vectors.

6. The method as described in claim 5, characterized in that, In step S3, the data processing procedure for time series prediction by the hybrid multivariate LSTM-PINN model includes: The input one-dimensional feature vector is standardized through a data normalization layer; The standardized sequence data is fed into a multi-layer bidirectional LSTM unit. By utilizing its forward and backward memory mechanisms, the influence of historical states on the current moment and the implicit constraints of future context on the current state in the time series are modeled simultaneously, capturing the dynamic dependency of NDVI and its driving factors over a long time span. The high-dimensional time series representation output by the multi-layer bidirectional LSTM unit is passed to the multi-head self-attention layer. By calculating the correlation weights between each time step and between different environmental variable dimensions, the key time periods and key driving factors that play a dominant role in NDVI prediction are dynamically identified, and selective enhancement and redundancy suppression of the information channel are performed. The feature sequences weighted and fused by the self-attention mechanism are sequentially passed through a Dropout layer and a fully connected layer. The Dropout layer randomly masks some neuron connections to prevent overfitting. The fully connected layer maps high-dimensional features into single-step or multi-step NDVI prediction values, completing the end-to-end transformation from multi-source heterogeneous temporal input to target vegetation index output.

7. The method as described in claim 6, characterized in that, In step S3, the multi-layer bidirectional LSTM unit includes multiple stacked bidirectional LSTM layers; each bidirectional LSTM layer is composed of a forward LSTM sub-network and a backward LSTM sub-network in parallel, which are used to process forward and backward time series information, respectively. Each subnetwork updates its hidden state and cell state at each time step via LSTM units.

8. The method as described in claim 7, characterized in that, Each LSTM unit includes a cell state, a forget gate, an input gate, and an output gate; The process of updating the hidden state and cell state through LSTM units at each time step specifically includes four interrelated computational steps: 1) Forget gate calculation; Based on the current time step t Input feature vector The hidden state of the previous time step Calculate the output of the forget gate : in, It is the Sigmoid activation function. and These are the learnable weight matrix and bias term of the forget gate, respectively; 2) Input gate and candidate cell state generation; Update the input gates that are written to the cell state in the current input. : Based on the current time step t New memory content added, candidate cell status updated. : in, and These are the learnable weight matrix and bias term of the input gate, respectively; and These are the learnable weight matrix and bias term for the candidate cell state, respectively; tanh For activation functions; 3) Cell state renewal; Update the current time step by combining the forget gate, input gate, and candidate cell states. t Cellular state : in, This represents the cell state at the previous time step. 4) Generation of output gates and hidden states; Calculate the output gate of the current cell state. : calculate t Hidden state of time : in, and These are the learnable weight matrix and bias term of the output gate, respectively.

9. The method as described in claim 8, characterized in that, In step S3, the physical constraint loss function includes: partial differential equation constraint loss L. PDE Initial condition loss L IN and boundary condition loss L BC ; Wherein, the L IN This includes physical range constraints for vegetation NDVI, minimum greenness constraints, and prohibition of negative photosynthesis; the L BC This includes constraints on vegetation growth rate, water stress response, temperature response curves, photosynthetic efficiency, vegetation index consistency, and seasonality of vegetation growth; the L PDE This includes cumulative biomass constraints.

Images