Method for carrying out natural warm grassland mowing space mapping by using NDVI time sequence

Abstract

The invention provides a method for drawing a natural warm grassland mowing grassland space by using an NDVI (normalized difference vegetation index) time sequence, and relates to a natural warm grassland mowing grassland space identification method based on high spatial resolution and dense time sequence NDVI. The method comprises four steps of time sequence NDVI establishment, machine learning-based grass mowing field identification, Grad-CAM-based model visualization and classification result post-processing. The accuracy of the space range of the produced grass mowing field meets the requirements of research and production, and resource checking of the grass mowing field in the pasturing area in China is facilitated.

Description

Technical Field

[0001] The invention provides a method for spatial mapping of natural temperate grassland grazing grounds by utilizing NDVI time series, and belongs 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 and using it as livestock feed after drying. Mowing is very common in grasslands in the temperate maritime climate zone of Western Europe. It also occurs in grasslands in the temperate continental climate zone of the eastern edge of Eurasia, the temperate continental climate zone of North America, and the tropical savanna climate zone of Australia. At present, some researchers have carried out grassland mowing activity monitoring based on remote sensing technology, such as monitoring the start time and frequency of grassland mowing. These studies have experienced a transition from single sensor, low resolution (>=250m), time section (<=1 year) monitoring to multi-source sensor (such as optical, SAR data fusion), medium and high resolution (10~30m), and 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 in Germany, where the most abundant results are obtained. Monitoring methods for the use of natural temperate grasslands in northern China are very rare. The main reason may be that there are significant differences in grassland management methods, policies and systems, and the importance of managers and farmers and herdsmen to the accumulation of historical data between different regions. However, as China pays more attention to the sustainable development of grasslands, especially the benefits of ecological measures such as the "three pastoral" and ecological projects, it is particularly important and urgent to establish a normalized and high-precision grassland use database.

[0003] Due to the huge differences in agricultural policies and climatic conditions between regions, the grassland cutting monitoring methods previously used in temperate maritime climate zones are difficult to apply to natural temperate grassland areas in northern China. Natural temperate grasslands in northern China are mainly distributed in arid and semi-arid areas and belong to temperate continental climates. Grassland cutting mainly occurs in grasslands with relatively good water and heat conditions, and is usually harvested only once a year, which is different from grasslands in temperate maritime climate zones that can be harvested multiple times a year. Therefore, in addition to focusing on the spatial distribution of grassland cutting, research on grassland management in temperate maritime climate zones also focuses on the earliest start time and frequency of grassland cutting, which is part of grassland utilization intensity. Previous research on grassland cutting in China focused mainly on spatial distribution. Fast, intelligent, and accurate monitoring methods are more urgent and practical than visual interpretation methods in routine grassland utilization monitoring. Summary of the invention

[0004] The present invention is dedicated to developing a natural grassland monitoring and spatial mapping method in typical temperate grasslands in northern China, and it is expected that this method can be applied to the inventory of natural grassland resources in northern China. The following problems have been solved: (1) To achieve spatial mapping of temperate natural grasslands in northern China; (2) There is no need to define any classification thresholds during the grazing ground spatial mapping process. By integrating deep learning methods, accurate, intelligent, end-to-end, and fully automated grazing ground detection can be achieved.

[0005] The present invention proposes a natural temperate grassland grazing ground spatial identification method based on high spatial resolution and dense time series NDVI. The method includes four steps: establishing time series NDVI, grazing ground identification based on machine learning, model visualization based on Grad-CAM, and post-processing of classification results.

[0006] A method for spatial mapping of natural temperate grassland grazing areas using NDVI time series mainly includes four steps: time series generation, classification of grazing areas and non-grazing areas, model interpretation and model reasoning.

[0007] (1) Time series generation; Collect Harmonized Landsat Sentinel-2 and HLS data for natural grassland detection.

[0008] HLS data is a seamless surface reflectance database generated by fusing and coordinating the land imagers and multi-spectrometers carried by Landsat-8 / 9 and Sentinel-2A / B remote sensing satellites respectively.

