A method for spatially mapping hay meadows in natural temperate grasslands using NDVI time series

Abstract

The present invention proposes a method for spatially mapping natural warm-season grassland mowing fields using NDVI time series, a method for spatially identifying natural warm-season grassland mowing fields based on high-spatial-resolution and dense-time-series NDVI. This method includes four steps: establishing time-series NDVI, identifying mowing fields based on machine learning, visualizing the model based on Grad-CAM, and post-processing classification results. The accuracy of the spatial range of the mowing fields produced by the present invention meets the needs of research and production, and helps China to conduct an inventory of mowing field resources in pastoral areas.

Description

Technical Field

[0001] The present invention provides a method for spatially mapping natural warm-season grassland hayfields using NDVI time series, belonging to the technical field of ecological management. Background Art

[0002] Grassland mowing is one of the main ways of grassland utilization and management. Grassland mowing refers to harvesting grassland when it grows well, drying it and using it as livestock feed. Mowing is very common in grasslands in the temperate maritime climate region of Western Europe. In addition, it also occurs in grasslands in the temperate continental climate region on the eastern edge of the Eurasian continent, the temperate continental climate region of North America, and the tropical savannah climate region of Australia. At present, some researchers have carried out monitoring of grassland mowing activities based on remote sensing technology, such as monitoring the start time and frequency of grassland mowing. These studies have experienced a transformation from single-sensor, low-resolution (>=250m), time-section (<=1 year) monitoring to multi-source sensors (such as optical and SAR data fusion), medium-high resolution (10~30m), multi-year (>=3 years) monitoring, providing data and method references for grassland management and monitoring research. However, these studies generally focus on permanent grasslands in Europe, especially the research results in Germany as the research area are the most abundant. There are very few monitoring methods for the utilization methods of natural warm-season grasslands in northern China. The main reason may be that there are significant differences in grassland management methods, policy systems, and the degree of attention paid by managers and farmers and herdsmen to the accumulation of historical data between different regions. However, with the increasing attention of China to the sustainable development of grasslands, especially the emergence of the benefits of ecological measures and ecological projects such as "three types of livestock" in grasslands, it is particularly important and urgent to establish a normalized and high-precision grassland utilization database.

[0003] Due to the huge differences in agricultural policies and climatic conditions between regions, the previous grassland mowing monitoring methods in the temperate maritime climate region are difficult to apply to the natural warm-season grassland areas in northern China. The natural warm-season grasslands in northern China are mainly distributed in arid and semi-arid areas, belonging to the temperate continental climate. Grassland mowing mainly occurs in grasslands with relatively good hydrothermal conditions, and usually only one harvest is made in a year, which is different from the grasslands in the temperate maritime climate region where multiple harvests can be made in a year. Therefore, the research on grassland management in the temperate maritime climate region not only focuses on the spatial distribution of grassland mowing, but also focuses on the earliest start time and frequency of grassland mowing - which is part of the grassland utilization intensity. And the previous grassland mowing research in China mainly focused on the spatial distribution. A fast, intelligent, and accurate monitoring method is more urgent and practical than the visual interpretation method in the routine grassland utilization monitoring. Summary of the Invention

[0004] The present invention is committed to developing a method for monitoring and spatially mapping natural hayfields in typical warm-season grasslands in northern China, and hopes that this method can be applied to the inventory of natural hayfield resources in northern China. The following problems have been solved:

[0005] (1)Realize the spatial mapping of hay meadows in temperate natural grasslands in northern China;

[0006] (2)During the spatial mapping of hay meadows, there is no need to artificially define any classification thresholds. By integrating deep learning methods, accurate, intelligent, end-to-end, and fully automated hay meadow detection can be achieved.

[0007] The present invention proposes a method for spatially identifying hay meadows in natural temperate grasslands based on high-spatial-resolution and dense-temporal NDVI. This method includes four steps: establishing temporal NDVI, hay meadow identification based on machine learning, model visualization based on Grad-CAM, and post-processing of classification results.

[0008] A method for spatially mapping hay meadows in natural temperate grasslands using NDVI time series mainly includes four steps: time series generation, classification of hay meadows and non-hay meadows, model interpretation, and model inference.

