Fire passing area extraction method fusing time sequence spectral features and multi-source fire products

Abstract

The invention relates to the technical field of an over-fire area extraction method, in particular to an over-fire area extraction method fusing time sequence spectral features and a multi-source fire product, which comprises the following steps of: S1, acquiring Sentinel-2 data, preprocessing and constructing a median image; s2, constructing a sample library of the overfire area and a sample library of the non-overfire area; s3, combining the training sample and the median image to screen and identify a spectral feature set of a fire passing area, and extracting an initial fire passing area; s4, the initial fire passing area is corrected by combining the multi-source fire point product and the enhanced combustion index, and a precise fire passing area is obtained; s5, merging the MCD64A1 product and the precise fire passing area, creating an envelope rectangle, and performing threshold segmentation on the month-by-month median image of the optimal spectral index by adopting an Otsu algorithm to obtain final fire passing area data; and S6, evaluating the precision of the fire passing area by using the verification sample, and drawing. Based on Sentinel-2 images and multi-source fire products, the invention provides a method for quickly identifying a large-scale and long-time-sequence fine fire passing area.

Description

Technical Field

[0001] The present invention relates to the technical field of burned area extraction methods, and particularly relates to a burned area extraction method that integrates temporal spectral features and multi-source fire products. Background Art

[0002] Forest fires are important disturbance factors in forest ecosystems, which can change the community structure, species composition, and ecosystem functions of forests, and thus have important impacts on aspects such as global climate change, carbon cycling, and biodiversity. Remote sensing satellites, with their advantages of high efficiency, accuracy, comprehensiveness, and near real-time nature, can accurately locate the position and scope of burned areas and quickly realize post-fire loss assessment and vegetation recovery monitoring. Currently, large-scale burned area products are mainly based on coarse-resolution optical satellite images such as MODIS. Limited by the spatial resolution, it is often difficult for these products to accurately identify small-area burned areas. The Sentinel-2 satellite has high spatial resolution, temporal resolution, and spectral resolution, and has the advantages of capturing the characteristics of burned areas and providing spatial context information. It is currently the best optical image data source for producing large-scale fine burned area products.

[0003] Currently, large-scale burned area detection mainly includes three steps: feature-based change detection, burned area detail correction, and post-processing. The accuracy of burned area products is improved through layer-by-layer detection and correction. Feature-based change detection includes generating an initial burned area based on spectral index thresholds and generating pixel-level combustion probabilities using methods such as random forest regression. Since the data produced in this step is vulnerable to the influence of cloud shadows and atmospheric correction errors, this data still needs to be further processed. Burned area detail correction reduces misclassification errors caused by non-fire-related spectral changes through temporal consistency checks, including validating with published fire point products and validating using feature indices to set detection criteria. The post-processing step mainly improves the integrity of burned area detection by using algorithms such as region growing, thereby reducing the missed classification error of burned areas. Currently, the methods for identifying burned areas face many challenges. First, due to the large dimensionality of data features and the complexity of image processing, the detection accuracy of burned areas is not high. Second, burned areas often have weak features and are difficult to be accurately identified and detected. Therefore, it is necessary to study how to select representative features and balance the detection accuracy of burned areas and the overall working efficiency. Summary of the Invention

[0004] In view of the above problems, the present invention provides a method for extracting burned areas by fusing temporal spectral features and multi-source fire products, realizing the extraction of large-scale and long-temporal fine burned areas. Based on Sentinel-2 temporal data, the present invention uses single-temporal and double-temporal scales to fully exploit spectral features, which can effectively reduce the influence of background heterogeneity and spectral diversity in the area; uses multi-source fire point products and enhanced combustion index to reduce the misclassification error in the initial burned area; combines with the MCD64A1 product and uses the Otsu algorithm to reduce the missed classification error in the fine burned area. The relevant research methods of the present invention can provide reference for the production of burned area products at the national, continental or even global scales.

[0005] A method for extracting burned areas by fusing temporal spectral features and multi-source fire products according to the present invention includes the following steps:

[0006] Step S1: Acquisition, preprocessing of Sentinel-2 data and construction of median image;

