A method for quickly identifying human-disturbed land plots of water and soil loss based on multi-temporal sentinel-2

Abstract

The application provides a method for quickly identifying human disturbance plots of water and soil loss based on multi-temporal Sentinel-2, and is based on medium-high spatial resolution images and time series Sentinel-2 data, greatly improves the efficiency of spatio-temporal distribution monitoring of human disturbance regions and the accuracy of extracted results, and the model approximation after feature optimization and optimal parameter selection can not only ensure that the data is not over-fitted, but also can obtain a better performance prediction model; at the same time, combined with the random forest algorithm with good robustness and strong generalization ability, the bottleneck problems such as insufficient high spatio-temporal data, low extraction accuracy and few effective methods for efficiently and accurately extracting human disturbance regions faced by human disturbance region monitoring are solved. The above method realizes the fine extraction and identification of the spatio-temporal distribution of human disturbance regions under the condition of ensuring the accuracy.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of rapid, efficient and high-precision remote sensing of human disturbance area monitoring, and in particular to a method for quickly identifying human-disturbed land of water and soil loss based on multi-temporal Sentinel-2. BACKGROUND

[0002] How to fully exploit remote sensing data for efficient identification of human disturbance areas, and strengthen the supervision of production and construction projects, so as to reduce human-induced soil erosion, is a problem that needs to be solved at present.

[0003] Remote sensing technology provides a certain solution to the problem of human disturbance area extraction. Medium and high resolution remote sensing images are widely used in specific human disturbance area mapping due to their relatively high spatial resolution and rich texture information. However, as of now, the application of this aspect is limited to small-scale study areas and specific disturbance types. When faced with large-scale areas and multiple disturbance types, the amount of medium and high resolution image data is huge, and the combination of object-oriented methods also directly increases the time-consuming of image segmentation, sample training and type determination. It is not friendly to quickly extract human disturbance in large-scale areas, and it is difficult to meet the needs of efficient and real-time supervision of production and construction projects in actual work.

[0004] In contrast, medium and low resolution images have more advantages in rapid extraction of large-scale areas, but the conditions such as image temporal resolution and spatial resolution make it difficult to accurately complete the extraction of human disturbance. For example, the spatial resolution of MODIS is limited, which makes it difficult to meet the requirements of fine mapping in fragmented areas; although Landsat images have higher spatial resolution, their temporal resolution is low, making it difficult to obtain accurate disturbance information through complete time series data. The emergence of the Sentinel series of satellites of the European Space Agency provides an opportunity for large-scale, long-term, high-precision human disturbance monitoring. Sentinel-2 data has a multi-spectral band resolution of 10m and a dual-satellite 2-5d return period, providing similar spectral and spatial information as land satellites, which can greatly improve the operation efficiency while maintaining high extraction accuracy, effectively solving the problems encountered by medium and high resolution satellites, MODIS and Landsat in rapid extraction of human disturbance areas. At present, there are few application researches on human disturbance extraction based on Sentinel-2 time series images, and there is still a lack of efficient extraction methods for human disturbance with wide coverage, high precision and comprehensive types. Therefore, it is necessary to explore the mapping method and mapping potential of Sentinel-2 in this problem.

[0005] Therefore, the current human disturbance area identification and extraction is faced with the problems of relying on manual interpretation, less research on multi-type disturbance extraction method in large area, low operation efficiency of high-resolution image combined with object-oriented method, and the like. Based on Sentinel-2 time series data, by selecting the optimal classification feature combination and the optimal classification parameter, the random forest method is used to realize high-precision, low-cost and large-scale human disturbance area identification and extraction year by year, which can provide theoretical support and method reference for subsequent human disturbance area extraction and soil erosion related problems. SUMMARY

[0006] The present application aims to provide a human disturbance area spatio-temporal distribution information extraction method based on random forest and feature optimization, which can efficiently, quickly and accurately monitor human disturbance.

[0007] The method for quickly identifying human disturbance land of soil erosion based on multi-temporal Sentinel-2 provided by the present application comprises the following steps:

[0008] S1, obtaining single-temporal medium-high spatial resolution Sentinel-2 remote sensing images and performing orthorectification, radiation correction and geometric registration data preprocessing to obtain single-temporal medium-high spatial resolution Sentinel-2 remote sensing images containing a research area;

[0009] S2, forming L2A level multi-temporal time series Sentinel-2 images from multiple single-temporal medium-high spatial resolution Sentinel-2 remote sensing images containing the research area, and performing format conversion output and batch cutting processing on the L2A level multi-temporal time series Sentinel-2 images through the European Space Agency SNAP platform to obtain medium-high spatial resolution time series data of the research area;

[0010] S3, for the medium-high spatial resolution time series data of the research area obtained in step S2, calculating the normalized vegetation index, the normalized water index, the normalized building index, the ratio of residential area index, the NDVI difference between sparse and vigorous vegetation, the contrast, the variance and the entropy, thereby generating a medium-high spatial resolution index time series data set, and combining the red edge band R, the green edge band G, the blue edge band B and the near-infrared band Nir in the time series data to generate a medium-high spatial resolution candidate feature data set;

[0011] S4, according to the medium-high spatial resolution candidate feature data set obtained in step S3, combining sample points of each land object sampled in the field, constructing a classification sample library and performing stratified sampling, thereby forming a training set and a validation set;

[0012] S5, forming training data and test data based on the training set, using the training data to first use the random forest algorithm for model training, and then determining the optimal model parameters and optimal feature combination through the performance of the test data;

[0013] S6, applying the optimal model parameters and optimal feature combination to the study area to realize the identification of vegetation, water body, cultivated land, impervious layer and human disturbance land in the study area, and complete the extraction of the land use type in the study area;

[0014] S7, using the verification set to evaluate the extraction result of step S6, if the accuracy evaluation result does not meet the preset expected accuracy, returning to step S5, if the accuracy evaluation result meets the preset expected accuracy, completing the land use type identification, extracting the human disturbance area, and finally completing the extraction of the spatio-temporal distribution of the human disturbance area.

[0015] Further, in step S1, when performing geometric registration, relying on open source data, manually selecting a plurality of control points for geometric correction, thereby realizing geometric registration.

[0016] Further, before the L2A level multi-temporal time series Sentinel-2 image is converted and output by the ESA SNAP platform, and batch cropping processing is performed, the image is also filtered based on cloud cover: the data is filtered by less than 5% of the whole cloud cover or less than 2% of the cloud cover in the study area; The multi-temporal time series Sentinel-2 image used is composed of four bands of blue, green, red and near-infrared, and the spatial resolution is 10 meters.

[0017] Further, the calculation formulas of the normalized vegetation index, the normalized water index, the normalized building index, the ratio residential index, the NDVI difference between sparse and vigorous vegetation, the contrast, the variance and the entropy are as follows:

[0018] Normalized vegetation index NDVI:

[0019] NDVI=(Nir-Red) / (Nir+Red)

[0020] Normalized water index NDWI:

[0021] NDWI=(Green-Nir) / (Green+Nir)

[0022] Normalized building index NDBI:

[0023] NDBI=(Mir-Nir) / (Mir+Nir)

[0024] Ratio residential index RRI:

[0025] RRI = Blue / Nir

[0026] dNDVI = NDVI

[0027] dNDVI = NDVI 旺盛期 -NDVI 稀疏期

[0028] Contrast:

[0029]

[0030] Variance:

[0031]

[0032] Entropy:

[0033]

[0034] In the formula, Nir, Red, Green, Mir, and Blue are the near-infrared band, the red band, the green band, the mid-infrared band, and the blue band of the Sentinel image respectively; i and j are the row and column coordinates of the pixel in the image; P(i, j) is the gray level joint probability matrix; and μ is the pixel mean.

[0035] Further, the step S5 specifically includes:

[0036] The 12 features in the medium-high spatial resolution candidate feature data set are optimized, the relative importance score of the candidate features is calculated by using the training data and the CART-based random forest method, and the 12 candidate features are sorted according to the relative importance; the features ranked in the top 1 / 3 according to the importance are selected as the basic features, and are combined to obtain a feature basic combination to be verified; one feature with higher importance is added to the previous combination according to the importance from large to small, and finally 9 groups of feature combinations to be verified are obtained; the out-of-bag data generated by the random forest algorithm is used to evaluate the OOB error of each group of feature combinations to be verified, and internal cross-validation is performed; the OOB error of each feature combination is obtained, and the combination with the smallest error is the preliminary optimal feature combination under the corresponding model parameters;

[0037] The trained random forest model under different model parameters is verified by using test data to verify the accuracy, and the optimal model parameters are obtained, and the preliminary optimal feature combination corresponding to the optimal model parameters is the final optimal feature combination.

[0038] Further, the model parameters include: the number of random trees, the number of leaf nodes, and the maximum leaf depth.

[0039] Furthermore, step S7 involves accuracy evaluation, specifically using four indicators: overall accuracy, user accuracy, producer accuracy, and Kappa coefficient. The specific calculation formulas are as follows:

[0040] accuracy=(TP+TN) / (TP+FN+FP+TN)

[0041] precision = TP / (TP + FP)

[0042] recall = TP / (TP + FN)

[0043] Kappa = (pp e ) / (1-p e )

[0044] in,

[0045] p e = (a1*b1+a2*b2+….+a i* b i ) / (n*n)

[0046] In the formula, TP represents the number of samples that were actually positive but were predicted as positive, FP represents the number of samples that were actually negative but were predicted as positive, TN represents the number of samples that were actually negative but were predicted as negative, and FN represents the number of samples that were actually positive but were predicted as negative; p o Represents the overall classification accuracy, a i b represents the number of actual samples of land cover type i. i represents the predicted value of the sample, n represents the total number of samples, and p represents the probability of correctly classifying all samples.

[0047] The beneficial effects of the technical solution provided by this invention are as follows: Based on medium-to-high spatial resolution imagery and time-series Sentinel-2 data, this invention significantly improves the efficiency of monitoring the spatiotemporal distribution of man-made disturbance areas and the accuracy of extraction results. The model approximation after feature optimization and optimal parameter selection not only ensures that the data is not overfitted but also yields a better-performing prediction model. Simultaneously, by combining the robust and highly generalizable random forest algorithm, it solves the bottleneck problems faced in monitoring man-made disturbance areas, such as insufficient high-spatiotemporal data, low extraction accuracy, and a lack of effective methods for efficiently and accurately extracting man-made disturbance areas. The above methods achieve refined extraction and identification of the spatiotemporal distribution of man-made disturbance areas while ensuring accuracy. Detailed Implementation

[0048] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments of the present invention will be further described below.

[0049] The method for quickly identifying human disturbance plots of water and soil loss based on multi-temporal Sentinel-2 of the application comprises the following steps:

[0050] S1, obtain single-temporal medium-high spatial resolution Sentinel-2 remote sensing images and perform orthorectification, radiation correction and geometric registration data preprocessing to obtain single-temporal medium-high spatial resolution Sentinel-2 remote sensing images containing a study area. The radiation correction is completed through the calibration coefficient of different sensors; the orthorectification is completed with the aid of global 30-meter digital elevation model data; and the geometric registration is completed manually by selecting ground control points based on Google Earth, Tianditu and other medium-high resolution non-offset auxiliary open source data.

[0051] S2, form multi-temporal time series Sentinel-2 images of L2A level from multiple said single-temporal medium-high spatial resolution Sentinel-2 remote sensing images containing a study area, and perform format conversion output and batch cropping processing on the multi-temporal time series Sentinel-2 images of L2A level through the European Space Agency SNAP platform to obtain medium-high spatial resolution time series data of the study area.