[0009] (1) Time series generation;

[0010] Collect Harmonized Landsat Sentinel-2 (HLS) data for natural hay meadow detection.

[0011] HLS data is generated by fusing and coordinating seamless surface reflectance databases from the Operational Land Imager (OLI) on Landsat-8 / 9 and the Multispectral Instrument (MSI) on Sentinel-2A / B remote sensing satellites.

[0012] For the obtained remote sensing data, first use its FMask band to remove cloud pixels, cloud shadow pixels, and cloud adjacent pixels. At the same time, download the ESA WorldCover global land use dataset, which aims to provide high-resolution global surface cover data. These data are based on Sentinel-1 and Sentinel-2 satellite data and provide surface cover data information at a resolution of 10 meters, including 11 different surface cover classes: forest, grassland, shrub, cropland, built-up, bare soil / sparse vegetation, ice / snow, water, wetland, mangrove, and moss. Extract the grassland area;

[0013] After obtaining cloud-free HLS data within the grassland area, calculate NDVI according to formula (1):

[0014] (1)

[0015] Where, and respectively represent the near-infrared band and the red band. Subsequently, linear interpolation and Savitzky-Golay filtering are used to obtain temporally continuous and spatially seamless NDVI data. First, the linear interpolation method is applied to obtain the initial interpolation of the missing NDVI values. Define the NDVI at the moment to be interpolated x m as y m , define the two observation points closest to the data to be interpolated as ( x a , y a ) and ( x b , y b ), where x a , x b represent the time points, y a , y b represent the observed NDVI values at the corresponding times. Then the initial interpolation at the moment to be interpolated x m is:

[0016] (2)

[0017] After initially filling the missing values in the NDVI sequence through the linear interpolation method, the Savitzky-Golay filtering is used to optimize the initial temporal NDVI. The Savitzky-Golay filtering retains the high-frequency features of the data by fitting polynomials of the data points, thereby reducing the NDVI data noise caused by cloud occlusion, shadows, or high atmospheric aerosol concentrations. The final NDVI sequence is generated using the savgol_filter method in SciPy;

[0018] (2) Classification of mowing pastures and non-mowing pastures;

[0019] The training set sample blocks and validation set sample blocks are obtained through visual interpretation, and the test set sample blocks are obtained through on-site inspections. Subsequently, the above vector samples are converted into raster sample points using the "vector to raster" tool. During the subsequent training, validation, and prediction processes of the model, the labels of the mowing pasture pixels are set to "1", and the labels of the non-mowing pasture pixels are set to "0".

[0020] Based on the above dataset, compare the performance of four classification models (random forest, multi-layer perceptron, 1D-CNN, and LSTM-FCN) in the classification of mowing pastures and non-mowing pastures.

[0021] Use cross-entropy loss to guide the model training. Define the loss function is:

[0022] (3)

[0023] wherein, represents the number of training samples; represents the number of classes; represents the true label of the n th sample, represents the predicted label of the n th sample. During the model training process, the Adam optimizer is used, and the initial learning rate is set to 0.001. A total of 100 rounds of model training are set. After each round of training, the validation set is used to verify the model, and the classification accuracy on the validation set is calculated. Finally, the model parameters with the highest accuracy on the validation set are used to test the classification accuracy on the test set and compare different models.

[0024] The overall accuracy OA, F1 accuracy, recall, and precision are used to evaluate the classification accuracy, and the calculation method is:

[0025] (4)

[0026] (5)

[0027] (6)

[0028] (7)

[0029] wherein, TP represents the true positive example, that is, the number of grassland pixels correctly predicted by the model; TN represents the true negative example, that is, the number of non-grassland pixels correctly predicted by the model; FP represents the false positive example, that is, the number of non-grassland pixels misclassified as grassland pixels by the model; FN represents the false negative example, that is, the number of grassland pixels misclassified as non-grassland pixels by the model.

[0030] (3) Model interpretability based on Grad-CAM of gradient-weighted mapping of class activation;

[0031] Grad-CAM calculates the gradient of the target class for the output of a specific convolutional layer. After average pooling these gradients, they are multiplied by the feature map of the convolutional layer to obtain the class activation map. Grad-CAM is applied to visually interpret the LSTM-FCN model. For the input time series data, the moments that contribute greatly to the prediction results are highlighted to qualitatively explain the model decision-making process.