[0009] The remote sensing data obtained first uses its FMask band to remove cloud pixels, cloud shadow pixels and cloud adjacent pixels. At the same time, download the ESAWorldCover global land use dataset, which aims to provide high-resolution global land cover data. These data are based on Sentinel-1 and Sentinel-2 satellite data, providing land cover data information at a resolution of 10 meters, including 11 different land cover categories: forest land, grassland, shrubs, cultivated land, buildings, bare soil / sparse vegetation areas, ice and snow, water areas, wetlands, mangroves and mosses. Extract the grassland range; After obtaining HLS data of cloud-free and grassland areas, NDVI is calculated according to formula (1): (1) in, and Represent the near infrared band and the red light band respectively. Subsequently, linear interpolation and Savitzky-Golay filtering are used to obtain time-series 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 time to be interpolated x mThe NDVI is y m , define the two observation points closest to the data to be interpolated as ( x a , y a )and( x b , y b ),in x a , x b Indicates the time point, y a , y b represents the observed NDVI value at the corresponding time. x m The initial interpolation value is: (2) After initially filling the missing values ​​in the NDVI series by linear interpolation, the Savitzky-Golay filter is used to optimize the initial time series NDVI. The Savitzky-Golay filter retains the high-frequency characteristics of the data by fitting a polynomial of the data points, thereby reducing the NDVI data noise caused by cloud cover, shadows, or high atmospheric aerosol concentrations. The savgol_filter method in SciPy is used to generate the final NDVI series; (2) classification of pastures into grazing and non-grazing areas; The training set sample blocks and validation set sample blocks were obtained through visual interpretation, and the test set sample blocks were obtained through field investigation. The vector samples were then converted into raster sample points using the "vector to raster" tool. In the subsequent model training, validation, and prediction processes, the labels of grassland pixels were set to "1" and the labels of non-grazing grassland pixels were set to "0".

[0010] Based on the above dataset, the effects of four classification models (random forest, multi-layer perceptron, 1D-CNN and LSTM-FCN) in the classification of grazing land and non-grazing land were compared.

[0011] Use cross entropy loss to guide model training. Define the loss function for: (3) in, Indicates the number of training samples; Indicates the number of categories; Indicates n The true labels of samples, Indicates nThe predicted labels of samples. The Adam optimizer is used in the model training process, 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 model is verified using the validation set, 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 perform classification accuracy tests on the test set and compare different models.

[0012] The overall accuracy OA, F1 accuracy, recall rate recall and precision rate precision are used to evaluate the classification accuracy. The calculation method is: (4) (5) (6) (7) Among them, TP represents true positive examples, that is, the number of grassland pixels correctly predicted by the model; TN represents true negative examples, that is, the number of non-grazing grassland pixels correctly predicted by the model; FP represents false negative examples, that is, the number of non-grazing grassland pixels mistakenly classified as grassland pixels by the model; FN represents false negative examples, that is, the number of grassland pixels mistakenly classified as non-grazing grassland pixels by the model.

[0013] (3) Model interpretability based on class activation gradient weighted mapping Grad-CAM; Grad-CAM calculates the gradient of the target category for the output of a specific convolutional layer, averages and pools these gradients, and then multiplies them with the feature map of the convolutional layer to obtain the class activation map. Grad-CAM is applied to the visual interpretation of the LSTM-FCN model. For the input time series data, the moments that contribute most to the prediction results are highlighted, and the model decision process is explained in a qualitative way.

[0014] (4) Model reasoning and post-processing; Model inference refers to the use of the best model parameters to predict the range of regional grasslands based on time-series remote sensing images. The NDVI images within the range of "T50TMP" and "T50TMQ" are synthesized in time series bands, and the best model is selected based on the accuracy comparison results on the test set to calculate the spatial distribution map of grassland pixels in the study area. In order to enhance the classification effect and reduce noise, the median filter is used to post-process the prediction result map.