[0052] The format conversion output can output a tif file. The multi-temporal time series Sentinel-2 images refer to remote sensing images obtained in different months within a year period, and the phase range should consider the growth phenology time of different crop types in the study area. Before performing format conversion output and batch cropping processing on the multi-temporal time series Sentinel-2 images of L2A level through the European Space Agency SNAP platform, the images are also screened based on cloud cover: the data is screened by selecting the whole scene cloud cover less than 5% or the cloud cover of the study area less than 2%; the multi-temporal time series Sentinel-2 images used are composed of four bands of blue, green, red and near-infrared, and the spatial resolution is 10 meters.

[0053] S3, for the medium-high spatial resolution time series data of the study area obtained in step S2, calculate the normalized vegetation index, the normalized water body index, the normalized building index, the ratio of residential area index, the NDVI difference between sparse and vigorous periods of vegetation (March and August), the contrast, the variance and the entropy, so as to generate a medium-high spatial resolution index time series data set, and combine the red edge band R, the green edge band G, the blue edge band B and the near-infrared band Nir in the time series data, a total of 12 features, to generate a medium-high spatial resolution selected feature data set.

[0054] The calculation formulas of the normalized vegetation index, the normalized water body index, the normalized building index, the ratio of residential area index, the NDVI difference between sparse and vigorous periods of vegetation, the contrast, the variance and the entropy are as follows:

[0055] Normalized Difference Vegetation Index NDVI:

[0056] NDVI = (Nir-Red) / (Nir+Red)

[0057] Normalized Difference Water Index NDWI:

[0058] NDWI = (Green-Nir) / (Green+Nir)

[0059] Normalized Difference Built-up Index NDBI:

[0060] NDBI = (Mir-Nir) / (Mir+Nir)

[0061] Ratio Resident Index RRI:

[0062] RRI = Blue / Nir

[0063] Difference of NDVI between sparse and flourishing vegetation period (March and August) dNDVI:

[0064] dNDVI = NDVI 旺盛期 -NDVI 稀疏期

[0065] Contrast:

[0066]

[0067] Variance:

[0068]

[0069] Entropy:

[0070]

[0071] In the formula, Nir, Red, Green, Mir and Blue are near-infrared band, red band, green band, mid-infrared band and blue band of the Sentinel image respectively; i and j are the row and column coordinates of the pixel in the image; P(i,j) is the gray level joint probability matrix; and μ is the pixel mean value.

[0072] S4, according to the medium-high spatial resolution feature data set obtained in step S3, a classification sample library is constructed in combination with sample points of each land object sampled in the field, and stratified sampling is performed according to a certain proportion, so as to form a training set and a verification set.

[0073] S5, training data and test data are formed based on the training set, the training data are used to perform model training by using a random forest algorithm, and then optimal model parameters and optimal feature combinations are determined by the performance of the test data.

[0074] Among them, based on the training set, 70% is extracted as training data with replacement, and 30% is extracted as test data.

[0075] The 12 features in the medium-high spatial resolution candidate feature data set are optimized, the relative importance score of the candidate features is calculated using the CART-based random forest method with the training data, and the 12 candidate features are sorted according to the relative importance; The features ranked in the top 1 / 3 according to the importance are selected as the basic features, and are combined to obtain the feature combination to be verified; According to the importance from large to small, one feature with higher importance is added on the basis of the previous combination, and finally 9 groups of feature combinations to be verified are obtained; The out-of-bag (OOB) data generated by the random forest algorithm, that is, the training data remaining after Bootstrap sampling, is used to evaluate the OOB error of each feature combination to be verified, and internal cross-validation is performed; The OOB error of each feature combination is obtained, and the combination with the smallest error is the preliminary optimal feature combination under the corresponding model parameters;

[0076] The trained random forest model under different model parameters is verified by the test data to obtain the optimal model parameters, and the preliminary optimal feature combination corresponding to the optimal model parameters is the final optimal feature combination. The model parameters include: the number of random trees, the number of leaf nodes, and the maximum leaf depth.

[0077] S6, apply the optimal model parameters and the optimal feature combination to the study area to realize the identification of vegetation, water body, cultivated land, impervious layer and human disturbance land in the study area, and complete the extraction of the land use types (vegetation, water body, cultivated land, impervious layer and human disturbance land) in the study area.