[0032] (4) Model inference and post-processing;

[0033] Model inference refers to predicting the range of hay meadows in a region based on time-series remote sensing images using the optimal model parameters. The NDVI images within the ranges of "T50TMP" and "T50TMQ" are synthesized in time-series bands. According to the accuracy comparison results on the test set, the optimal model is selected, and the spatial distribution map of hay meadow pixels in the study area is calculated. To enhance the classification effect and reduce noise, median filtering is used to post-process the prediction result map.

[0034] The technical effects of the present invention are as follows:

[0035] (1) Based on the natural temperate grasslands in northern China, the present invention proposes a method for mapping the spatial distribution of hay meadows based on time-series remote sensing data. Combining remote sensing data with high spatio-temporal resolution, the technical solution of the present invention can be applied to the mapping of hay meadows with long time-series and large spatial extents, providing a technical reference for preparing a list of hay meadow resources in China;

[0036] (2) Compared with the previous threshold-based methods, the technical solution of the present invention has better generalization and higher accuracy; the technical solution of the present invention has obtained an overall accuracy of 87.03% and an F1 accuracy of 89.02% in the study area, meeting the requirements of research and production;

[0037] (3) The technical solution of the present invention has a high degree of automation and intelligence, and is more time-saving and labor-saving than the methods of field investigation or visual interpretation for mapping the spatial distribution of hay meadows;

[0038] (4) The present invention only requires optical remote sensing images and does not require other auxiliary data (such as radar data, etc.) to achieve high mapping accuracy, and the data processing process is more efficient;

[0039] (5) Based on the class activation gradient weighted mapping method, the present invention visualizes the maximum contribution time step of the time-series data input into the deep learning model, which helps to understand the decision-making of the time-series model for specific inputs. Description of the Drawings

[0040] Figure 1 is the flow chart of the present invention;

[0041] Figure 2 is the schematic diagram of the LSTM-FCN structure. Detailed Embodiments

[0042] As Figure 1 shown, the present method mainly includes 4 steps: time-series generation, classification of hay meadows and non-hay meadows, model interpretation, and model inference.

[0043] (1) Time-series generation;

[0044] The present invention collects Harmonized Landsat Sentinel-2 (HLS) data for the detection of natural hayfields. The HLS data is generated by fusing and coordinating the seamless surface reflectance databases obtained from the Landsat Imager and Multispectral Instrument carried on the Landsat-8 / 9 and Sentinel-2A / B remote sensing satellites respectively. The production process of the HLS product includes: atmospheric correction, cloud and cloud shadow masking, geospatial registration, bidirectional reflectance distribution function normalization, and spectral channel adjustment. The HLS includes two groups of products, S30 and L30, both with a spatial resolution of 30m and a temporal resolution of 2 - 3 days.

[0045] The present invention obtains all remote sensing data with image tile numbers 'T50TMP' and 'T50TMQ' from March 1, 2023, to October 31, 2023, through the EarthData API. For the obtained remote sensing data, the present invention first uses its FMask band to remove cloud pixels, cloud shadow pixels, and cloud adjacent pixels. At the same time, the present invention downloads the ESA WorldCover global land use dataset, which aims to provide high-resolution global surface cover data. These data are based on Sentinel-1 and Sentinel-2 satellite data and provide surface cover data information for 2020 and 2021 at a resolution of 10 meters, including 11 different surface cover categories: forest, grassland, shrub, cultivated land, building, bare soil / sparse vegetation area, ice and snow, water area, wetland, mangrove, and moss. The present invention extracts the grassland area, and all subsequent technical processes are based on the grassland area.

[0046] After obtaining cloud-free HLS data within the grassland area, calculate NDVI according to formula (1):

[0047] (1)

[0048] where and represent the near-infrared band and the red band respectively. Subsequently, the present invention uses linear interpolation and Savitzky-Golay filtering to obtain temporally continuous and spatially seamless NDVI data. First, apply the linear interpolation method to obtain the initial interpolation of the missing NDVI values. Define the NDVI at the interpolation time x m as y m , and define the two observation points closest to the data to be interpolated as ( x a , y a ) and ( x b ,y b ), where x a , x b represents a time point, y a , y b represents the observed NDVI value at the corresponding time. Then the initial interpolation at the time to be interpolated x m is:

[0049] (2)

[0050] After initially filling the missing values in the NDVI sequence by linear interpolation, the initial temporal NDVI is optimized using Savitzky-Golay filtering. Savitzky-Golay filtering preserves the high-frequency features of the data by fitting polynomials to the data points, thereby reducing the NDVI data noise caused by cloud occlusion, shadow, or high atmospheric aerosol concentration. The present invention uses the savgol_filter method in SciPy to generate the final NDVI sequence for subsequent classification algorithms.

[0051] (2) Classification of mowing grasslands and non-mowing grasslands;

[0052] The present invention obtained 569 training set sample blocks (including 362 mowing grassland sample blocks and 207 non-mowing grassland sample blocks) and 183 validation set sample blocks (including 122 mowing grassland sample blocks and 61 non-mowing grassland sample blocks) through visual interpretation, and obtained 156 test set sample blocks (including 123 mowing grassland sample blocks and 33 non-mowing grassland sample blocks) through on-site investigation. Subsequently, the present invention converted the above vector samples into raster sample points through the "vector to raster" tool. Finally, there are 66,839 mowing grassland pixels and 104,428 non-mowing grassland pixels in the training set; 23,954 mowing grassland pixels and 29,754 non-mowing grassland pixels in the validation set; 15,096 mowing grassland pixels and 10,523 non-mowing grassland pixels in the test set. During the subsequent training, validation, and prediction processes of the model, the label of the mowing grassland pixels is set to "1", and the label of the non-mowing grassland pixels is set to "0".

[0053] Based on the above dataset, the present invention compares the effects of four classification models (random forest, multi-layer perceptron, 1D-CNN, and LSTM-FCN) in the classification of mowing grasslands and non-mowing grasslands.

[0054] During the process of generating numerous decision trees in random forest classification, random sampling is performed on the sample observations and feature variables of the modeling dataset respectively. Each sampling result is a tree, and each tree generates rules and classification results that conform to its own attributes. Finally, the forest integrates the rules and classification results of all decision trees to achieve the classification of the random forest algorithm. In the present invention, NDVI at different times is regarded as the input feature of the random forest model. The number of trees in the random forest is set to 500 in the present invention, and the maximum feature utilization rate is set to 0.8.

[0055] A multi-layer perceptron is a feedforward artificial neural network with an input layer, a hidden layer, and an output layer, and each layer is connected to the layer above it. The multi-layer perceptron is the simplest deep network and can be used for tasks such as image recognition and text recognition. In the present invention, a multi-layer perceptron containing three linear layers is used for model classification. The output vector lengths of the first two linear layers are 500, and the last linear layer is the output layer with an output dimension of 2.

[0056] 1D-CNN is a model that applies a convolutional neural network to one-dimensional vectors and can be used for time series classification tasks. In the present invention, a 1D-CNN model containing two hidden layers and a classification layer is designed. Each hidden layer contains a 1D convolutional layer, a ReLU activation function, and a batch normalization; the classification layer contains two linear layers, and the two linear layers are connected by a ReLU activation function in the middle. The kernel size of the convolutional layer in the first hidden layer is 3, the padding length is 1, the number of input features is 1, and the number of output features is 128; the kernel size of the convolutional layer in the second hidden layer is set to 5, the padding length is 2, the number of input features is 128, and the number of output features is 256. The input dimensions of the two linear layers in the classification layer are 256 and 128 respectively, and the output dimensions are 128 and 2 respectively.

[0057] LSTM-FCN combines the fully convolutional network (FCN) with the long short term memory network (LSTM) and applies it to the deep feature learning of time series data. As Figure 2, the present invention inputs time series data into an LSTM-FCN model. This model consists of two branches. One branch is an LSTM, which is used to extract time information, and its hidden layer scale is 128. The other branch is a fully convolutional block as a feature extractor. The fully convolutional block contains three temporal convolutional blocks. Each convolutional block contains a one-dimensional convolutional layer, a batch normalization, and an activation function. Between every two convolutional blocks, there is an SE (Squeeze and Excite) block connected to adaptively calibrate the weight values of the depth feature map. The output dimension of the first convolutional block is 128, the output dimension of the second convolutional block is 256, and the output dimension of the third convolutional block is 128. A global pooling is set after the fully convolutional block to adjust the feature dimension. The features extracted by the LSTM and the features extracted by the fully convolutional block are concatenated. A SoftMax classifier is used for classification, and finally the classification result is output.