[0015] The technical effects of the present invention are as follows: (1) Based on the natural temperate grassland in northern China, this paper proposes a method for spatial mapping of grassland based on time-series remote sensing data. Combined with remote sensing data with high temporal and spatial resolution, the technical solution of this paper can be applied to grassland mapping applications with long time series and large spatial range, providing a technical reference for the preparation of grassland resource inventories in my country. (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 achieved an overall accuracy of 87.03% and an F1 accuracy of 89.02% in the study area, which meets the needs of research and production. (3) The technical solution of the present invention has a high degree of automation and intelligence, which is more time-saving and labor-saving than the method of field investigation or visual interpretation and mapping of the spatial distribution of pastures; (4) The present invention only requires optical remote sensing images, without other auxiliary data (such as radar data, etc.) to achieve higher mapping accuracy, and the data processing process is more efficient; (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 of the time series model for specific inputs. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 is a flow chart of the present invention; Figure 2 This is a schematic diagram of the LSTM-FCN structure. DETAILED DESCRIPTION

[0017] like Figure 1 As shown in the figure, this method mainly includes four steps: time series generation, classification of pasture and non-pasture, model interpretation and model reasoning.

[0018] (1) Time series generation; The present invention collects Landsat-Sentinel-2 harmonized data (Harmonized Landsat Sentinel-2, HLS) data for natural grassland detection. HLS data is a seamless surface reflectance database generated by fusing and coordinating the land imagers and multi-spectrometers carried by Landsat-8 / 9 and Sentinel-2A / B remote sensing satellites respectively. The HLS product production process includes: atmospheric correction, cloud and cloud shadow masking, geospatial alignment, bidirectional reflectance distribution function normalization and spectral channel adjustment. HLS includes two groups of products, S30 and L30, both of which have a spatial resolution of 30m and a temporal resolution of 2 to 3 days.

[0019] The present invention obtains all remote sensing data from March 1, 2023 to October 31, 2023, with image block numbers 'T50TMP' and 'T50TMQ' through the EarthData API. For the remote sensing data obtained, 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 ESAWorldCover 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: woodland, grassland, shrubs, cultivated land, buildings, bare soil / sparse vegetation areas, ice and snow, waters, wetlands, mangroves and mosses. The present invention extracts the grassland range, and the subsequent technical processes are all based on the grassland range.

[0020] After obtaining HLS data of cloud-free and grassland areas, NDVI is calculated according to formula (1): (1) in, and Represent the near infrared band and the red light band respectively. Subsequently, the present invention uses linear interpolation and Savitzky-Golay filtering to obtain time-series continuous and spatially seamless NDVI data. First, the linear interpolation method is applied to obtain the initial interpolation of the missing NDVI value. Define the time to be interpolated x m The NDVI is y m , define the two observation points closest to the data to be interpolated as ( x a , y a )and( x b , y b ),in x a , x b Indicates the time point, y a , y b represents the observed NDVI value at the corresponding time. x m The initial interpolation value is: (2) After initially filling the missing values ​​in the NDVI sequence by linear interpolation, the initial time series NDVI is optimized by Savitzky-Golay filtering. Savitzky-Golay filtering retains the high-frequency characteristics of the data by fitting a polynomial of the data points, thereby reducing the NDVI data noise caused by cloud cover, shadows, or high atmospheric aerosol concentrations. The present invention uses the savgol_filter method in SciPy to generate the final NDVI sequence for subsequent classification algorithms.

[0021] (2) classification of pastures into grazing and non-grazing areas; The present invention obtained 569 training set sample blocks (including 362 pasture sample blocks and 207 non-pasture sample blocks) and 183 validation set sample blocks (including 122 pasture sample blocks and 61 non-pasture sample blocks) through visual interpretation, and 156 test set sample blocks (including 123 pasture sample blocks and 33 non-pasture sample blocks) through field 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 pasture pixels and 104,428 non-pasture pixels in the training set; 23,954 pasture pixels and 29,754 non-pasture pixels in the validation set; and 15,096 pasture pixels and 10,523 non-pasture pixels in the test set. In the subsequent model training, validation and prediction process, the label of the pasture pixel is set to "1" and the label of the non-pasture pixel is set to "0".

[0022] Based on the above data set, the present invention compares the effects of four classification models (random forest, multi-layer perceptron, 1D-CNN and LSTM-FCN) in the classification of grazing land and non-grazing land.

[0023] In the process of generating many decision trees, random forest classification randomly samples the sample observations and feature variables of the modeling data set. Each sampling result is a tree, and each tree will generate rules and classification results that conform to its own attributes. The forest finally integrates the rules and classification results of all decision trees to achieve the classification of the random forest algorithm. The present invention regards NDVI at different times as the input features of the random forest model. The present invention sets the number of trees in the random forest to 500 and the maximum feature utilization rate to 0.8.

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

[0025] 1D-CNN is a model that applies convolutional neural networks to one-dimensional vectors and can be used for time series classification tasks. The present invention designs a 1D-CNN model comprising two hidden layers and a classification layer, each of which comprises a 1D convolutional layer, a ReLU activation function and a batch normalization; the classification layer comprises two linear layers, and the two linear layers are connected by a ReLU activation function. 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.

[0026] LSTM-FCN combines the Fully Convolutional Network (FCN) with the Long Short Term Memory Network (LSTM) and applies it to deep feature learning of time series data. Figure 2 The present invention inputs the time series data into the LSTM-FCN model, which includes two branches, one branch is LSTM, which is used to extract time information, and its hidden layer scale is 128; the other branch is a full convolution block as a feature extractor, and the full convolution block includes three time convolution blocks, each of which includes a one-dimensional convolution layer, a batch normalization and an activation function. Every two convolution blocks are connected by a SE (Squeeze and Excite) block to adaptively calibrate the weight value of the deep feature map. The output dimension of the first convolution block is 128, the output dimension of the second convolution block is 256, and the output dimension of the third convolution block is 128. A global pooling is set after the full convolution block to adjust the feature dimension, and the features extracted by the LSTM are concatenated with the features extracted by the full convolution block, and the SoftMax classifier is used for classification, and the classification result is finally output.

[0027] The present invention uses cross entropy loss to guide model training. Define the loss function for: (3) in, Indicates the number of training samples; Indicates the number of categories; Indicates n The true labels of samples, Indicates nThe predicted labels of samples. The present invention uses the Adam optimizer in the model training process, 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 model is verified using the validation set, and the classification accuracy on the validation set is calculated. Finally, the present invention takes the model parameters with the highest accuracy on the validation set to perform classification accuracy tests on the test set and compare different models.

[0028] The present invention uses overall accuracy (OA), F1 accuracy, recall rate and precision rate to evaluate the classification accuracy, and the calculation method is: (4) (5) (6) (7) Among them, TP represents true positive examples, that is, the number of grassland pixels correctly predicted by the model; TN represents true negative examples, that is, the number of non-grazing grassland pixels correctly predicted by the model; FP represents false negative examples, that is, the number of non-grazing grassland pixels mistakenly classified as grassland pixels by the model; FN represents false negative examples, that is, the number of grassland pixels mistakenly classified as non-grazing grassland pixels by the model.

[0029] (3) Model interpretability based on Gradient-weighted Class Activation Mapping (Gradient-CAM); Grad-CAM calculates the gradient of the target category for the output of a specific convolutional layer, averages and pools these gradients, and then multiplies them with the feature map of the convolutional layer to obtain a class activation map. 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 a visual interpretation for any network based on the CNN model. The present invention applies Grad-CAM to the visual interpretation of the LSTM-FCN model. For the input time series data, the moments that contribute most to the prediction result are highlighted, and the model decision process is explained in a qualitative way.

[0030] (4) Model reasoning and post-processing; Model reasoning refers to the use of optimal model parameters to predict the range of regional pastures based on time-series remote sensing images. The present invention performs time-series band synthesis on the NDVI images within the range of "T50TMP" and "T50TMQ". For example, the original NDVI image set at 81 moments is synthesized into a remote sensing image containing 81 bands according to the chronological order. According to the aforementioned accuracy comparison results on the test set, the best model is selected to calculate the spatial distribution map of pasture pixels in the study area. In order to enhance the classification effect and reduce noise, the present invention uses median filtering to post-process the prediction result map.

Claims

1. A method for spatial mapping of natural temperate grassland using NDVI time series, characterized in that: It includes 4 steps: time series generation, classification of grazing and non-grazing pastures, model interpretation and model inference; (1) Time series generation; Collect Harmonized Landsat Sentinel-2 data and HLS data for natural grassland detection; HLS data is a seamless surface reflectance database generated by fusing and coordinating the land imagers and multi-spectrometers carried on Landsat-8 / 9 and Sentinel-2A / B remote sensing satellites respectively; The remote sensing data obtained was first used to remove cloud pixels, cloud shadow pixels and cloud adjacent pixels using its FMask band; at the same time, the ESAWorldCover global land use dataset was downloaded. ESA WorldCover aims to provide high-resolution global land cover data; The data is based on Sentinel-1 and Sentinel-2 satellite data, providing land cover data information at a resolution of 10 meters, including 11 different land cover categories: forest land, grassland, shrubs, cultivated land, buildings, bare soil / sparse vegetation areas, ice and snow, water areas, wetlands, mangroves and mosses; grassland range is extracted; After obtaining HLS data of cloud-free and grassland areas, calculate NDVI: Subsequently, linear interpolation and Savitzky-Golay filtering were used to obtain temporally continuous and spatially seamless NDVI data; (2) classification of pastures into grazing and non-grazing areas; The training set sample blocks and validation set sample blocks were obtained through visual interpretation, and the test set sample blocks were obtained through field investigation. The sample blocks were then converted into raster sample points through the "Vector to Raster" tool. In the subsequent model training, validation, and prediction processes, the labels of grassland pixels were set to "1" and the labels of non-grazing grassland pixels were set to "0". Based on the above data set, the effects of four classification models in the classification of threshed pasture and non-threshed pasture were compared; Use cross entropy loss to guide model training; define loss function for: in, Indicates the number of training samples; Indicates the number of categories; Indicates n The true labels of samples, Indicates n The predicted labels of samples were obtained; the Adam optimizer was used in the model training process, and the initial learning rate was set to 0.001; a total of 100 rounds of model training were set; after each round of training, the model was verified using the validation set, and the classification accuracy on the validation set was calculated; finally, the model parameters with the highest accuracy on the validation set were used to perform classification accuracy tests on the test set and compare different models; The overall accuracy OA, F1 accuracy, recall rate recall and precision rate precision are used to evaluate the classification accuracy; (3) Model interpretability based on class activation gradient weighted mapping Grad-CAM; Grad-CAM calculates the gradient of the target category for the output of a specific convolutional layer, averages and pools these gradients, and then multiplies them with the feature map of the convolutional layer to obtain the class activation map. Grad-CAM is applied to the visual interpretation of the LSTM-FCN model. For the input time series data, the moments that contribute most to the prediction results are highlighted, and the model decision process is explained in a qualitative way. (4) Model reasoning and post-processing; Model reasoning refers to the use of optimal model parameters to predict the scope of regional pastures based on time-series remote sensing images; the NDVI images within the "T50TMP" and "T50TMQ" ranges are synthesized into time-series bands, and the optimal model is selected based on the accuracy comparison results on the test set to calculate the spatial distribution map of pasture pixels in the study area; in order to enhance the classification effect and reduce noise, the median filter is used to post-process the prediction result map.

2. The method for spatial mapping of natural temperate grassland using NDVI time series according to claim 1, characterized in that: In step (1), NDVI is calculated according to the formula: in, and They represent the near-infrared band and the red light band respectively.

3. The method for spatial mapping of natural temperate grassland using NDVI time series according to claim 1, characterized in that: In step (1), the specific method of using linear interpolation and Savitzky-Golay filtering to obtain temporally continuous and spatially seamless NDVI data is as follows; Apply linear interpolation method to obtain initial interpolation of missing NDVI values; define the time to be interpolated x m The NDVI is y m , define the two observation points closest to the data to be interpolated as ( x a , y a )and( x b , y b ),in x a , x b Indicates the time point, y a , y b Represents the observed NDVI value at the corresponding time; then the interpolation time x m The initial interpolation value is: After initially filling the missing values ​​in the NDVI sequence by linear interpolation, the initial time series NDVI was optimized using Savitzky-Golay filtering. The Savitzky-Golay filter retains the high-frequency characteristics of the data by fitting a polynomial of the data points, thereby reducing the NDVI data noise caused by cloud cover, shadows, or high atmospheric aerosol concentrations. The savgol_filter method in SciPy was used to generate the final NDVI sequence.

4. The method for spatial mapping of natural temperate grassland using NDVI time series according to claim 1, characterized in that: In step (2), the four classification models are: random forest, multi-layer perceptron, 1D-CNN and LSTM-FCN.

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

Images