[0078] S7, use the verification set to evaluate the precision of the extraction result of step S6. If the precision evaluation result does not meet the preset expected precision, return to step S5. If the precision evaluation result meets the preset expected precision, complete the land use type identification and extract the human disturbance area, and finally complete the extraction of the spatiotemporal distribution of the human disturbance area.

[0079] The accuracy, user precision, producer precision and Kappa coefficient are used to evaluate the accuracy, user precision, producer precision and Kappa coefficient. When the four indexes meet the preset expected precision, the land use type identification is completed. The specific calculation formula is as follows:

[0080] accuracy=(TP+TN) / (TP+FN+FP+TN)

[0081] precision = TP / (TP + FP)

[0082] recall = TP / (TP + FN)

[0083] Kappa = (p - p e ) / (1 - p e )

[0084] wherein,

[0085] p e = (a1*b1 + a2*b2 + … + a i* b i ) / (n*n)

[0086] In the formula, TP is the number of samples actually being positive samples and predicted as positive samples, FP is the number of samples actually being negative samples and predicted as positive samples, TN is the number of samples actually being negative samples and predicted as negative samples, and FN is the number of samples actually being positive samples and predicted as negative samples; p o represents the overall classification accuracy, a i represents the number of real samples of i type ground objects, b i represents the predicted value of the sample, n represents the total number of samples, and p represents the probability of correct classification of all samples.

[0087] In this document, it should be understood that the use of orientation words should not limit the scope of the application claimed.

[0088] In the case of no conflict, the above embodiments and features in the embodiments can be combined with each other.

[0089] The above only describes the preferred embodiments of the present application and should not be used to limit the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application should be included in the protection scope of the present application.

Claims

1. A method for rapidly identifying anthropogenically disturbed plots of land affected by soil erosion based on multi-temporal Sentinel-2, characterized in that, Includes the following steps: S1. Acquire single-temporal high spatial resolution Sentinel-2 remote sensing images and perform orthorectification, radiometric correction and geometric registration data preprocessing to obtain single-temporal high spatial resolution Sentinel-2 remote sensing images containing the study area. S2. Multiple single-temporal, medium-high spatial resolution Sentinel-2 remote sensing images containing the study area are used to form L2A-level multi-temporal time-series Sentinel-2 images. The L2A-level multi-temporal time-series Sentinel-2 images are converted and output through the ESA SNAP platform and batch cropped to obtain medium-high spatial resolution time-series data of the study area. S3. For the medium-to-high spatial resolution time series data of the study area obtained in step S2, calculate the normalized vegetation index, normalized water index, normalized building index, ratio residential area index, NDVI difference between sparse and lush vegetation periods, contrast, variance, and entropy to generate a medium-to-high spatial resolution index time series dataset. Combine this with the red-edge band R, green-edge band G, blue-edge band B, and near-infrared band Nir in the time series data to generate a medium-to-high spatial resolution candidate feature dataset. S4. Based on the medium-high spatial resolution candidate feature dataset obtained in step S3, and combined with the sample points of various objects sampled in the field, construct a classification sample library and perform stratified sampling to form a training set and a validation set. S5. Based on the training set, training data and test data are formed. The model is first trained using the random forest algorithm using the training data, and then the optimal model parameters and optimal feature combination are determined by the performance of the test data. S6. Apply the optimal model parameters and optimal feature combination to the study area to identify vegetation, water bodies, cultivated land, impermeable layers and human-disturbed land parcels in the study area, and complete the extraction of land use types in the study area. S7. Use the validation set to evaluate the accuracy of the extraction results in step S6. If the accuracy evaluation result does not meet the preset expected accuracy, return to step S5. If the accuracy evaluation result meets the preset expected accuracy, complete the land use type identification, extract the human disturbance area, and finally complete the spatiotemporal distribution extraction of the human disturbance area. The random forest in step S5 specifically includes: Twelve features from the medium-high spatial resolution candidate feature dataset are optimized. The relative importance scores of the candidate features are calculated using the CART-based random forest method with the training data, and the twelve candidate features are ranked according to their relative importance. The features ranked in the top 1 / 3 according to their importance are selected as the basic features and combined to obtain the basic combination of features to be verified. Based on the importance of each feature combination from highest to lowest, a more important feature is added to each combination, resulting in nine feature combinations to be validated. Then, using out-of-bag (OOB) data generated by the random forest algorithm, the OOB error of each feature combination is evaluated, and internal cross-validation is performed. The OOB error of each feature combination is obtained, and the combination with the smallest error is the preliminary optimal feature combination under the corresponding model parameters. The accuracy of the trained random forest model under different model parameters is validated using test data to obtain the optimal model parameters. The preliminary optimal feature combination corresponding to the optimal model parameters is the final optimal feature combination.