[0007] Step S2: Construction of sample libraries for burned areas and non-burned areas;

[0008] Step S3: Combining training samples and median image to screen the spectral feature set for identifying burned areas and extracting the initial burned area;

[0009] Step S4: Combining multi-source fire point products and enhanced combustion index to correct the initial burned area to obtain the accurate burned area;

[0010] Step S5: Merging the MCD64A1 product and the accurate burned area, creating an envelope rectangle, and using the Otsu algorithm to perform threshold segmentation on the monthly median image of the best spectral index to obtain the final burned area data;

[0011] Step S6: Using verification samples to evaluate the accuracy of the burned area and making a map.

[0012] Preferably, the specific steps of step S1 are as follows:

[0013] Step S11: Obtain the atmospherically corrected Sentinel-2 MSI image based on the Google Earth Engine platform, screen the images with less cloud cover, use the quality control band (QA60) to mask the pixels affected by clouds and cloud shadows in each image, and perform band resampling;

[0014] Step S12: Use the median composite method and the nearest date image interpolation method to construct the Sentinel-2 monthly median composite image.

[0015] Preferably, the specific steps of step S2 are as follows: Based on the Sentinel-2 monthly median composite image, manually select the burned area samples and non-burned area samples, construct a sample library, and divide the training samples and validation samples.

[0016] Preferably, the specific steps of step S3 are as follows:

[0017] Step S31: Based on the monthly median composite image, construct multiple spectral features, including 10 MSI bands (B2 - B8, B 8A and B 11 -B 12 ), and spectral features related to NIR and SWIR. The detailed information of each spectral feature is shown in Table 1: Table 1 Spectral Features Based on Sentinel-2 Images

[0018] Step S32: Using the constructed spectral features, calculate the monthly image median, extract the feature values corresponding to the training samples (burned area and non-burned area), and calculate the feature values of the training samples for single-temporal (current month) and double-temporal (current month and previous month);

[0019] Step S33: The separability index (M b ) is used to distinguish the ability of the spectral features of the Sentinel-2 image to distinguish between the burned area and the non-burned area. This index determines the degree of separability between them by comparing the spectral differences between objects and their internal spectral similarities. The calculation formula is as follows: where μb and μub are the pixel average values of the double-temporal differences in the burned area and non-burned area, σb and σub are the corresponding standard deviations. The higher the M b value, the better the separability. Combining the spectral feature values of the training samples, take the average value to calculate the M b values for the burned area and non-burned area under different features;

[0020] Step S34: Random forest is an ensemble learning model improved and developed based on multiple decision trees, which can quantitatively measure the contribution degree of features to classification. Among them, the out-of-bag (OOB) method calculates the importance of features by randomly swapping the values of a certain feature in the out-of-bag dataset, re-predicting, and measuring the degree of decrease in accuracy;

[0021] Step S35: A random forest model with 500 decision trees was constructed. Each spectral feature was used as an independent variable, and whether it was a burned area was used as the dependent variable. 80% of the training samples in Step S2 were used for training, 20% of the training samples were used for validation, and ten-fold cross-validation was performed.

[0022] Step S36: Considering M b and the results of the random forest model comprehensively, the spectral features with higher M b in single-temporal and double-temporal phases and large feature contribution degrees were selected as the spectral feature set for extracting the burned area.

[0023] Step S37: The single-temporal and double-temporal segmentation thresholds of each spectral feature were determined in combination with the training samples. The single-temporal and double-temporal burned areas corresponding to each feature were extracted according to the thresholds. On this basis, the common area of the single and double-temporal phases of each feature was calculated month by month as the initial burned area.

[0024] Preferably, the specific steps of Step S4 are as follows:

[0025] Step S41: The fire point products detected by MODIS and VIIRS were used to verify the large-area initial burned area. When the fire point fell within the patch range, this patch was considered a burned area patch.

