Method for generating high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning

Abstract

The present invention discloses a method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning, including: establishing an XGBoost regression model with Landsat surface reflectance and solar zenith angle as inputs and GLASS 250-meter NDVI as the output to estimate the NDVI under Landsat 30-meter cloud-free pixels; for Landsat pixels contaminated by clouds, finding similar pixels within the spatio-temporal window by constructing a seasonal curve, and filling them by weighting the similar pixels and the GLASS NDVI value to generate a temporally continuous Landsat NDVI curve; using the EEMD method to perform temporal decomposition on the Landsat NDVI curve to obtain the intrinsic mode function and the residual curve; constructing a Bi TCN model to perform temporal prediction on the intrinsic mode function and the residual curve, and combining the prediction results to obtain the inversion result of the high spatio-temporal resolution NDVI; it can accurately estimate the NDVI time series with a resolution of 30 meters, providing important data support for applications such as vegetation growth monitoring, crop yield estimation, and ecosystem change assessment.

Description

Technical Field

[0001] The present invention relates to the technical field of quantitative remote sensing parameter inversion of vegetation, in particular to a method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning. Background Art

[0002] The Normalized Difference Vegetation Index (NDVI) is an important parameter for characterizing the surface vegetation cover and growth status, and has wide applications in the fields of crop monitoring, ecological environment assessment, climate change research, etc. At present, obtaining NDVI based on satellite remote sensing data has become the main means for obtaining NDVI at the regional scale or even the global scale. However, there are still many limitations in the spatio-temporal resolution and accuracy of existing NDVI products. Existing NDVI products mainly include two categories: one is the products generated based on medium and low spatial resolution satellite data such as MODIS, such as the GLASS NDVI product with a resolution of 250 meters, the MOD13 product with a resolution of 500 meters, etc. Although such products have a high temporal resolution and good temporal continuity, their rough spatial resolution is difficult to depict the fine spatial heterogeneity of surface vegetation, restricting their application in fine-scale vegetation monitoring. The other is the products generated based on high spatial resolution satellite data such as Landsat. Although the spatial resolution of 30 meters can better reflect the spatial details of surface vegetation, due to the influence of factors such as clouds and fog, the effective observation data is sparse and the temporal continuity is poor, making it difficult to meet the requirements of dynamic monitoring.

[0003] To improve the spatio-temporal resolution of NDVI products, researchers have proposed various methods: (1) Methods based on data fusion, such as algorithms like STARFM and ESTARFM, generate NDVI products with relatively high spatio-temporal resolution by fusing remote sensing data with high temporal resolution and high spatial resolution. However, such methods require obtaining high-quality reference images and have relatively low prediction accuracy during periods of drastic surface changes; (2) Methods based on time interpolation, such as Whittaker smoother and Savitzky-Golay filtering, reconstruct continuous NDVI time series by interpolating discrete NDVI observations. However, such methods often neglect the spatial correlation of NDVI and have poor reconstruction effects during periods with sparse observations; (3) Methods based on machine learning, such as random forest and deep learning, predict by establishing relationships between NDVI and other remote sensing features. However, such methods often require a large number of high-quality training samples, and it is difficult to ensure the temporal continuity of prediction results. In addition, NDVI time series have obvious multi-scale characteristics, including information on multiple time scales such as seasonal changes, interannual fluctuations, and long-term trends. Existing methods often model the NDVI time series as a whole, making it difficult to fully utilize the change characteristics of different time scales and affecting the accuracy of prediction results.

[0004] Therefore, there is an urgent need to develop new methods that can not only make full use of the advantages of multi-source remote sensing data but also effectively extract and utilize the multi-scale characteristics of NDVI time series to generate NDVI products with high spatio-temporal resolution and good temporal continuity, providing data support for refined vegetation monitoring. Summary of the Invention

[0005] To solve the problems existing in the prior art, the object of the present invention is to provide a method for generating a normalized difference vegetation index with high spatio-temporal resolution based on modal decomposition and deep learning. The present invention can accurately estimate the NDVI time series at a resolution of 30 meters, providing important data support for applications such as vegetation growth monitoring, crop yield estimation, and ecosystem change assessment.