[0058] The present invention uses cross-entropy loss to guide the model training. Define the loss function as:

[0059] (3)

[0060] where, represents the number of training samples; represents the number of classes; represents the true label of the n th sample, represents the predicted label of the n th sample. In the process of model training, the present invention uses the Adam optimizer, and the initial learning rate is set to 0.001. The present invention sets a total of 100 rounds of model training. After each round of training, the validation set is used to verify the model, and the classification accuracy on the validation set is calculated. Finally, the model parameters with the highest accuracy on the validation set are taken to conduct classification accuracy tests and comparisons of different models on the test set.

[0061] The present invention uses the overall accuracy (OA), F1 accuracy, recall, and precision to evaluate the classification accuracy. The calculation methods are as follows:

[0062] (4)

[0063] (5)

[0064] (6)

[0065] (7)

[0066] Among them, TP represents the true positive, that is, the number of grassland pixels correctly predicted by the model; TN represents the true negative, that is, the number of non-grassland pixels correctly predicted by the model; FP represents the false positive, that is, the number of non-grassland pixels misclassified as grassland pixels by the model; FN represents the false negative, that is, the number of grassland pixels misclassified as non-grassland pixels by the model.

[0067] (3) Model interpretability based on Gradient-weighted Class Activation Mapping (Grad-CAM);

[0068] Grad-CAM calculates the gradient of the target class with respect to the output of a specific convolutional layer. After average pooling these gradients, they are multiplied by the feature map of the convolutional layer to obtain the class activation mapping. This process can be understood as weighting the feature map, so that the feature map that contributes more to the prediction result gets a higher weight in the mapping. Grad-CAM can generate visual explanations for any network based on the CNN model. In the present invention, Grad-CAM is applied to the visual interpretation of the LSTM-FCN model. For the input time-series data, the moments that contribute greatly to the prediction result are highlighted, and the model decision-making process is explained in a qualitative way.

[0069] (4) Model inference and post-processing;

[0070] Model inference refers to predicting the grassland range of the region based on the time-series remote sensing images using the optimal model parameters. In the present invention, the NDVI images within the "T50TMP" and "T50TMQ" ranges are synthesized into time-series bands. For example, the original NDVI image set at 81 moments is synthesized into a single remote sensing image containing 81 bands according to the chronological order. According to the accuracy comparison results on the test set mentioned above, the optimal model is selected, and the spatial distribution map of the grassland pixels in the study area is calculated. To enhance the classification effect and reduce noise, the present invention uses median filtering for post-processing of the prediction result map.

Claims