[0026] Step S42: The fine burned area was identified by constructing the Enhanced Burned area index (EBAI). The identification standard for the small-area burned area in the current month was as follows: compared with the previous and next months of the current year, the spectral curve fluctuated most significantly in this month, that is, when the spectral value of this month was the maximum or minimum value; compared with the same month of the previous year, the spectral value fluctuated the most in the current month; compared with the difference between the monthly images of the previous year, the difference between the spectral values of this month and the previous month in the current year was the most significant. By repeatedly analyzing the segmentation threshold of the burned area image to be verified each month, the threshold of EBAI was determined as m, and the burned area was extracted.

[0027] Step S43: The large-area burned areas and small-area burned areas that had been verified in each month were combined as the corrected accurate burned area.

[0028] Preferably, the specific steps of Step S5 are as follows:

[0029] The Otsu algorithm is an unsupervised and non-parametric automatic global threshold determination method. It realizes the binary classification of the image by statistically analyzing the histogram of the global image. The segmentation threshold of this algorithm is the threshold corresponding to the maximum inter-class variance after iteration, which is taken as the optimal threshold.

[0030] Step S52: Combine the MCD64A1 product and the data of the accurately burned area, create a 200m buffer for the combined result, and obtain the minimum bounding rectangle of each burned patch. Use the Otsu algorithm to perform threshold segmentation on the monthly median image of the optimal spectral index to obtain the final monthly burned area data.

[0031] Preferably, the specific steps of step S6 are as follows:

[0032] Step S61: Use verification samples (burned area and unburned area) to evaluate the accuracy of extracting the burned area from Sentinel-2 images, and calculate the omission error (OE) and commission error (CE) to evaluate the extraction accuracy of the burned area;

[0033] Step S62: Use ArcGIS software to create a map of the burned area.

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

[0035] Based on Sentinel-2 time series data, the present invention uses single-temporal and double-temporal scales to fully exploit spectral features, which can effectively reduce the influence of regional background heterogeneity and spectral diversity; uses multi-source fire point products and enhanced combustion index to reduce the commission error in the initial burned area; combines the MCD64A1 product and uses the Otsu algorithm to reduce the omission error in the fine burned area. The related research methods of the present invention can provide reference for the production of burned area products at the national, continental or even global scales. Description of the Drawings

[0036] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following will briefly introduce the drawings required for the description of the embodiments. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.

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

[0038] Figure 2 is the M b value and the random forest importance score for identifying the burned area by each spectral index in Embodiment 2 of the present invention.

[0039] Figure 3 is the spatial distribution map of the burned area in Embodiment 2 of the present invention. Detailed Embodiments

[0040] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the following will be combined with the embodiments of the present invention. Describe the technical solutions in the embodiments of the present invention clearly and completely. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0041] Embodiment 1

[0042] Refer to Figure 1 , an embodiment of the present invention provides a method for extracting burned areas by fusing temporal spectral features and multi-source fire products, including the following steps:

[0043] Step S1: Acquisition, preprocessing of Sentinel-2 data and construction of monthly median images. The specific steps of step S1 are as follows:

[0044] Step S11: Obtain the atmospherically corrected Sentinel-2 MSI images based on the Google Earth Engine platform, screen the images with less cloud cover, use the quality control band (QA60) to mask the pixels affected by clouds and cloud shadows in each image, and perform band resampling.

[0045] Step S12: Based on the Google Earth Engine platform, use the median composite method and the nearest date image interpolation method to construct the Sentinel-2 monthly median composite images.

[0046] Step S2: Perform difference calculation (i.e., subtract the image of the previous month from the image of the current month) on the monthly median composite Sentinel-2 images through the Google Earth Engine platform to obtain monthly difference images. To reduce the influence of factors such as terrain, corresponding non-burned area sample points are selected near each burned area sample point. Through visual discrimination, burned area samples and non-burned area samples are selected to form a sample library.

[0047] Step S3: Use the separability index (M b ) and the random forest model, combine the training samples and the median images to screen and identify the spectral feature set of the burned area, and extract the initial burned area. The specific steps of step S3 are as follows:

[0048] Step S31: Based on the monthly median composite images, construct multiple spectral features, including 10 MSI bands (B2 - B8, B 8A and B 11 - B 12 ) and the spectral features related to NIR and SWIR. The detailed information of each spectral feature is shown in Table 1: Table 1 Spectral Features Based on Sentinel-2 Images

[0049] Step S32: Using the constructed spectral features, calculate the median value of the monthly images, extract the eigenvalue corresponding to the training samples (burned area and non-burned area), and calculate the eigenvalue of the training samples for single-temporal (current month) and double-temporal (current month and previous month).

[0050] Step S33: The separability index (M b ) is used to distinguish the ability of the spectral features of Sentinel-2 images to distinguish between burned areas and non-burned areas. This index determines the degree of separability between them by comparing the spectral differences between ground objects and their internal spectral similarities. The calculation formula is as follows: where μb and μub are the pixel average values of the double-temporal differences in the burned area and non-burned area, and σb and σub are the corresponding standard deviations. When the value of M b is less than 1, the separability is poor; when the value of M b is greater than 1, the separability is good. The higher the value of M b , the better the separability between the burned area and the unburned area. Combining the spectral eigenvalue of the training samples, taking the average value to calculate the M b value of the burned area and non-burned area under different features.

[0051] Step S34: Random forest is an ensemble learning model improved and developed based on multiple decision trees, which can quantitatively measure the contribution degree of features to classification. Among them, the out-of-bag (OOB) method calculates the importance of features by randomly exchanging the value of a certain feature in the out-of-bag dataset, re-predicting, and measuring the degree of accuracy decline.

[0052] Step S35: Construct a random forest model with 500 decision trees. Use each spectral feature as an independent variable and whether it is a burned area as a dependent variable. Use 80% of the training samples in Step S2 for training and 20% of the training samples for verification, and perform ten-fold cross-validation.

[0053] Step S36: Considering the results of M b and the random forest model comprehensively, select the spectral features with higher M b in single-temporal and double-temporal and large feature contribution degree in the random forest model as the spectral feature set for extracting the burned area.

[0054] Step S37: Determine the single-temporal and double-temporal segmentation thresholds for each spectral feature in combination with the training samples. Extract the single-temporal and double-temporal burned areas corresponding to each feature according to the thresholds. On this basis, obtain the common area of the single and double temporal phases of each feature month by month as the initial burned area.

[0055] Step S4: Although multiple typical burned areas are integrated to determine the segmentation thresholds for each feature, due to the influence of terrain, vegetation coverage, and pixel detection algorithms, there are likely to be a large number of noises in the initial burned area, which in turn reduces the detection accuracy of the initial burned area. Therefore, in order to produce high-quality burned area products, it is necessary to further correct the initial burned area. To reduce the misclassification error, the initial burned area is corrected by combining multi-source fire point products and the enhanced combustion index to obtain the accurate burned area. The specific steps of the said Step S4 are as follows:

[0056] Step S41: At present, the official has released a variety of high-quality fire point products that can be obtained for free, such as the fire point products detected by MODIS and VIIRS. Limited by the satellite spatial resolution and revisit cycle, these fire points often coincide with large fires with relatively long durations. Therefore, the released fire point products can be used to verify large areas of the initial burned area. The standard for verifying large areas of the burned area based on the fire point products is that when a fire point falls within the patch range, this patch is considered a burned area patch, and the un-verified patches are mostly small-area burned areas.

[0057] Step S42: Limited by the spatial resolution and temporal resolution of the fire point products, it is often difficult to identify burned patches with short combustion periods and small areas. The spectral changes of vegetation caused by fires can be reflected in the pre-fire and post-fire images of the same year, as well as in the images of the current year (post-fire) and the same period of the previous year. Therefore, in order to better identify such small-area burned patches, the enhanced burned area index (EBAI) can be constructed to identify the fine burned area. The identification standard for small-area burned areas in the current month is as follows: compared with the previous month and the next month of the current year, the spectral curve fluctuates most significantly in this month, that is, when the spectral value of this month is the maximum or minimum value; compared with the same month of the previous year, the spectral value fluctuates the most in the current month; compared with the differences in monthly images of the previous year, the spectral value difference between this month and the previous month of the current year is the most significant. By repeatedly analyzing the segmentation thresholds of the burned area images to be verified each month, the threshold of EBAI is determined as m, and the burned area is extracted.