[0006] To achieve the above object, the technical solution adopted by the present invention is: A method for generating a normalized difference vegetation index with high spatio-temporal resolution, comprising the following steps:

[0007] S1: Establish an XGBoost regression model with Landsat surface reflectance and solar zenith angle as inputs and GLASS 250-meter NDVI as the output to estimate the NDVI under Landsat 30-meter clear-sky pixels;

[0008] S2: For Landsat pixels contaminated by clouds, similar pixels are searched within a spatio-temporal window by constructing seasonal curves, and the weighted similar pixels and GLASS NDVI values are used for filling to generate a temporally continuous Landsat NDVI curve;

[0009] S3: The EEMD method is used to decompose the Landsat NDVI curve in time series to obtain the intrinsic mode functions and the residual curve;

[0010] S4: A Bi TCN model is constructed to perform time series prediction on the intrinsic mode functions and the residual curve, and the combined prediction results are used to obtain the NDVI inversion result with high spatio-temporal resolution.

[0011] As a further improvement of the present invention, in step S1, when selecting representative samples for training the XGBoost regression model, the Kmeans method is used to cluster the time series GLASS NDVI, and sampling points are uniformly selected globally according to the clustering results, and at the same time, it is ensured that all 30-meter pixels under a 250-meter GLASS pixel belong to the same class to weaken the scale effect.

[0012] As a further improvement of the present invention, in step S2, a 3km×3km spatial window centered on the target pixel is established, and for all pixels within the spatial window, the mean Landsat NDVI of the previous and next 5 years is calculated to obtain the seasonal NDVI curve, and the correlation coefficient between the seasonal NDVI curves of other pixels and the target pixel is calculated, and the top 50 pixels with the highest correlation coefficient are selected as similar pixels.

[0013] As a further improvement of the present invention, step S3 is specifically as follows:

[0014] The noise amplitude during EEMD decomposition is set to 0.2 times the standard deviation of the time series. The obtained IMF components are sorted in descending order of frequency to obtain the intrinsic mode functions, and the remaining trend term is used as the residual curve.

[0015] As a further improvement of the present invention, in step S4, Bi TCN models for predicting the intrinsic mode functions and the residual curve are respectively constructed. The input time series length is 5 years, and the output time series length is 1 year;

[0016] The convolutional kernel size of the time convolutional network is set to 3, the dilation rate is [1, 2, 4, 8], the hidden layer dimension is 64, and the hidden layer dimension of the bidirectional gated recurrent unit is set to 128; the root mean square error is used as the loss function, and the Adam optimizer is used for model optimization, and the initial learning rate is set to 0.001; when the validation set loss does not decrease for 5 consecutive epochs, the learning rate is reduced to 0.1 times the original, and when the validation set loss does not decrease for 10 consecutive epochs, the training is stopped.

[0017] As a further improvement of the present invention, the training and application process of the Bi TCN model further includes:

[0018] Dividing the training samples into three parts, that is, 70% of the sample set is used to train the Bi TCN model, 20% of the sample set is used for model verification, and 10% of the sample set is used for independent testing of the model;

[0019] Using the Bi TCN model to estimate the intrinsic mode function and residual curve of the time series, and summing all the intrinsic mode functions and residual curves to obtain the time series NDVI estimation result.

[0020] The beneficial effects of the present invention are:

[0021] 1. The present invention integrates time series modal decomposition and machine learning technologies, makes full use of the advantages of multi-source remote sensing data, and effectively solves the problems of low accuracy and poor time continuity in the traditional method for generating high spatio-temporal resolution NDVI products.

[0022] 2. The present invention selects representative training samples through Kmeans clustering and adopts a spatio-temporal filling strategy based on similarity, effectively weakening the influence of scale effect and improving the accuracy of NDVI estimation.

