High-temporal-spatial-resolution normalized vegetation difference index generation method based on modal decomposition and deep learning

Abstract

The invention discloses a high-temporal-spatial-resolution normalized vegetation difference index generation method based on modal decomposition and deep learning, and the method comprises the steps: building an XGBoost regression model with Landsat surface reflectance and solar zenith angle as input and GLASS 250-meter NDVI as output, and estimating the NDVI under Landsat 30-meter clear sky pixels; for the Landsat pixels polluted by the cloud, similar pixels are searched in a space-time window by constructing a seasonal curve, the similar pixels are weighted and filled with a GLASS NDVI value, and a Landsat NDVI curve with continuous time is generated; performing time sequence decomposition on the Landsat NDVI curve by adopting an EEMD (Ensemble Empirical Mode Decomposition) method to obtain an intrinsic mode function and a residual curve; constructing a Bi TCN model to perform time sequence prediction on the intrinsic mode function and the residual curve, and combining prediction results to obtain a high-temporal-spatial-resolution NDVI inversion result; the NDVI time sequence with the 30-meter resolution ratio can be accurately estimated, and important data support is provided for applications such as vegetation growth monitoring, crop yield estimation and ecological system change estimation.

Description

Technical Field

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

[0002] Normalized Difference Vegetation Index (NDVI) is an important parameter to characterize the coverage and growth status of surface vegetation, and is widely used in crop monitoring, ecological environment assessment, climate change research and other fields. At present, obtaining NDVI based on satellite remote sensing data has become the main means of obtaining NDVI at regional and even global scales. However, the existing NDVI products still have many limitations in terms of temporal and spatial resolution and accuracy. The existing NDVI products mainly include two categories: one is products generated based on medium and low spatial resolution satellite data such as MODI S, such as GLASS NDVI products with a resolution of 250 meters and MOD13 products with a resolution of 500 meters. Although this type of product has a high temporal resolution and good temporal continuity, its coarse spatial resolution makes it difficult to depict the fine spatial heterogeneity of surface vegetation, which limits its application in fine-scale vegetation monitoring. The other type is 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, which makes it difficult to meet the needs of dynamic monitoring.

[0003] In order to improve the spatiotemporal resolution of NDVI products, researchers have proposed a variety of methods: (1) Data fusion-based methods, such as STARFM and ESTARFM algorithms, generate NDVI products with high spatiotemporal resolution by fusing high temporal resolution and high spatial resolution remote sensing data. However, these methods require high-quality reference images and have low prediction accuracy during periods of drastic surface changes; (2) Time interpolation-based methods, such as Whittaker smoother and Savitzky-Golay filter, reconstruct continuous NDVI time series by interpolating discrete NDVI observations. However, these methods often ignore the spatial correlation of NDVI and have poor reconstruction effects during periods of sparse observations; (3) Machine learning-based methods, such as random forests and deep learning, predict by establishing the relationship between NDVI and other remote sensing features. However, these methods often require a large number of high-quality training samples, and the temporal continuity of the prediction results is difficult to guarantee. In addition, NDVI time series has obvious multi-scale characteristics, containing 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, which makes it difficult to fully utilize the variation characteristics of different time scales, affecting the accuracy of the prediction results.

[0004] Therefore, it is urgent to develop new methods that can fully utilize the advantages of multi-source remote sensing data and effectively extract and utilize the multi-scale characteristics of NDVI time series to generate NDVI products with high spatiotemporal resolution and good temporal continuity, thus providing data support for refined vegetation monitoring. Summary of the invention

[0005] In order to solve the problems existing in the prior art, the purpose of the present invention is to provide a method for generating a high spatiotemporal resolution normalized vegetation difference index based on modal decomposition and deep learning. The present invention can accurately estimate the NDVI time series with a resolution of 30 meters, and provide 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 vegetation difference index with high temporal and spatial resolution, comprising the following steps:

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

[0008] S2: For Landsat pixels polluted by clouds, similar pixels are found in the spatiotemporal window by constructing seasonal curves, and weighted similar pixels are filled with GLASS NDVI values ​​to generate a time-continuous Landsat NDVI curve;

[0009] S3: The Landsat NDVI curve is decomposed into time series using the EEMD method to obtain the intrinsic mode function and residual curve;

[0010] S4: A Bi TCN model is constructed to perform time series prediction of the intrinsic mode function and residual curve, and the prediction results are combined to obtain high temporal and spatial resolution NDVI inversion results.