[0058] Step S43: Combine the verified large-area burned areas and verified small-area burned areas of each month as the corrected accurate burned area.

[0059] Step S5: Applying a fixed threshold on Sentinel-2 images tends to omit some slightly burned fire patches and split a single large burned area into multiple small patches due to the unevenness of the underlying surface. Therefore, there may be missing classification in the accurate burned area data in the above steps. Sentinel-2 imagery can detect more refined burned area information, while the coarser spatial resolution MODIS satellite also has certain advantages in detecting large-scale burned areas. The published MCD64A1 burned area product can be used as auxiliary data to improve the accuracy of burned area detection to reduce the missing classification error in local areas. Merge the MCD64A1 product and the accurate burned area, create an envelope rectangle, and use the Otsu algorithm to perform threshold segmentation on the monthly median image of the best spectral index to obtain the final burned area data. The specific steps of the said step S5 are as follows:

[0060] Step S51: The Otsu algorithm is an unsupervised, non-parametric automatic global threshold determination method that realizes binary classification of an image by statistically analyzing the global image histogram. The segmentation threshold of this algorithm is, after iteration, the threshold corresponding to the maximized between-class variance, which is taken as the optimal threshold;

[0061] Step S52: Merge the MCD64A1 product and the accurate burned area data corresponding to each month, create a 200m buffer for the merged result, and obtain the minimum envelope rectangle of each fire patch. Use the Otsu algorithm to perform threshold segmentation on the monthly median image of the best spectral index to obtain the final monthly burned area data.

[0062] Step S6: Use validation samples to evaluate the accuracy of the burned area and create a map. The specific steps of the said step S6 are as follows:

[0063] Step S61: Based on the burned area sample points and non-burned area sample points in the training samples, evaluate the extraction accuracy of the burned area by calculating the omission error (OE) and the commission error (CE).

[0064] Step S62: Use ArcGIS software to create a map of the burned area.

[0065] Example 2

[0066] Refer to Figures 2-3 , which is the spatial distribution information of forest burned areas in the Heilongjiang River Basin from 2020 to 2023 extracted by the method of Example 1, including the following steps:

[0067] Step S1: Remote sensing data acquisition, preprocessing, and median image construction. The specific steps of the said step S1 are as follows:

[0068] Step S11: To accurately extract the burned areas in the forest regions of the Heilongjiang River Basin, the study used the global 30-meter HGFC data to mask out non-forest areas. Based on the Google Earth Engine platform, Sentinel-2 MSI images of the forest range in the Heilongjiang River Basin from 2019 to 2023 were obtained. Images with cloud cover less than 30% were screened, and the quality control band (QA60) of the Sentinel-2 images was used to mask the pixels affected by clouds and cloud shadows in each image. The B5, B6, B7, B8A, B11, and B12 bands with a spatial resolution of 20 meters were resampled to 10 meters.

[0069] Step S12: Using the time synthesis function provided by the Google Earth Engine platform, the median of the observations in the monthly time series stack for each observation was calculated to construct a monthly dense time series stack of Sentinel-2 images. To improve the utilization rate of the Sentinel-2 images, the nearest-date images were used to fill in the missing values.

[0070] Step S2: Based on the Google Earth Engine platform, difference calculations were performed on the monthly median composite Sentinel-2 images (i.e., subtracting the image of the previous month from the image of the current month) to obtain monthly difference images. To reduce the influence of factors such as terrain, corresponding non-burned sample points were selected near each burned sample point. Through visual discrimination, the samples were distributed as evenly as possible. A total of 600 burned area samples and 600 adjacent non-burned area samples were selected to form a sample library. Among them, 300 burned area samples and 300 adjacent non-burned area samples were randomly selected as training samples, and the remaining data were used as validation samples.

[0071] Step S3: Using the separability index (M b ) and the random forest model, combined with the training samples and the median image, the spectral feature set for identifying burned areas was screened, and the initial burned areas were extracted. The specific steps of Step S3 are as follows:

[0072] Step S31: Based on the monthly median composite images, multiple spectral features were constructed, including 10 MSI bands (B2 - B8, B 8A and B 11 -B 12 ) and spectral features related to NIR and SWIR. The detailed information of each spectral feature is shown in Table 1: Table 1 Spectral Features Based on Sentinel-2 Images

[0073] Step S32: Using the constructed spectral features, calculate the median of the monthly images, extract the eigenvalue corresponding to the training samples (burned areas and non-burned areas), and calculate the eigenvalue of the training samples for single-temporal (current month) and double-temporal (current month and previous month).

[0074] Step S33: The separability index (M b ) is used to distinguish the ability of the spectral features of Sentinel-2 images to distinguish between burned areas and non-burned areas. This index determines the degree of separability between them by comparing the spectral differences between ground objects and their internal spectral similarities. Its calculation formula is as follows: where μb and μub are the pixel averages of the double-temporal differences in the burned area and non-burned area, and σb and σub are the corresponding standard deviations. When the value of M b is less than 1, the separability is poor; when the value of M b is greater than 1, the separability is good. The higher the value of M b , the better the separability between the burned area and the unburned area. Combining the spectral eigenvalues of the training samples, take the average value to calculate the M b value for the burned area and non-burned area under different features;

[0075] Step S34: Random forest is an ensemble learning model improved and developed based on multiple decision trees, which can quantitatively measure the contribution degree of features to classification. Among them, the out-of-bag (OOB) method calculates the importance of features by randomly exchanging the value of a certain feature in the out-of-bag dataset, re-predicting, and measuring the degree of decrease in accuracy. The random forest importance assessment based on OOB includes two indicators: "MeanDecreaseAccuracy" and "MeanDecreaseGini". When the value of "MeanDecreaseAccuracy" or "MeanDecreaseGini" of a certain variable is larger, it indicates that the importance of this variable is greater.

[0076] Step S35: A random forest model with 500 decision trees is constructed. Using each spectral feature as an independent variable and whether it is a burned area as a dependent variable, 80% of the training samples in Step S2 are used for training, 20% of the training samples are used for verification, and ten-fold cross-validation is performed.

[0077] Step S36: Considering the results of M b and the random forest model comprehensively, select the spectral features with higher M b in single-temporal and double-temporal and with a large feature contribution degree in the random forest model as the spectral feature set for extracting the burned area ( Figure 2) After selection, the spectral features are MIRBI, NBR2, and VI1112.

[0078] Step S37: Synthesize the monthly median images corresponding to MIRBI, NBR2, and VI1112 respectively. Perform threshold segmentation on the single-temporal and double-temporal images respectively. The thresholds for each feature are determined by visual interpretation in combination with typical burned areas in the training samples. The single-temporal thresholds corresponding to MIRBI, NBR2, and VI1112 are 1.62, 0.13, and 1.31 respectively; the double-temporal thresholds corresponding to MIRBI, NBR2, and VI1112 are 0.27, -0.03, and -0.12 respectively. Use the threshold method to obtain the detected single-temporal and double-temporal burned area ranges respectively, and take the common area as the initial burned area recognized by the Sentinel-2 image.

[0079] Step S4: To reduce the misclassification error in the initial burned area, combine multi-source fire point products and the enhanced combustion index to correct the initial burned area and obtain the accurate burned area. The specific steps of step S4 are as follows:

[0080] Step S41: Use MODIS and VIIRS fire point products. Preprocess the forest area fire point products using ArcGIS software, including operations such as projection conversion, removing non-forest fire points, and attribute segmentation, to generate the monthly fire point products of the Heilongjiang River Basin from 2020 to 2023. Verify the large-scale burned areas through the temporal and spatial attributes of the fire points.