[0023] 3. The present invention uses the EEMD method to perform multi-scale decomposition on the NDVI time series, and combines the Bi TCN deep learning model to model the change characteristics of different time scales respectively, which can better capture the time dynamic change law of NDVI.

[0024] 4. The present invention can efficiently generate NDVI products with a 30-meter resolution for a long time series, providing important data support for applications such as vegetation growth monitoring, crop yield estimation, and ecosystem change assessment. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 It is a flowchart of an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0026] The embodiments of the present invention will be described in detail below with reference to the drawings.

[0027] Embodiment 1

[0028] As Figure 1 shown, a method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning includes:

[0029] S1: Establish an XGBoost regression model with Landsat surface reflectance and solar zenith angle as inputs and GLASS 250m NDVI as the output to estimate the NDVI of Landsat 30m clear-sky pixels.

[0030] S2: For Landsat pixels contaminated by clouds, find similar pixels within the spatio-temporal window by constructing seasonal curves, and fill them with the weighted similar pixels and GLASS NDVI values to generate a temporally continuous Landsat NDVI curve.

[0031] S3: Use the EEMD method to decompose the Landsat NDVI curve in time series to obtain the intrinsic mode functions and the residual curve.

[0032] S4: Construct a Bi TCN model to perform time series prediction on the intrinsic mode functions and the residual curve, and combine the prediction results to obtain the NDVI inversion result with high spatio-temporal resolution.

[0033] In one implementation, step S1 further includes:

[0034] Taking the southwest region as an example, first use the Kmeans clustering method to classify the GLASS NDVI time series. By analyzing the temporal variation characteristics of GLASS NDVI from 2000 to 2023, determine the optimal number of clusters. On this basis, select training samples globally according to the clustering results. To reduce the uncertainty brought by scale conversion, when selecting 250m GLASS pixels, use Global Land 30 data to ensure that all 30m Landsat pixels under this pixel have the same land cover type. Finally, about 1 million representative training sample points are selected globally.

[0035] For each training sample point, extract the corresponding Landsat red-band and near-infrared-band surface reflectance data, and at the same time obtain the corresponding solar zenith angle information. Use these data as input features and GLASS NDVI as the output label to train the XGBoost regression model. The model parameters are set as follows: the learning rate is 0.1, the maximum tree depth is 6, the minimum number of samples in leaf nodes is 3, and the regularization parameter is 0.3. Determine the optimal parameter combination through cross-validation, and finally construct a regression model that can accurately estimate the NDVI at 30m resolution.

[0036] In one implementation, step S2 includes:

[0037] For Landsat pixels contaminated by clouds, first, a 3 km × 3 km spatial window is established with the target pixel as the center. Within this window, the observation data for the five years before and after all pixels are extracted, and the seasonal variation characteristics of their NDVI are calculated. Specifically, the multi-year average of the NDVI time series for each pixel is calculated to obtain a characteristic curve representing the seasonal variation of vegetation growth.

[0038] By calculating the correlation coefficients between the seasonal curves of the target pixel and other pixels within the window, pixels with similar vegetation phenological characteristics are identified. The top 50 pixels with the highest correlation coefficients are selected as similar pixels. Then, according to the spatial distances between these similar pixels and the target pixel, the corresponding weight coefficients are calculated. These weight coefficients and the GLASS NDVI product are jointly used to determine the NDVI value of the target pixel.

[0039] In one implementation, step S3 includes:

[0040] The obtained NDVI time series is decomposed by the EEMD method. The noise amplitude is set to 0.2 times the standard deviation of the time series, and each time series is subjected to 100 times of ensemble empirical mode decomposition. EEMD decomposition can decompose the NDVI time series into several IMF components and a residual term.

[0041] Frequency analysis is performed on the decomposed IMF components and sorted in descending order of frequency. Generally, high-frequency IMF components (with a period less than 3 months) reflect short-term fluctuations, medium-frequency IMF components (with a period of 3 - 12 months) reflect seasonal changes, low-frequency IMF components (with a period greater than 12 months) reflect interannual changes, and the final residual term reflects the long-term change trend. The IMF components of different frequencies reveal the variation characteristics of NDVI at different time scales.