[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 sample points are uniformly selected globally based on the clustering results, while ensuring that all 30-meter pixels under a 250-meter GLASS pixel belong to the same class to reduce 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 in the spatial window, the average Landsat NDVI values ​​of the past five years are calculated to obtain a seasonal NDVI curve, and the correlation coefficients of the seasonal NDVI curves of other pixels and the target pixel are calculated, and the top 50 pixels with the highest correlation coefficients 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 IMF components obtained by decomposition are sorted from high to low according to frequency to obtain the intrinsic mode function, and the remaining trend term is used as the residual curve.

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

[0016] The convolution kernel size of the temporal convolutional network is set to 3, the expansion 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, 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] The training samples are divided into three parts, namely, 70% of the sample set is used to train the BiTCN model, 20% of the sample set is used to verify the model, and 10% of the sample set is used for independent testing of the model;

[0019] The BiTCN model is used to estimate the time series intrinsic mode functions and residual curves, and the time series NDVI estimation results are obtained by summing up all the intrinsic mode functions and residual curves.

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

[0021] 1. The present invention integrates time series modal decomposition and machine learning technology, fully utilizes the advantages of multi-source remote sensing data, and effectively solves the problems of low accuracy and poor temporal continuity of traditional methods when generating high spatiotemporal resolution NDVI products.

[0022] 2. The present invention selects representative training samples through Kmeans clustering and adopts a spatiotemporal filling strategy based on similarity, which effectively weakens the impact of scale effect and improves the accuracy of NDVI estimation.

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

[0024] 4. The present invention can efficiently generate long-term series NDVI products with a resolution of 30 meters, 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 is a flow chart of an embodiment of the present invention. DETAILED DESCRIPTION

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

[0027] Example 1

[0028] like Figure 1 As shown, a method for generating a high spatiotemporal resolution normalized vegetation difference 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 input and GLASS 250m NDVI as output to estimate the NDVI of Landsat 30m clear sky pixels;

[0030] S2: For Landsat pixels polluted by clouds, similar pixels are found in the spatiotemporal window by constructing seasonal curves, and weighted similar pixels are filled with GLASS NDVI values ​​to generate a time-continuous Landsat NDVI curve;

[0031] S3: The Landsat NDVI curve is decomposed into time series using the EEMD method to obtain the intrinsic mode function and residual curve;

[0032] S4: A Bi TCN model is constructed to perform time series prediction of the intrinsic mode function and residual curve, and the prediction results are combined to obtain high temporal and spatial resolution NDVI inversion results.

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

[0034] Taking the southwest region as an example, the Kmeans clustering method was first used to classify the GLASS NDVI time series. By analyzing the temporal variation characteristics of GLASS NDVI from 2000 to 2023, the optimal number of clusters was determined. On this basis, training samples were selected globally based on the clustering results. In order to reduce the uncertainty caused by scale conversion, when selecting 250-meter GLASS pixels, GLOBALL and 30 data was used to ensure that all 30-meter Landsat pixels under the pixel had the same surface cover type. Finally, about 1 million representative training sample points were selected worldwide.

[0035] For each training sample point, the corresponding Landsat red band and near-infrared band surface reflectance data were extracted, and the corresponding solar zenith angle information was obtained. Using these data as input features and GLASS NDVI as output labels, the XGBoost regression model was trained. The model parameters were set as follows: learning rate was 0.1, maximum tree depth was 6, minimum number of leaf node samples was 3, and regularization parameter was 0.3. The optimal parameter combination was determined by cross-validation, and finally a regression model that can accurately estimate 30-meter resolution NDVI was constructed.

[0036] In one embodiment, step S2 includes:

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

[0038] By calculating the correlation coefficient between the target pixel and the seasonal curve of other pixels in the window, pixels with similar vegetation phenology characteristics are identified. The top 50 pixels with the highest correlation coefficients are selected as similar pixels. Then, the corresponding weight coefficients are calculated based on the spatial distance between these similar pixels and the target pixel. These weight coefficients are used together with the GLASS NDVI product to determine the NDVI value of the target pixel.

[0039] In one embodiment, 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 100 ensemble empirical mode decompositions are performed on each time series. EEMD decomposition can decompose the NDVI time series into several IMF components and a residual term.

[0041] The decomposed IMF components are subjected to frequency analysis and sorted from high to low frequency. Generally, high-frequency IMF components (period less than 3 months) reflect short-term fluctuations, medium-frequency IMF components (period 3-12 months) reflect seasonal changes, low-frequency IMF components (period greater than 12 months) reflect interannual changes, and the final residual term reflects the long-term trend. IMF components of different frequencies reveal the variation characteristics of NDVI on different time scales.

[0042] In one embodiment, step S4 includes:

[0043] According to the different characteristics of the intrinsic mode function and the residual curve, a special Bi TCN model is constructed respectively. The input of each model is the time series of the past 5 years, and the output is the forecast sequence of the next year. Specifically, the temporal convolutional network adopts a convolution kernel size of 3 and an expansion rate of [1,2,4,8]. Through multi-layer convolution operations, the receptive field is gradually expanded to capture time features of different scales. The hidden layer dimension 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 accuracy of the prediction. During the model training process, the root mean square error is used as the loss function, the Adam optimizer is used for parameter optimization, and the initial learning rate is set to 0.001.

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