[0081] Step S42: For the small-scale burned patches to be verified, select the MIRBI index with the best performance to construct the enhanced combustion index. Taking the small-scale burned area in 2020 as an example, synthesize the monthly median images of the best spectral features in 2020 and extract the image medians to form a set a1 of 12-month image medians. Synthesize the monthly median images in 2019 based on Sentinel-2 data to form a set a2 of 12-month image medians. Subtract the monthly images, that is, subtract the previous month's image median from the current month's image median, to obtain the set b1 of monthly image median differences in 2020 and the set b2 of monthly image median differences in 2019. Obtain the spectral curves c1 and c2 through a1 - a2 and b1 - b2. The specific formula for EBAI is EBAI = b1 + c1 + c2. The best threshold of EBAI is manually adjusted by comprehensively considering the extraction effects of burned areas in each month. Finally, determine the threshold m of EBAI to be -1.6. Extract all the burned patches with an EBAI threshold less than -1.6.

[0082] Step S43: Combine the large-scale burned areas and small-scale burned areas verified in each month as the corrected accurate burned area.

[0083] Step S5: To reduce the omission error, the MCD64A1 product and the accurately burned area are merged. A 200m buffer is created for the monthly burned area data after merging, and its minimum bounding rectangle is obtained. The Otsu algorithm is used to perform threshold segmentation on the monthly MIRBI median image to obtain the final burned area data.

[0084] Step S6: Use the validation samples to evaluate the accuracy of the burned area and create a burned area map. The specific steps of Step S6 are as follows:

[0085] Step S61: Based on the time attribute and spatial attribute in the validation samples, the omission error (OE) and commission error (CE) are calculated to be 0.08 and 0.10 respectively.

[0086] Step S62: Use ArcGIS software to create a burned area map ( Figure 3 ). The above has generally described the present invention in detail. However, based on the present invention, some modifications or improvements can be made, which are obvious to those of ordinary skill in the technical field. Therefore, modifications or improvements made without departing from the spirit and idea of the present invention are within the protection scope of the present invention.

Claims

1. A method for extracting the burned area by integrating time series spectral features and multi-source fire products, characterized in that: The steps include: Step S1: Sentinel-2 data acquisition, preprocessing and median image construction; Step S2: constructing a sample library of the burned area and the non-burned area; Step S3: combining the training samples and the median image to screen and identify the spectral feature set of the burned area, and extracting the initial burned area; Step S4: Merge the MCD64A1 product and the precise burned area to create an envelope rectangle, use the Otsu algorithm to perform threshold segmentation on the monthly median image of the best spectral index, and obtain the final burned area data; combine the multi-source fire point product and the enhanced combustion index to correct the initial burned area and obtain the precise burned area; Step S5: Merge the MCD64A1 product and the precise burned area, create an envelope rectangle, and use the Otsu algorithm to perform threshold segmentation on the monthly median image of the best spectral index to obtain the final burned area data; Step S6: Use the validation samples to assess the accuracy of the burned area and make a map.

2. The method for extracting the burned area by integrating the time series spectrum characteristics and multi-source fire products according to claim 1 is characterized in that: The specific steps of step S1 are: Step S11: Based on the Google Earth Engine platform, the atmospherically corrected Sentinel-2 MSI images are obtained, the images with smaller cloud cover are screened, and the pixels in each image that are disturbed by clouds and cloud shadows are masked using the quality control band (QA60), and band resampling is performed; Step S12: Use the median synthesis method and the nearest date image interpolation method to construct the Sentinel-2 monthly median synthesis image.

3. The method for extracting the burned area by integrating the time series spectrum characteristics and multi-source fire products according to claim 2 is characterized in that: The specific steps of step S2 are: Based on the monthly median composite images of Sentinel-2, samples of burned areas and non-burned areas were manually selected to build a sample library, and the samples were divided into training samples and validation samples.