2. The method for rapid identification of anthropogenically disturbed land parcels affected by soil erosion based on multi-temporal Sentinel-2 as described in claim 1, is characterized in that, In step S1, during geometric registration, several control points are manually selected based on open-source data for geometric correction, thereby achieving geometric registration.

3. The method for rapidly identifying anthropogenically disturbed land parcels prone to soil erosion based on multi-temporal Sentinel-2 as described in claim 1, characterized in that, Before converting and batch cropping the L2A-level multi-temporal time-series Sentinel-2 images through the ESA SNAP platform, the images were also screened based on cloud cover: the cloud cover of the entire scene was less than 5% or the cloud cover of the study area was less than 2%. The multi-temporal time-series Sentinel-2 images used consist of four bands: blue, green, red, and near-infrared, with a spatial resolution of 10 meters.

4. The method for rapidly identifying anthropogenically disturbed land parcels prone to soil erosion based on multi-temporal Sentinel-2 as described in claim 1, characterized in that, The formulas for calculating the Normalized Difference Vegetation Index (NDVI), Normalized Water Index (NDVI), Normalized Building Index (NDVI), Ratio Residential Index, NDVI Difference between Sparse and Abundant Vegetation Periods, Contrast Ratio, Variance, and Entropy are as follows: Normalized Difference Vegetation Index (NDVI): NDVI = (Nir - Red) / (Nir + Red) Normalized Difference Water Index (NDWI): NDWI=(Green-Nir) / (Green+Nir) Normalized Building Index (NDBI): NDBI = (Mir - Nir) / (Mir + Nir) Ratio Residential Index (RRI): RRI = Blue / Nir dNDVI difference between sparse and lush vegetation periods: dNDVI=NDVI 旺盛期 -NDVI 稀疏期 Contrast: Variance: Entropy: In the formula, Nir, Red, Green, Mir, and Blue represent the near-infrared, red, green, mid-infrared, and blue bands of the sentinel image, respectively; i and j are the row and column coordinates of the pixel in the image; P(i,j) is the gray-level joint probability matrix; and μ is the pixel mean.

5. The method for rapidly identifying anthropogenically disturbed land parcels prone to soil erosion based on multi-temporal Sentinel-2 as described in claim 1, characterized in that, The model parameters include: the number of random trees, the number of leaf nodes, and the maximum leaf depth.

6. The method for rapidly identifying anthropogenically disturbed land parcels prone to soil erosion based on multi-temporal Sentinel-2 as described in claim 1, characterized in that, Step S7 involves accuracy evaluation, specifically using four indicators: overall accuracy, user accuracy, producer accuracy, and Kappa coefficient. The specific calculation formulas are as follows: accuracy=(TP+TN) / (TP+FN+FP+TN) precision = TP / (TP + FP) recall = TP / (TP + FN) Kappa=(pp e ) / (1-p e ) in, p e =(a1*b1+a2*b2+….+ai*bi) / (n*n) In the formula, TP represents the number of samples that were actually positive but were predicted as positive, FP represents the number of samples that were actually negative but were predicted as positive, TN represents the number of samples that were actually negative but were predicted as negative, and FN represents the number of samples that were actually positive but were predicted as negative; p o Represents the overall classification accuracy, a i b represents the number of actual samples of land cover type i. i represents the predicted value of the sample, n represents the total number of samples, and p represents the probability of correctly classifying all samples.

Images