[0042] In one implementation, step S4 includes:

[0043] Specific Bi TCN models are constructed respectively for the different characteristics of the intrinsic mode functions and residual curves. The input of each model is the time series of the past 5 years, and the output is the prediction series for the next year. Specifically, the time convolutional network adopts a setting with a convolutional kernel size of 3 and dilation rates of [1, 2, 4, 8], and gradually expands the receptive field through multi-layer convolutional operations to capture time characteristics at different scales. The dimension of the hidden layer of the network is set to 64 to balance the expressive power and computational efficiency of the model.

[0044] In the bidirectional gated recurrent unit part, the hidden layer dimension is set to 128. This bidirectional structure can simultaneously consider the forward and backward dependencies of the time series, which helps to improve the prediction accuracy. During the model training process, the root mean square error is used as the loss function, and the Adam optimizer is used for parameter optimization with the initial learning rate set to 0.001.

[0045] To prevent the model from overfitting, a dynamic learning rate adjustment strategy is adopted: when the validation set loss has not improved for 5 consecutive epochs, the learning rate is reduced to 0.1 times the original; when the validation set loss has not improved for 10 consecutive epochs, the training process is stopped. At the same time, all samples are divided into a training set, a validation set, and a test set according to the ratio of 7:2:1.

[0046] In the application stage, first, the trained BiTCN model is used to predict the future changes of each IMF component and the residual term respectively. Then, all the prediction results are reconstructed to obtain the complete NDVI time series. To verify the performance of the model, we selected a representative study area to conduct experiments. The results show that this method can not only accurately predict the temporal changes of NDVI, but also the prediction results have good temporal continuity.

[0047] Finally, in the southwestern region, the Landsat NDVI products with a 30-meter spatial resolution and an 8-day temporal resolution from 2000 to 2022 are obtained using steps 1 to 4. Through the comparative analysis with ground observation data and existing NDVI products, the effectiveness and reliability of this method are verified.

[0048] The present invention can integrate the advantages of time series modal decomposition and machine learning techniques, solve the problems of low accuracy and poor temporal continuity in the traditional method for generating high spatio-temporal resolution NDVI products, generate high-precision, long-time series 30-meter resolution NDVI products, and provide important data support for vegetation dynamic monitoring, crop growth assessment, and ecosystem change research; integrate time series modal decomposition and machine learning techniques, make full use of the advantages of multi-source remote sensing data and advanced algorithms, consider the multi-scale change characteristics and spatial heterogeneity of the NDVI time series, and generate a normalized difference vegetation index product with high spatio-temporal resolution and good temporal continuity.

[0049] Example 2

[0050] A method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning, specifically including the following steps:

[0051] Step 1: Using the surface reflectance data obtained from Landsat satellites and the solar zenith angle as input variables, selecting the GLASS 250-meter resolution NDVI product as the training target, constructing a regression model based on the XGBoost machine learning algorithm, and achieving accurate estimation of the 30-meter resolution NDVI of Landsat clear-sky pixels.

[0052] Taking the southwest region of China as an example, this area has complex terrain and rich and diverse vegetation types, including various surface cover types such as subtropical evergreen broad-leaved forests, temperate deciduous broad-leaved forests, alpine meadows, and farmlands. First, perform Kmeans clustering analysis on the GLASS NDVI time series from 2000 to 2023. By calculating the silhouette coefficient and the clustering distance, the optimal number of clusters is finally determined to be 12, which better reflects the phenological characteristics of different vegetation types.

[0053] In each clustering category, select representative 250-meter GLASS pixels in combination with the Global Land 30 land cover data to ensure that all 30-meter Landsat pixels under this pixel have the same surface cover type, thereby reducing the uncertainty caused by scale conversion. Finally, about 1 million training sample points are screened in the southwest region.