1. A method for spatially mapping hayfields in natural temperate grasslands using NDVI time series, characterized in that, It includes four steps: time series generation, classification of mowing fields and non-mowing fields, model interpretation, and model inference; (1) Time series generation; Collect Harmonized Landsat Sentinel-2 (HLS) data for natural mowing field detection; HLS data is generated by fusing and coordinating the seamless surface reflectance databases generated by the Landsat Imager and Multispectral Imager carried on Landsat-8 / 9 and Sentinel-2A / B remote sensing satellites respectively; For the obtained remote sensing data, first use its FMask band to remove cloud pixels, cloud shadow pixels, and cloud adjacent pixels; at the same time, download the ESA WorldCover global land use dataset, which aims to provide high-resolution global surface cover data; These data are based on Sentinel-1 and Sentinel-2 satellite data and provide surface cover data information at a resolution of 10 meters, including 11 different surface cover classes: forest, grassland, shrub, cultivated land, building, bare soil / sparse vegetation area, ice and snow, water area, wetland, mangrove, and moss; extract the grassland range; After obtaining cloud-free HLS data within the grassland range, calculate NDVI; Subsequently, use linear interpolation and Savitzky-Golay filtering to obtain temporally continuous and spatially seamless NDVI data; (2) Classification of mowing fields and non-mowing fields; Obtain training set sample blocks and validation set sample blocks through visual interpretation, and obtain test set sample blocks through on-site investigation; then convert the sample blocks into raster sample points through the "vector to raster" tool; during the subsequent model training, validation, and prediction processes, the label of mowing field pixels is set to "1", and the label of non-mowing field pixels is set to "0"; Based on the above dataset, compare the performance of four classification models in the classification of mowing fields and non-mowing fields; Use cross-entropy loss to guide model training; define the loss function as follows: Among them, represents the number of training samples; represents the number of categories; represents the n true label of the th sample, n represents the predicted label of the n th sample; During the model training process, the Adam optimizer is used, and the initial learning rate is set to 0.001; A total of 100 rounds of model training are set; After each round of training, the validation set is used to verify the model, and the classification accuracy on the validation set is calculated; Finally, the model parameters with the highest validation set accuracy are used to test the classification accuracy on the test set and compare different models; Use overall accuracy (OA), F1 accuracy, recall, and precision for classification accuracy evaluation; (3) Model interpretability based on Gradient-weighted Class Activation Mapping (Grad-CAM); Grad-CAM calculates the gradient of the target class for the output of a specific convolutional layer, averages these gradients after pooling, and then multiplies them by the feature map of the convolutional layer to obtain the class activation mapping; apply Grad-CAM to visually interpret the LSTM-FCN model. For the input time series data, highlight the moments that contribute greatly to the prediction results, and qualitatively explain the model decision-making process; (4) Model inference and post-processing; Model inference refers to using the optimal model parameters to predict the mowing field range in the area based on time series remote sensing images; perform temporal band synthesis on the NDVI images within the "T50TMP" and "T50TMQ" ranges, select the best model according to the accuracy comparison results on the test set, and calculate the spatial distribution map of mowing field pixels in the study area; to enhance the classification effect and reduce noise, use median filtering for post-processing of the prediction result map.

2. A method for spatially mapping mowing fields of natural temperate grasslands using NDVI time series according to claim 1, characterized in that, In step (1), calculate NDVI according to the formula: Among them, and respectively represent the near-infrared band and the red light band.

3. A method for spatially mapping the mowing fields of natural warm steppes using NDVI time series according to claim 1, characterized in that, In step (1), the specific method for obtaining temporally continuous and spatially seamless NDVI data by using linear interpolation and Savitzky-Golay filtering is as follows: The initial interpolation of the missing NDVI value is obtained by using the linear interpolation method; define the moment to be interpolated x m of NDVI as y m , define the two observation points closest to the data to be interpolated as ([[]] x a , y a ) and ([[]] x b , y b ), where x a , x b represents the time point, y a , y b represents the observed NDVI value at the corresponding time; then the initial interpolation of the moment to be interpolated x m is: After the missing values in the NDVI sequence are initially filled by the linear interpolation method, the initial temporal NDVI is optimized using the Savitzky-Golay filter; the Savitzky-Golay filter retains the high-frequency features of the data by fitting polynomials to the data points, thereby reducing the NDVI data noise caused by cloud occlusion, shadows, or high atmospheric aerosol concentration; the final NDVI sequence is generated using the savgol_filter method in SciPy.

4. A method for spatially mapping the mowing fields of natural warm steppes using NDVI time series according to claim 1, characterized in that, In step (2), the four classification models are: Random Forest, Multilayer Perceptron, 1D-CNN, and LSTM-FCN.

5. A method for spatially mapping hayfields in natural warm steppes using NDVI time series according to claim 1, characterized in that, In step (2), the overall accuracy OA, F1 accuracy, recall, and precision are used to evaluate the classification accuracy, and the calculation method is: Among them, TP represents the true positive, that is, the number of grassland pixels correctly predicted by the model; TN represents the true negative, that is, the number of non-grassland pixels correctly predicted by the model; FP represents the false positive, that is, the number of non-grassland pixels misclassified as grassland pixels by the model; FN represents the false negative, that is, the number of grassland pixels misclassified as non-grassland pixels by the model.

Images