4. The method for extracting the burned area by integrating the time series spectrum characteristics and multi-source fire products according to claim 3 is characterized in that: The specific steps of step S3 are: Step S31: Based on the monthly median composite image, multiple spectral features are constructed, including 10 MSI bands (B2-B8, B 8A and B 11 -B 12 ) and spectral features related to NIR and SWIR. Detailed information on each spectral feature is shown in Table 1: Table 1 Spectral characteristics based on Sentinel-2 images Step S32: using the constructed spectral features, calculating the median of the images for each month, extracting the feature values ​​corresponding to the training samples (burned areas and non-burned areas), and calculating the feature values ​​of the training samples in a single phase (current month) and in two phases (current month and previous month); Step S33: Separability index (M b ) is used to distinguish the ability of the spectral characteristics of Sentinel-2 images in burned areas and non-burned areas. The index determines the degree of separability between objects by comparing the spectral differences between them and the spectral similarities within them. The calculation formula is as follows: Among them, μb and μub are the pixel averages of the dual-phase difference between the burned area and the non-burned area, σb and σub are the corresponding standard deviations, and M b The higher the value, the better the separability. Combined with the spectral feature values ​​of the training samples, the average value is taken to calculate the M of the burned area and the non-burned area under different characteristics. b value; Step S34: Random forest is an integrated learning model improved and developed on the basis of multiple decision trees, which can quantitatively measure the contribution of features to classification. Among them, the out-of-bag data (OOB) method randomly exchanges the value of a feature in the out-of-bag data set, re-predicts, and calculates the importance of the feature by measuring the degree of accuracy reduction; Step S35: A random forest model with 500 decision trees was constructed, with each spectral feature as an independent variable and whether it was a burned area as a dependent variable. 80% of the training samples in step S2 were used for training, and 20% of the training samples were used for validation, and dry-fold cross validation was performed; Step S36: Comprehensively consider M b The results of the random forest model were used to select M in single-phase and double-phase b The spectral features with higher contribution are used as the spectral feature set for extracting the burned area; Step S37: Determine the segmentation thresholds of the single-phase and dual-phase of each spectral feature in combination with the training samples, extract the single-phase and dual-phase burned areas corresponding to each feature according to the thresholds, and on this basis, obtain the common areas of the single-phase and dual-phase of each feature as the initial burned areas on a monthly basis.

5. The method for extracting the burned area by integrating the time series spectrum characteristics and multi-source fire products according to claim 4 is characterized in that: The specific steps of step S4 are: Step S41: Use the fire point products detected by MODIS and VIIRS to verify the large-scale initial burned area. When the fire point falls within the patch range, the patch is considered to be a burned area patch; Step S42: Identify fine burned areas by constructing an Enhanced Burned Area Index (EBAI). The identification criteria for small burned areas in the current month are: compared with the previous month and the next month of the current year, the spectral curve of this month fluctuates most significantly, that is, the spectral value of this month is the maximum value; compared with the same month of the previous year, the spectral value of the current month fluctuates the most; compared with the monthly image difference of the previous year, the spectral value difference between this month of the current year and the previous month is the most significant. By repeatedly analyzing the segmentation threshold of the burned area image to be verified each month, the threshold of EBAI is determined to be m, and the burned area is extracted; Step S43: Merge the verified large-area burned areas and verified small-area burned areas in each month as the corrected accurate burned areas.

6. The method for extracting the burned area by integrating the time series spectrum characteristics and multi-source fire products according to claim 5 is characterized in that: The specific steps of step S5 are: Step S51: Otsu algorithm is an unsupervised, non-parametric automatic global threshold determination method, which realizes the binary classification of the image by statistically analyzing the histogram of the global image. The segmentation threshold of the algorithm is the threshold corresponding to the maximum inter-class variance after iteration, which is taken as the optimal threshold; Step S52: Merge the MCD64A1 product and the precise burned area data, make a 200m buffer zone for the merged result, and obtain the minimum envelope rectangle of each burned patch. Use the Otsu algorithm to perform threshold segmentation on the monthly median image of the best spectral index to obtain the final monthly burned area data.

7. The method for extracting the burned area by integrating the time series spectrum characteristics and multi-source fire products according to claim 6 is characterized in that: The specific steps of step S6 are: Step S61: Use verification samples (burned areas and non-burned areas) to evaluate the accuracy of extracting burned areas from Sentinel-2 images, and calculate the omission error (OE) and commission error (CE) to evaluate the accuracy of extracting burned areas; Step S62: Mapping the fire area using ArcGIS software.

Images