[0054] For each training sample point, extract the surface reflectance data of the red band (0.63 - 0.69 μm) and the near-infrared band (0.76 - 0.90 μm) observed by Landsat satellites and the solar zenith angle information. The XGBoost model parameters are set as follows: the learning rate is 0.1, the maximum tree depth is 6, the minimum number of samples in leaf nodes is 3, and the regularization parameter λ is set to 0.3. Determine the optimal parameter combination through 5-fold cross-validation.

[0055] Step 2: For Landsat pixels where the NDVI value cannot be directly obtained due to cloud contamination, by analyzing the seasonal change characteristics of the NDVI of the target pixel and its surrounding pixels over the years, find reference pixels with similar vegetation growth rhythms in the time and space dimensions, and combine the observed values of these similar pixels and the GLASS NDVI product to fill the cloud-contaminated pixels in a weighted average manner to generate a continuous and complete Landsat NDVI time series.

[0056] In the study of the southwest region, first establish a 3km×3km search window centered on the target pixel. The selection of this window size is based on the fragmented terrain characteristics of this region, which should not only ensure that enough potential similar pixels are included but also consider the impact of terrain changes on vegetation characteristics. For each pixel in the window, extract its observation records for the five years before and after the period from 2000 to 2023.

[0057] To obtain stable seasonal characteristics, the 23-year observations of each pixel are grouped and superimposed according to the date, and the multi-year average NDVI seasonal change curve is calculated. By calculating the Pearson correlation coefficient between the target pixel and the seasonal curves of other pixels within the window, the top 50 pixels with the highest correlation coefficient are selected as similar pixels. Finally, considering the impact of the terrain undulation in the southwestern region on vegetation growth, both spatial distance and elevation difference are taken into account when calculating the weights.

[0058] Step 3: Use the ensemble empirical mode decomposition (EEMD) technique to decompose the obtained Landsat NDVI time series, adaptively separating the time series into intrinsic mode functions (IMFs) that reflect the change characteristics at different time scales and a residual curve that represents the long-term change trend, providing a more detailed modeling basis for subsequent time series prediction.

[0059] In the practice of the southwestern region, the noise amplitude of the EEMD method is set to 0.2 times the standard deviation of the time series, and each time series is decomposed by the ensemble empirical mode decomposition 100 times. The decomposed IMF components can be mainly divided into three categories: (1) high-frequency components (period < 3 months), reflecting short-term fluctuations; (2) medium-frequency components (period 3 - 12 months), corresponding to seasonal changes; (3) low-frequency components (period > 12 months), characterizing interannual changes. The final residual term reflects the overall change trend of NDVI.

[0060] Through the analysis of a large number of samples, it is found that the NDVI time series of different vegetation types in the southwestern region can usually be decomposed into 4 - 6 meaningful IMF components. Among them, the amplitude of the seasonal component in the farmland area is relatively large, while the seasonal fluctuations in the evergreen forest area are relatively small. This decomposition method can effectively separate the change characteristics at different time scales.

[0061] Step 4: Based on the obtained intrinsic mode functions and residual curve, construct a deep learning model (BiTCN) that integrates a temporal convolutional network and a bidirectional gated recurrent unit, simultaneously using convolutional operations to capture local temporal features and a recurrent structure to model long-term dependencies. Through the independent prediction and reconstruction integration of features at different time scales, a high-precision NDVI time series estimation result is finally obtained.

[0062] According to the characteristics of different frequency components in the southwestern region, specialized BiTCN models are constructed respectively. The input of the model is the sequence of the past 5 years, and the output is the prediction result for the next year. The temporal convolutional network adopts a multi-layer causal convolutional structure, with a convolutional kernel size of 3 and dilation rates of [1, 2, 4, 8]. The dimension of the network hidden layer is set to 64, and the dimension of the hidden layer of the bidirectional gated recurrent unit is 128. The root mean square error is used as the loss function, and the Adam optimizer is used for parameter optimization with an initial learning rate of 0.001.

[0063] To improve the generalization ability of the model, the samples are divided into a training set, a validation set, and a test set at a ratio of 7:2:1. At the same time, considering the differences in the phenological characteristics of vegetation at different altitudes in the southwestern region, a sampling strategy based on elevation stratification is added during data augmentation. When the validation set loss has not improved for 5 consecutive epochs, the learning rate is reduced to 0.1 times the original; when it has not improved for 10 consecutive epochs, the training stops.