[0046] In the application stage, the future changes of each IMF component and residual term are first predicted by the trained BiTCN model. Then, all the prediction results are reconstructed to obtain a complete NDVI time series. In order 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, the Landsat NDVI products with a spatial resolution of 30 meters and a temporal resolution of 8 days were obtained in the southwest region using steps 1 to 4 from 2000 to 2022. The effectiveness and reliability of this method were verified by comparative analysis with ground observation data and existing NDVI products.

[0048] The present invention can integrate the advantages of time series modal decomposition and machine learning technology, solve the problems of low accuracy and poor temporal continuity in generating high spatiotemporal resolution NDVI products by traditional methods, and generate high-precision, long-time series 30-meter resolution NDVI products, providing important data support for vegetation dynamic monitoring, crop growth assessment and ecosystem change research; integrate time series modal decomposition and machine learning technology, 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 NDVI time series, and generate normalized vegetation difference index products with high spatiotemporal resolution and good temporal continuity.

[0049] Example 2

[0050] A method for generating a high spatiotemporal resolution normalized vegetation difference index based on modal decomposition and deep learning specifically includes the following steps:

[0051] Step 1: Using the surface reflectance data and solar zenith angle obtained by the Landsat satellite as input variables, the GLASS 250-meter resolution NDVI product is selected as the training target, and a regression model based on the XGBoost machine learning algorithm is constructed to accurately estimate the 30-meter resolution NDVI of Landsat clear sky pixels;

[0052] Taking the southwest region of my country as an example, the terrain in this region is complex and the vegetation types are rich and diverse, including subtropical evergreen broad-leaved forests, temperate deciduous broad-leaved forests, plateau meadows, farmland and other surface cover types. First, Kmeans cluster analysis was performed on the GLASSNDVI time series from 2000 to 2023. By calculating the silhouette coefficient and cluster spacing, the optimal number of clusters was finally determined to be 12, which better reflects the phenological characteristics of different vegetation types.

[0053] In each cluster category, a representative 250-meter GLASS pixel was selected in combination with the GLOBAL Land 30 land cover data to ensure that all 30-meter Landsat pixels under the pixel had the same surface cover type, thereby reducing the uncertainty caused by scale conversion. Finally, about 1 million training sample points were screened in the southwest region.

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

[0055] Step 2: For Landsat pixels that cannot directly obtain NDVI values ​​due to cloud pollution, the seasonal variation characteristics of NDVI of the target pixel and its surrounding pixels over the years are analyzed to find reference pixels with similar vegetation growth rhythms in time and space dimensions. The observation values ​​of these similar pixels are combined with the GLASS NDVI product, and the cloud-polluted pixels are filled in by weighted averaging to generate a Landsat NDVI time series that is continuous and complete in time.

[0056] In the study of the southwest region, a 3km×3km search window centered on the target pixel was first established. The window size was selected based on the characteristics of the fragmented terrain in the region, ensuring that enough potential similar pixels were included and considering the impact of terrain changes on vegetation characteristics. For each pixel in the window, the observation records of the five years before and after the period of 2000-2023 were extracted.

[0057] In order to obtain stable seasonal characteristics, the 23-year observation values ​​of each pixel were grouped and superimposed by date, and the multi-year average NDVI seasonal change curve was calculated. By calculating the Pearson correlation coefficient between the target pixel and the seasonal curves of other pixels in the window, the top 50 pixels with the highest correlation coefficient were selected as similar pixels. Finally, considering the impact of terrain undulations on vegetation growth in the southwest region, spatial distance and elevation differences were considered simultaneously when calculating weights.

[0058] Step 3: Use the ensemble empirical mode decomposition technique to decompose the obtained Landsat NDVI time series, adaptively separate the time series into intrinsic mode functions reflecting the change characteristics of different time scales, and residual curves representing long-term change trends, providing a more detailed modeling basis for subsequent time series prediction;

[0059] In the practice in the southwest region, the noise amplitude of the EEMD method is set to 0.2 times the standard deviation of the time series, and 100 ensemble empirical mode decompositions are performed on each time series. The IMF components obtained by decomposition 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), representing interannual changes. The final residual term reflects the overall 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 southwest region can usually be decomposed into 4-6 meaningful IMF components. Among them, the seasonal component amplitude of the farmland area is large, while the seasonal fluctuation of the evergreen forest area is relatively small. This decomposition method can effectively separate the change characteristics on different time scales.