[0064] After obtaining the prediction results of each IMF component and the residual term, the complete NDVI time series is reconstructed by summation. The 30-meter resolution NDVI products from 2000 to 2022 generated by this method are compared and verified with ground observation data. The results show that: under the complex terrain and diverse vegetation conditions in the southwestern region, this method can still maintain a high estimation accuracy (R 2 > 0.85, RMSE <0.1) and good temporal continuity.

[0065] Based on the above steps, the Landsat NDVI products with a resolution of 30 meters and an 8-day interval from 2000 to 2022 in the southwestern region are finally generated.

[0066] The above-described embodiments only represent the specific implementation manners of the present invention, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the patent of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several modifications and improvements can still be made, and these all belong to the protection scope of the present invention.

Claims

1. A method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning, characterized in that It includes the following steps: S1: Establish an XGBoost regression model with Landsat surface reflectance and solar zenith angle as inputs and GLASS 250m NDVI as the output to estimate the NDVI of Landsat 30m clear-sky pixels; S2: For Landsat pixels contaminated by clouds, find similar pixels within the spatio-temporal window by constructing seasonal curves, and fill them with the weighted similar pixels and GLASS NDVI values to generate a temporally continuous Landsat NDVI curve; S3: Use the EEMD method to perform temporal decomposition on the Landsat NDVI curve to obtain the intrinsic mode functions and the residual curve; S4: Construct a BiTCN model to perform temporal prediction on the intrinsic mode functions and the residual curve, and combine the prediction results to obtain the NDVI inversion result with high spatio-temporal resolution.

2. The method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning according to claim 1, wherein In step S1, when selecting representative samples for training the XGBoost regression model, use the Kmeans method to cluster the time-series GLASS NDVI, uniformly select sample points globally according to the clustering results, and at the same time ensure that all 30m pixels under a 250m GLASS pixel belong to the same class to weaken the scale effect.

3. The method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning according to claim 1, characterized in that In step S2, establish a 3km×3km spatial window centered on the target pixel. For all pixels within the spatial window, calculate the mean Landsat NDVI in the previous and next 5 years to obtain the seasonal NDVI curve, calculate the correlation coefficient between the seasonal NDVI curves of other pixels and the target pixel, and select the top 50 pixels with the highest correlation coefficient as similar pixels.

4. The method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning according to claim 1, characterized in that, Step S3 is specifically as follows: Set the noise amplitude during EEMD decomposition to 0.2 times the standard deviation of the time series. Sort the decomposed IMF components from high to low frequency to obtain the intrinsic mode functions, and use the remaining trend term as the residual curve.

5. The method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning according to claim 1, wherein In step S4, construct BiTCN models for predicting the intrinsic mode functions and the residual curve respectively. The input time series length is 5 years, and the output time series length is 1 year; Set the convolutional kernel size of the time convolutional network to 3, the dilation rate to [1, 2, 4, 8], the hidden layer dimension to 64, and set the hidden layer dimension of the bidirectional gated recurrent unit to 128; use the root mean square error as the loss function, and use the Adam optimizer to optimize the model. Set the initial learning rate to 0.001; when the validation set loss does not decrease for 5 consecutive epochs, reduce the learning rate to 0.1 times the original value, and when the validation set loss does not decrease for 10 consecutive epochs, stop training.

6. The method for generating a high spatio-temporal resolution normalized difference vegetation index based on modal decomposition and deep learning according to claim 5, wherein The training and application process of the BiTCN model also includes: Divide the training samples into three parts, that is, 70% of the sample set is used to train the BiTCN model, 20% of the sample set is used for model validation, and 10% of the sample set is used for independent testing of the model; Use the BiTCN model to estimate the time-series intrinsic mode functions and the residual curve, and sum all the intrinsic mode functions and the residual curve to obtain the time-series NDVI estimation result.

Images