[0061] Step 4: Based on the intrinsic mode function and residual curve obtained by decomposition, a deep learning model integrating temporal convolutional network and bidirectional gated recurrent unit (Bi TCN) is constructed. At the same time, the convolution operation is used to capture the local temporal characteristics and the cyclic structure is used to model the long-term dependencies. Through the independent prediction and reconstruction integration of the characteristics of 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 southwest region, a special Bi TCN model is constructed. The model input is the past 5-year sequence, and the output is the prediction result for the next year. The temporal convolutional network adopts a multi-layer causal convolution structure, with a convolution kernel size of 3 and an expansion rate of [1, 2, 4, 8]. The hidden layer dimension of the network is set to 64, and the hidden layer dimension 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 were divided into training set, validation set and test set in a ratio of 7:2:1. At the same time, considering the differences in vegetation phenology characteristics at different altitudes in the southwest region, a sampling strategy based on elevation stratification was added during data enhancement. When the validation set loss did not improve for 5 consecutive epochs, the learning rate was reduced to 0.1 times the original; when it did not improve for 10 consecutive epochs, the training was stopped.

[0064] After obtaining the prediction results of each IMF component and residual term, the complete NDVI time series is reconstructed by summing up. The 30-meter resolution NDVI products generated by this method from 2000 to 2022 were compared with ground observation data for verification. The results show that under the complex terrain and diverse vegetation conditions in the southwest 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 30-meter, 8-day Landsat NDVI product for the southwest region from 2000 to 2022 was finally generated.

[0066] The above-mentioned embodiments only express the specific implementation of the present invention, and the description thereof is relatively specific and detailed, but it cannot be understood as limiting the scope of the present invention. It should be pointed out that, for ordinary technicians in this field, several variations and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A method for generating a high spatiotemporal resolution normalized vegetation difference index based on modal decomposition and deep learning, characterized in that: The following steps are involved: S1: Establish an XGBoost regression model with Landsat surface reflectance and solar zenith angle as input and GLASS 250m NDVI as output to estimate the NDVI of Landsat 30m clear sky pixels; S2: For Landsat pixels polluted by clouds, similar pixels are found in the spatiotemporal window by constructing seasonal curves, and weighted similar pixels are filled with GLASS NDVI values ​​to generate a time-continuous Landsat NDVI curve; S3: The Landsat NDVI curve is decomposed into time series using the EEMD method to obtain the intrinsic mode function and residual curve; S4: A BiTCN model is constructed to perform time series prediction of the intrinsic mode function and residual curve, and the prediction results are combined to obtain high temporal and spatial resolution NDVI inversion results.

2. The method for generating a high spatiotemporal resolution normalized vegetation difference index based on modal decomposition and deep learning according to claim 1, characterized in that: 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. Sample points are selected evenly around the world based on the clustering results, and at the same time, all 30-meter pixels under a 250-meter GLASS pixel are guaranteed to belong to the same class to reduce the scale effect.

3. The method for generating a high spatiotemporal resolution normalized vegetation difference index based on modal decomposition and deep learning according to claim 1, characterized in that: In step S2, a 3 km × 3 km spatial window centered on the target pixel is established. For all pixels in the spatial window, the average Landsat NDVI values ​​of the past and present five years are calculated to obtain the seasonal NDVI curve. The correlation coefficients of the seasonal NDVI curves of other pixels and the target pixel are calculated, and the top 50 pixels with the highest correlation coefficients are selected as similar pixels.

4. The method for generating a high spatiotemporal resolution normalized vegetation difference index based on modal decomposition and deep learning according to claim 1, characterized in that: Step S3 is as follows: The noise amplitude during EEMD decomposition is set to 0.2 times the standard deviation of the time series. The IMF components obtained by decomposition are sorted from high to low according to frequency to obtain the intrinsic mode function, and the remaining trend term is used as the residual curve.

5. The method for generating a high spatiotemporal resolution normalized vegetation difference index based on modal decomposition and deep learning according to claim 1, characterized in that: In step S4, BiTCN models for the prediction of intrinsic mode functions and residual curves are constructed respectively, with the input time series length being 5 years and the output time series length being 1 year; The convolution kernel size of the temporal convolutional network is set to 3, the expansion 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, 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.

6. The method for generating a high spatiotemporal resolution normalized vegetation difference index based on modal decomposition and deep learning according to claim 5, characterized in that: The training and application process of the BiTCN model also includes: The training samples are divided into three parts, namely, 70% of the sample set is used to train the BiTCN model, 20% of the sample set is used to verify the model, and 10% of the sample set is used for independent testing of the model; The BiTCN model is used to estimate the time series intrinsic mode functions and residual curves, and the time series NDVI estimation results are obtained by summing up all the intrinsic mode functions and residual curves.

Images