Method for determining manual recovery priority of burned area of coniferous forest in cold temperate zone

Abstract

The invention discloses a manual recovery priority determination method for a coniferous forest burned area in a cold temperate zone, and belongs to the technical field of remote sensing. The method solves the problem that when coniferous forests in different sub-regions in a burned region are recovered by an existing method, the natural updating capacity difference of the coniferous forests in the sub-regions is not considered. The method comprises the following steps: firstly, calculating a difference value normalized burning index of each pixel in a target area of a post-fire image, then determining an initial natural updating capability value according to the difference value normalized burning index, then modifying the natural updating capability of each pixel according to the altitude, the gradient and the slope direction value, and further determining the type of each pixel. And finally, according to the type of the plaque and the size of the area of the plaque, determining the manual recovery priority of each plaque. The method can be applied to manual recovery of the coniferous forest burned area in the cold temperate zone.

Description

Technical Field

[0001] The invention belongs to the field of remote sensing technology, and in particular relates to a method for determining the priority of artificial restoration of a burnt area in a cold temperate coniferous forest. Background Art

[0002] Due to global warming, the frequency of lightning-induced fires has increased year by year, especially in the northern part of Greater Khingan Range. Cold temperate coniferous forests are the main forest types in Greater Khingan Range, mainly including Larix gmelinii, Pinus pumila, Picea koraiensis and Pinus sylvestris, among which Pinus pumila accounts for more than 80% in high-altitude areas where lightning fires are frequent. Pinus pumila is a highly flammable tree species because it is rich in oil and has a clumping or clustered distribution. Therefore, once a fire occurs, cold temperate coniferous forests will suffer huge losses. In the early post-fire period, due to the influence of fire severity and terrain, there are significant differences in the natural regeneration capacity of vegetation in different regions. However, the traditional artificial restoration method of burned areas does not take into account the differences in the natural regeneration capacity of coniferous forests in each sub-region. The afforestation design and implementation are relatively extensive, requiring a lot of manpower, material and financial resources. By considering the natural regeneration capacity of different regions and giving priority to restoring areas with low natural regeneration capacity, we can reduce restoration costs while effectively protecting forest resources. Summary of the invention

[0003] The purpose of the present invention is to solve the problem that the existing methods do not take into account the differences in the natural regeneration capacity of coniferous forests in different sub-areas of a burnt area when restoring coniferous forests in different sub-areas of a burnt area, and to propose a method for determining the priority of artificial restoration of burnt areas of cold-temperate coniferous forests.

[0004] The technical solution adopted by the present invention to solve the above technical problems is: a method for determining the priority of artificial restoration of burnt areas in cold temperate coniferous forests, the method specifically comprising the following steps:

[0005] Step 1: determine the burned area to be restored, and obtain pre-fire remote sensing images and post-fire remote sensing images of the burned area to be restored;

[0006] Preprocessing the acquired remote sensing images respectively to obtain preprocessed remote sensing images before the fire and preprocessed remote sensing images after the fire;

[0007] Step 2: determine the burnt area boundary according to the pre-processed pre-fire remote sensing image and the pre-processed post-fire remote sensing image, and then cut out the pre-fire image and the post-fire image of the burnt area from the pre-processed pre-fire remote sensing image and the pre-processed post-fire remote sensing image according to the burnt area boundary;

[0008] Step 3: taking the coniferous forest area in the burned area as the target area for restoration, generating a vector boundary of the target area, and then cutting out the pre-fire image and the post-fire image of the target area from the pre-fire image and the post-fire image of the burned area according to the boundary of the target area;

[0009] Step 4: Calculate the difference normalized burn index of each pixel in the target area post-fire image according to the reflectivity of each pixel in the near-infrared band and the infrared band in the pre-fire image and post-fire image of the target area cropped in step 3, and assign a value to the natural renewal capacity of each pixel according to the difference normalized burn index of each pixel, that is, determine the initial natural renewal capacity value of each pixel in the target area post-fire image before artificial restoration;

[0010] Step 5: correcting the initial natural renewal capacity value of each pixel in the post-fire image of the target area according to the altitude, slope and aspect information of each pixel in the post-fire image of the target area to obtain a corrected natural renewal capacity value;

[0011] Step 6: Determine the type of each pixel in the post-fire image of the target area according to the corrected natural renewal ability value, fuse adjacent pixels of the same type into a patch, and calculate the area of ​​each fused patch;

[0012] Step 7: Determine the priority of artificial restoration of different patches in the target area according to the patch type and patch area.

[0013] Furthermore, the acquired remote sensing images are preprocessed respectively, and the preprocessing method is:

[0014] Step 1, performing histogram equalization and contrast enhancement on the remote sensing image in sequence to obtain a remote sensing image with enhanced contrast;

[0015] Step 1 and 2: The contrast-enhanced remote sensing image is then subjected to band combination to obtain a preprocessed remote sensing image.

[0016] Furthermore, the cloud coverage of the pre-fire remote sensing image and the post-fire remote sensing image of the burned area to be restored during collection needs to be less than a set cloud coverage threshold.

[0017] Furthermore, the boundary of the burned area is determined based on the pre-processed remote sensing image before the fire and the pre-processed remote sensing image after the fire, specifically:

[0018] Taking the pre-processed remote sensing image before the fire as a benchmark, the burned area was manually drawn on the pre-processed remote sensing image after the fire through ENVI software and visual interpretation to generate a vector file of the burned area, that is, to obtain the boundary of the burned area.

[0019] Furthermore, the specific process of step 4 is as follows:

[0020] Step 4: Calculate the normalized burn index (NBR) of the i-th pixel in the pre-fire image of the target area according to the reflectivity of the i-th pixel in the near-infrared band and the reflectivity of the infrared band in the pre-fire image of the target area. pre,i :

[0021] NBR pre,i =(ρ pre,i,nir -ρ pre,i,swir ) / (ρ pre,i,nir +ρ pre,i,swir )

[0022] Among them, ρ pre,i,nir is the reflectivity of the i-th pixel in the pre-fire image of the target area in the near-infrared band, ρ pre,i,swir is the reflectivity of the ith pixel in the pre-fire image of the target area in the infrared band;

[0023] According to the reflectivity of the i-th pixel in the post-fire image of the target area in the near-infrared band and the reflectivity of the infrared band, the normalized burn index NBR of the i-th pixel in the post-fire image of the target area is calculated. post,i :

[0024] NBR post,i =(ρ post,i,nir -ρ post,i,swir ) / (ρ post,i,nir +ρ post,i,swir )

[0025] Among them, ρ post,i,nir is the reflectivity of the i-th pixel in the post-fire image of the target area in the near-infrared band, ρ post,i,swir is the reflectivity of the ith pixel in the post-fire image of the target area in the infrared band;

[0026] According to NBR pre,i and NBR post,i Calculate the difference normalized burn index dNBR of the i-th pixel in the post-fire image of the target area i ;

[0027] dNBR i =NBR pre,i -NBR post,i

[0028] Step 42: Determine the initial natural renewal capacity value of each pixel in the post-fire image of the target area according to the difference normalized burn index of each pixel;

[0029] If dNBR i <0.1, then the initial natural renewal capacity value R of the i-th pixel iAssign a value of 4;

[0030] If 0.1≤dNBR i <0.3, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 3;

[0031] If 0.3≤dNBR i <0.6, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 2;

[0032] If dNBR i ≥0.6, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 1.

[0033] Furthermore, the specific process of step five is as follows:

[0034] Step 51: Extract the altitude, slope and aspect information of each pixel in the post-fire image of the target area, and record the altitude information of the i-th pixel extracted as E i , the slope information of the extracted i-th pixel is recorded as S i , the extracted slope information of the i-th pixel is recorded as A i ;

[0035] Step 52: According to the altitude, slope and aspect information of each pixel, the initial natural renewal capacity value of each pixel in the post-fire image of the target area is corrected, specifically:

[0036] Step 521: Based on the altitude information E of the i-th pixel i The initial natural renewal capacity value R i To make a correction:

[0037] If the elevation information E of the i-th pixel i In the medium altitude area, the natural renewal capacity value of the i-th pixel is modified to R i ′=R i +1;

[0038] If the elevation information E of the i-th pixel i If the pixel is in a low altitude area or a high altitude area, the natural renewal capacity value of the i-th pixel is modified to R i ′=R i ;

[0039] Step 522: The natural update capability value R is calculated based on the slope information of the i-th pixel. i ' to correct:

[0040] If the slope information S of the i-th pixel i Satisfy: 6°≤S i<25°, the natural renewal capability value of the i-th pixel is modified to R i =R i ′+1;

[0041] If the slope information S of the i-th pixel i Satisfy: 25°≤S i or S i <6°, then the natural renewal capability value of the i-th pixel is modified to R i =R i ′;

[0042] Step 523: According to the slope information of the i-th pixel, the natural update capacity value R i ″To make correction:

[0043] If the slope information A of the i-th pixel i Satisfy: 112.5°≤A i <247.5°, then the natural renewal capability value of the i-th pixel is modified to R i ″′=R i ″+1;

[0044] If 0≤A i <112.5°or 247.5°≤A i <360°, the natural renewal capability value of the i-th pixel is modified to R i ″′=R i ″;

[0045] Step 53: R i ″′ is taken as the corrected natural update capability value of the i-th pixel.

[0046] Furthermore, the specific process of step six is ​​as follows:

[0047] Step 6. Obtain the median R of the corrected natural update capability value of each pixel median ;

[0048] If R i ″′>R median , then set the type of the i-th pixel to A;

[0049] If R i ″′<R median , then set the type of the i-th pixel to B;

[0050] If R i ″′=R median , then the type of the i-th pixel is C;

[0051] Step 6.2: Merge adjacent pixels of the same type into one patch, and then calculate the area of ​​each patch separately.j , j = 1, 2,…, J, J represents the number of patches after fusion.

[0052] Furthermore, the specific process of step seven is as follows:

[0053] (1) For plaque type B:

[0054] If the area of ​​the jth patch of type B is Area j ≥1ha, ha means hectare, then the priority of artificial restoration of the jth type B patch is high;

[0055] If the area of ​​the jth type B patch satisfies: 0.1ha≤Area j <1ha, then the priority of artificial restoration of the jth type B patch is low;

[0056] If the area of ​​the jth patch of type B is j <0.1ha, it is determined that the jth type B patch does not need artificial restoration;

[0057] (2) For plaque type C:

[0058] If the area of ​​the jth patch of type C is Then the priority of artificial restoration of the jth type C plaque is medium, where n represents the number of type C plaques;

[0059] If the area of ​​the jth patch of type C is Then it is determined that the jth plaque of type C does not need to be manually restored;

[0060] (3) For plaque type A:

[0061] All type A plaques do not require artificial restoration.

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

[0063] The present invention first calculates the difference normalized burn index of each pixel in the target area of ​​the post-fire image, and then determines the initial natural regeneration capacity value according to the difference normalized burn index, and then modifies the natural regeneration capacity of each pixel according to the size of the altitude, slope and aspect value, and then determines the type of each pixel, and merges adjacent pixels of the same type into a patch, and finally, determines the artificial restoration priority of each patch according to the patch type and patch area. The method of the present invention can be adapted to local conditions and give priority to restoring areas with low natural regeneration capacity, which can not only reduce the workload of artificial restoration, but also effectively protect the cold temperate coniferous forest resources, and provide a reference for the post-fire restoration of cold temperate coniferous forests. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 It is a flow chart of a method for calculating the priority of artificial restoration of a burnt area of ​​a cold temperate coniferous forest according to the present invention;

[0065] Figure 2 It is the manual restoration priority distribution map of the target area. DETAILED DESCRIPTION

[0066] Specific implementation method 1: Combination Figure 1 The present embodiment describes a method for calculating the priority of artificial restoration of a burnt area of ​​a cold temperate coniferous forest, and the method specifically comprises the following steps:

[0067] Step 1: determine the burned area to be restored, and obtain pre-fire remote sensing images and post-fire remote sensing images of the burned area to be restored;

[0068] The remote sensing images taken the closest to the fire are selected as the post-fire remote sensing images, and the remote sensing images taken one year before the fire and in the same season as the fire are selected as the pre-fire remote sensing images. That is, the pre-fire remote sensing images and the post-fire remote sensing images should be remote sensing images of the burned area in the same season in different years to avoid the interference of seasonal factors on priority calculation;

[0069] Preprocessing the acquired remote sensing images respectively to obtain preprocessed remote sensing images before the fire and preprocessed remote sensing images after the fire;

[0070] Step 2: determine the burnt area boundary according to the pre-processed pre-fire remote sensing image and the pre-processed post-fire remote sensing image, and then cut out the pre-fire image and the post-fire image of the burnt area from the pre-processed pre-fire remote sensing image and the pre-processed post-fire remote sensing image according to the burnt area boundary;

[0071] Step 3: Determine the cold temperate coniferous forest area as the restoration target, take the coniferous forest area in the burned area as the target area for restoration, generate the vector boundary of the target area, and then cut out the pre-fire image and post-fire image of the target area from the pre-fire image and post-fire image of the burned area according to the boundary of the target area;

[0072] The method for generating the vector boundary of the target area is as follows: obtain the latest annual forest resource survey data covering the burned area to be restored, load the forest resource survey vector data in ArcGIS software, screen the obtained forest resource survey data, and screen the small classes in which coniferous trees account for more than 70% of the tree species composition, that is, obtain the boundary of the area where the main tree species is coniferous forest;

[0073] Step 4: Calculate the difference normalized burn ratio (dNBR) of each pixel in the target area post-fire image according to the reflectivity of each pixel in the near-infrared band and the infrared band in the pre-fire image and post-fire image of the target area cropped in step 3, and assign a value to the natural renewal capacity of each pixel according to the difference normalized burn ratio of each pixel, that is, determine the initial natural renewal capacity value of each pixel in the target area post-fire image before artificial restoration;

[0074] Step 5: According to the altitude, slope and aspect information of each pixel in the post-fire image of the target area (since the pixels in the pre-fire image and the post-fire image are corresponding, either the pre-fire image or the post-fire image can be used as a reference here), the initial natural renewal capacity value of each pixel in the post-fire image of the target area is corrected to obtain a corrected natural renewal capacity value;

[0075] Step 6: Determine the type of each pixel in the post-fire image of the target area according to the corrected natural renewal ability value, fuse adjacent pixels of the same type into a patch, and calculate the area of ​​each fused patch;

[0076] Step 7: Determine the priority of artificial restoration of different patches in the target area according to the patch type and patch area.

[0077] After a fire in a cold temperate coniferous forest, post-fire assessment and restoration is a complex process. If the entire burn site is surveyed manually, the workload will be very heavy. However, the present invention can quickly assess and formulate a burn site restoration plan based on the restoration priority of each sub-area. In the manual restoration process, each patch that needs to be restored manually is restored in order from high to low priority.

[0078] Specific implementation method 2: This implementation method is different from the specific implementation method 1 in that the acquired remote sensing images are preprocessed respectively, and the preprocessing method is:

[0079] Step 1, performing histogram equalization and contrast enhancement on the remote sensing image in sequence to obtain a remote sensing image with enhanced contrast;

[0080] Step 1 and 2: The contrast-enhanced remote sensing image is then subjected to band combination to obtain a preprocessed remote sensing image.

[0081] The remote sensing image acquired is Sentinel-2MSI L2A data. The processing methods are histogram equalization, contrast enhancement and band combination (that is, the visible light, near-infrared band (10m resolution) and short-wave infrared band (20m resolution) are combined, and the data is resampled to 10m resolution). The visual effect of the image can be improved through histogram equalization and contrast enhancement.

[0082] The other steps and parameters are the same as those in the first embodiment.

[0083] Specific implementation method three: This implementation method is different from specific implementation methods one or two in that the cloud coverage when collecting the pre-fire remote sensing images and the post-fire remote sensing images of the burned area to be restored needs to be less than a set cloud coverage threshold.

[0084] The present invention sets a cloud coverage threshold to limit the maximum cloud coverage when collecting remote sensing images, thereby ensuring that clear remote sensing images are screened.

[0085] The other steps and parameters are the same as those in the first or second embodiment.

[0086] Specific implementation method 4: This implementation method is different from the first to third specific implementation methods in that the boundary of the burned area is determined based on the pre-processed remote sensing image before the fire and the pre-processed remote sensing image after the fire, specifically:

[0087] Taking the pre-processed remote sensing image before the fire as a benchmark, the burned area was manually drawn on the pre-processed remote sensing image after the fire through ENVI software and visual interpretation to generate a vector file of the burned area, that is, to obtain the boundary of the burned area.

[0088] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0089] Specific implementation method 5: This implementation method is different from the specific implementation methods 1 to 4 in that the specific process of step 4 is as follows:

[0090] Step 4: Calculate the normalized burn index (NBR) of the i-th pixel in the pre-fire image of the target area according to the reflectivity of the i-th pixel in the near-infrared band and the reflectivity of the infrared band in the pre-fire image of the target area. pre,i :

[0091] NBR pre,i =(ρ pre,i,nir -ρ pre,i,swir ) / (ρ pre,i,nir +ρ pre,i,swir )

[0092] Among them, ρ pre,i,niris the reflectivity of the i-th pixel in the pre-fire image of the target area in the near-infrared band, ρ pre,i,swir is the reflectivity of the ith pixel in the pre-fire image of the target area in the infrared band;

[0093] According to the reflectivity of the i-th pixel in the post-fire image of the target area in the near-infrared band and the reflectivity of the infrared band, the normalized burn index NBR of the i-th pixel in the post-fire image of the target area is calculated. post,i :

[0094] NBR post,i =(ρ post,i,nir -ρ post,i,swir ) / (ρ post,i,nir +ρ post,i,swir )

[0095] Among them, ρ post,i,nir is the reflectivity of the i-th pixel in the post-fire image of the target area in the near-infrared band, ρ post,i,swir is the reflectivity of the ith pixel in the post-fire image of the target area in the infrared band;

[0096] According to NBR pre,i and NBR post,i Calculate the difference normalized burn index dNBR of the i-th pixel in the post-fire image of the target area i ;

[0097] dNBR i =NBR pre,i -NBR post,i

[0098] Step 42: Determine the initial natural renewal capacity value of each pixel in the post-fire image of the target area according to the difference normalized burn index of each pixel;

[0099] If dNBR i <0.1, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 4;

[0100] If 0.1≤dNBR i <0.3, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 3;

[0101] If 0.3≤dNBR i <0.6, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 2;

[0102] If dNBR i ≥0.6, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 1.

[0103] The present invention numbers the pixels in the pre-fire image of the target area to be restored and the post-fire image of the target area to be restored. The pixels at the same position need to have the same pixel number in the two images to facilitate subsequent calculations.

[0104] The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

[0105] Specific implementation method 6: This implementation method is different from the specific implementation methods 1 to 5 in that the specific process of step 5 is as follows:

[0106] Step 51: Using ArcGIS software, extract the elevation, slope and aspect information of each pixel in the post-fire image of the target area from the DEM data covering the target area, and record the elevation information of the i-th pixel extracted as E i , the slope information of the extracted i-th pixel is recorded as S i , the extracted slope information of the i-th pixel is recorded as A i ;

[0107] Step 52: Using the EVNI software, the initial natural renewal capacity value of each pixel in the post-fire image of the target area is corrected according to the altitude, slope and aspect information of each pixel. Specifically:

[0108] Step 521: Based on the altitude information E of the i-th pixel i The initial natural renewal capacity value R i To make a correction:

[0109] According to the minimum and maximum values ​​of the altitude information of each pixel in the target area, the altitude of the target area is divided into three parts, namely low altitude area, medium altitude area and high altitude area. Among them, the natural regeneration ability of vegetation in the medium altitude area is stronger.

[0110] If the altitude information E of the i-th pixel i In the medium altitude area, the natural renewal capacity value of the i-th pixel is modified to R i ′=R i +1;

[0111] If the altitude information E of the i-th pixel i If the pixel is in a low altitude area or a high altitude area, the natural renewal capacity value of the i-th pixel is modified to R i ′=R i ;

[0112] Step 522: The natural update capability value R is calculated based on the slope information of the i-th pixel. i ' to make corrections:

[0113] According to the slope classification standard, the natural regeneration capacity of vegetation on gentle slopes (6° to 25°) is higher than that on slight slopes (<6°) and steep slopes (≥25°);

[0114] If the slope information S of the i-th pixel i Satisfy: 6°≤S i <25°, the natural renewal capability value of the i-th pixel is modified to R i =R i ′+1;

[0115] If the slope information S of the i-th pixel i Satisfy: 25°≤S i or S i <6°, then the natural renewal capability value of the i-th pixel is modified to R i =R i ′;

[0116] Step 523: According to the slope information of the i-th pixel, the natural update capacity value R i ″To make correction:

[0117] According to the slope aspect classification standard, the natural regeneration capacity of vegetation on sunny slopes (112.5°~247.5°) is higher than that on shady slopes (0°~112.5° and 247.5°~360°);

[0118] If the slope information A of the i-th pixel i Satisfy: 112.5°≤A i <247.5°, then the natural renewal capability value of the i-th pixel is modified to R i ″′=R i ″+1;

[0119] If 0≤A i <112.5°or 247.5°≤A i <360°, the natural renewal capability value of the i-th pixel is modified to R i ″′=R i ″;

[0120] Step 53: R i ″′ is taken as the corrected natural update capability value of the i-th pixel.

[0121] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5.

[0122] Specific implementation method 7: This implementation method is different from one of the specific implementation methods 1 to 6 in that the specific process of step 6 is as follows:

[0123] Step 6. Obtain the median R of the corrected natural update capability value of each pixelmedian ;

[0124] If R i ″′>R median , then set the type of the i-th pixel to A;

[0125] If R i ″′<R median , then set the type of the i-th pixel to B;

[0126] If R i ″′=R median , then the type of the i-th pixel is C;

[0127] Step 6. Use the patch analyst module in ArcGIS software to merge adjacent pixels of the same type into one patch, and then calculate the area of ​​each patch separately. j (Arcgis software calculates the true area of ​​the patch based on the geographic coordinates), j = 1, 2, …, J, J represents the number of patches after fusion.

[0128] The other steps and parameters are the same as those in Specific Embodiments 1 to 6.

[0129] Step 62 is further explained as follows: For any pixel in the target area, assuming that the type of the pixel is type A, the pixels adjacent to the pixel (which can be above, below, left or right) and of type A belong to the same patch as the pixel, and so on, and then continue to expand outward from the boundary pixels of this patch until there are no pixels of type A in the pixels adjacent to the boundary pixels of this patch, and finally this patch is obtained. Similarly, for pixels of type B, the above method is also used to fuse pixels of type B into patches. Finally, the types of each patch are: a patch composed of pixels of type A, a patch composed of pixels of type B, and a patch composed of pixels of type C.

[0130] Specific implementation eight: This implementation is different from one of the specific implementations one to seven in that the specific process of step seven is as follows:

[0131] (1) For a patch of type B (i.e., all pixels in the patch are of type B):

[0132] If the area of ​​the jth patch of type B is j ≥1ha, ha means hectare, then the priority of artificial restoration of the jth type B patch is high;

[0133] If the area of ​​the jth type B patch satisfies: 0.1ha≤Area j<1ha, then the priority of artificial restoration of the jth type B patch is low;

[0134] If the area of ​​the jth patch of type B is Area j <0.1ha, it is determined that the jth type B patch does not need artificial restoration;

[0135] (2) For plaque type C:

[0136] If the area of ​​the jth type of patch is C Then the priority of artificial restoration of the jth type C plaque is medium, where n represents the number of type C plaques;

[0137] If the area of ​​the jth type of patch is C Then it is determined that the jth plaque of type C does not need to be manually restored;

[0138] (3) For plaque type A:

[0139] All type A plaques do not require artificial restoration.

[0140] The other steps and parameters are the same as those in Specific Embodiments 1 to 7.

[0141] Example

[0142] The present invention is further described in detail below with reference to the accompanying drawings and examples. The case was selected at 12:40 on June 23, 2017, when a lightning fire occurred in the No. 56 area of ​​the Daxigou Protection Area in Huzhong District, Daxinganling. The fire point was located at 123°23′25″E, 52°0′51″N. This fire is the largest lightning fire in the fire records of Huzhong Forestry Bureau from 2014 to 2024. The total area of ​​the fire is 212.4ha, and most of the fire area is mixed forest of larch and Pinus pumila, which is a typical cold temperate coniferous forest type. Data from the meteorological station closest to the fire showed that on the day of the fire, the fire temperature reached a maximum of 34.5°C, the average wind speed was 1.1m / s, the wind direction was northwest, and the gust speed reached 13.8m / s. The specific steps and processes are as follows:

[0143] (1) Remote sensing data acquisition

[0144] Sentinel-2 MSI images covering the burned area were downloaded from the official website of the European Space Agency (ESA). The images were taken on July 19, 2016 (before the fire) and July 19, 2017 (after the fire), with a cloud cover of 10% and a data level of L2A.

[0145] (2) Extraction of the burnt area

[0146] Using ENVI software, we manually drew the burned area on the 2017 fire image based on the 2016 pre-fire image to generate a vector file of the burn area. We cropped the two images according to the burn area boundary to obtain the two-phase multispectral images of the burn area.

[0147] (3) Recovery target area extraction

[0148] By superimposing the 2020 Class II forest resource survey data of Huzhong Forestry Bureau, it was found that the burned areas are all distributed in the cold temperate coniferous forest area, so the entire burned area is a target area for restoration.

[0149] (4) Assignment of initial natural regeneration capacity of the target area

[0150] Calculate the dNBR value of each pixel in the target area to restore, and classify it according to the threshold to obtain the initial natural renewal capacity R of each pixel. i .

[0151] (5) Correct and restore the natural regeneration capacity of the target area

[0152] The DEM data covering the restoration target area is obtained on the geospatial data cloud platform, and the elevation, slope and aspect information are extracted from the DEM data.

[0153] The altitude range of the restoration target area is between 750m and 1036m, and the altitude of the area is divided into three equal parts: low altitude area (750m-845m), medium altitude area (845m-940m) and high altitude area (940m-1036m). i In the medium altitude area (845m~940m), the natural renewal capacity of the ith pixel is R i +1, if E i In low or high altitude areas, the natural renewal capacity R of the i-th pixel i unchanged; the natural renewal capacity corrected according to the altitude information is recorded as R i ′;

[0154] The slope of the recovery target area ranges from 1° to 31°. If 6°≤S i <25°, then the natural update ability of the correction pixel i is R i ′+1, if S i <6° or S i ≥25, then the natural renewal capacity of the i-th pixel R i ′ remains unchanged; the natural renewal capacity corrected according to the slope information is recorded as R i ″;

[0155] The slope of the restoration target area ranges from 0° to 360°. If 112.5°≤A i <247.5°, then correct the natural renewal ability R of the i-th pixel i ″+1, if 0°≤A i <112.5°or 247.5°≤A i <360°, then the natural renewal capacity of the i-th pixel is R i "constant.

[0156] (6) Pixel Fusion

[0157] By calculation, the natural renewal capacity of each pixel in the target area image is assigned a value between 1 and 7, and the median is 4. Therefore, the type of pixels with a natural renewal capacity of 5 to 7 is A; the type of pixels with a natural renewal capacity of 4 is B; the type of pixels with a natural renewal capacity of 1 to 3 is C.

[0158] Using the patch analyst module in ArcGIS software, adjacent pixels of the same type are merged into one patch, and the area of ​​the fused patch is calculated. j .

[0159] (7) Priority ranking of manual restoration in the target area

[0160] According to the plaque type and plaque area, the artificial restoration priority of different patches in the target area is determined, and the following is obtained: Figure 2 The manual recovery priority distribution diagram of the recovery target area is shown.

[0161] Figure 2 In the figure, the red area is of type A, and Area j Areas composed of patches of ≥1ha, where the natural regeneration capacity of vegetation is poor and the area is large, and the priority of artificial restoration is high; orange areas are type B, and The natural regeneration capacity of vegetation in this area is average, but the area accounts for more than 50% of the total area of ​​type B patches, and the priority of artificial restoration is medium; the yellow area is type C, and 0.1ha≤Area j The area is composed of patches of less than 1ha. Although the natural regeneration ability of vegetation in this area is poor, the area is small, and the priority of artificial restoration is low. The green area is composed of patches with strong natural regeneration ability of vegetation or weak natural regeneration ability of vegetation but small area. This area does not require artificial restoration.

[0162] The above calculation examples of the present invention are only used to explain the calculation model and calculation process of the present invention in detail, and are not intended to limit the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still within the scope of protection of the present invention.

Claims

1. A method for determining the priority of artificial restoration of burnt areas in cold temperate coniferous forests, characterized in that: The method specifically comprises the following steps: Step 1: determine the burned area to be restored, and obtain pre-fire remote sensing images and post-fire remote sensing images of the burned area to be restored; Preprocessing the acquired remote sensing images respectively to obtain preprocessed remote sensing images before the fire and preprocessed remote sensing images after the fire; Step 2: determine the burnt area boundary according to the pre-processed pre-fire remote sensing image and the pre-processed post-fire remote sensing image, and then cut out the pre-fire image and the post-fire image of the burnt area from the pre-processed pre-fire remote sensing image and the pre-processed post-fire remote sensing image according to the burnt area boundary; Step 3: taking the coniferous forest area in the burned area as the target area for restoration, generating a vector boundary of the target area, and then cutting out the pre-fire image and the post-fire image of the target area from the pre-fire image and the post-fire image of the burned area according to the boundary of the target area; Step 4: Calculate the difference normalized burn index of each pixel in the target area post-fire image according to the reflectivity of each pixel in the near-infrared band and the infrared band in the pre-fire image and post-fire image of the target area cropped in step 3, and assign a value to the natural renewal capacity of each pixel according to the difference normalized burn index of each pixel, that is, determine the initial natural renewal capacity value of each pixel in the target area post-fire image before artificial restoration; Step 5: correcting the initial natural renewal capacity value of each pixel in the post-fire image of the target area according to the altitude, slope and aspect information of each pixel in the post-fire image of the target area to obtain a corrected natural renewal capacity value; Step 6: Determine the type of each pixel in the post-fire image of the target area according to the corrected natural renewal ability value, fuse adjacent pixels of the same type into a patch, and calculate the area of ​​each fused patch; Step 7: Determine the priority of artificial restoration of different patches in the target area according to the patch type and patch area.

2. The method for determining the priority of artificial restoration of burnt areas in cold temperate coniferous forests according to claim 1, characterized in that: The remote sensing images obtained are preprocessed respectively, and the preprocessing method is as follows: Step 1, performing histogram equalization and contrast enhancement on the remote sensing image in sequence to obtain a remote sensing image with enhanced contrast; Step 1 and 2: The contrast-enhanced remote sensing image is then subjected to band combination to obtain a preprocessed remote sensing image.

3. The method for determining the priority of artificial restoration of burnt areas in cold temperate coniferous forests according to claim 2, characterized in that: The cloud coverage of the pre-fire remote sensing image and the post-fire remote sensing image of the burned area to be restored when they are collected must be less than a set cloud coverage threshold.

4. The method for determining the priority of artificial restoration of burnt areas in cold temperate coniferous forests according to claim 3 is characterized in that: The method of determining the boundary of the burned area according to the pre-processed remote sensing image before the fire and the pre-processed remote sensing image after the fire is specifically as follows: Taking the pre-processed remote sensing image before the fire as a benchmark, the burned area was manually drawn on the pre-processed remote sensing image after the fire through ENVI software and visual interpretation to generate a vector file of the burned area, that is, to obtain the boundary of the burned area.

5. The method for determining the priority of artificial restoration of burnt areas in cold temperate coniferous forests according to claim 4, characterized in that: The specific process of step 4 is as follows: Step 4: Calculate the normalized burn index (NBR) of the i-th pixel in the pre-fire image of the target area according to the reflectivity of the i-th pixel in the near-infrared band and the reflectivity of the infrared band in the pre-fire image of the target area. pre,i : NBR pre,i =(ρ pre,i,nir -r pre,i,swir ) / (ρ pre,i,nir +r pre,i,swir ) Among them, ρ pre,i,nir is the reflectivity of the i-th pixel in the pre-fire image of the target area in the near-infrared band, ρ pre,i,swir is the reflectivity of the ith pixel in the pre-fire image of the target area in the infrared band; According to the reflectivity of the i-th pixel in the post-fire image of the target area in the near-infrared band and the reflectivity of the infrared band, the normalized burn index NBR of the i-th pixel in the post-fire image of the target area is calculated. post,i : NBR post,i =(ρ post,i,nir -r post,i,swir ) / (ρ post,i,nir +r post,i,swir ) Among them, ρ post,i,nir is the reflectivity of the i-th pixel in the post-fire image of the target area in the near-infrared band, ρ post,i,swir is the reflectivity of the ith pixel in the post-fire image of the target area in the infrared band; According to NBR pre,i and NBR post,i Calculate the difference normalized burn index dNBR of the i-th pixel in the post-fire image of the target area i ; dNBR i =NBR pre,i -NBR post,i Step 42: Determine the initial natural renewal capacity value of each pixel in the post-fire image of the target area according to the difference normalized burn index of each pixel; If dNBR i <0.1, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 4; If 0.1≤dNBR i <0.3, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 3; If 0.3≤dNBR i <0.6, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 2; If dNBR i ≥0.6, then the initial natural renewal capacity value R of the i-th pixel i Assign a value of 1.

6. A method for determining the priority of artificial restoration of burnt areas in cold temperate coniferous forests according to claim 5, characterized in that: The specific process of step five is as follows: Step 51: Extract the altitude, slope and aspect information of each pixel in the post-fire image of the target area, and record the altitude information of the i-th pixel extracted as E i , the slope information of the extracted i-th pixel is recorded as S i , the extracted slope information of the i-th pixel is recorded as A i ; Step 52: According to the altitude, slope and aspect information of each pixel, the initial natural renewal capacity value of each pixel in the post-fire image of the target area is corrected, specifically: Step 521: Based on the altitude information E of the i-th pixel i The initial natural renewal capacity value R i To make a correction: If the elevation information E of the i-th pixel i In the medium altitude area, the natural renewal capacity value of the i-th pixel is modified to R i ′=R i +1; If the elevation information E of the i-th pixel i If the pixel is in a low altitude area or a high altitude area, the natural renewal capacity value of the i-th pixel is modified to R i ′=R i ; Step 522: The natural update capability value R is calculated based on the slope information of the i-th pixel. i ' to correct: If the slope information S of the i-th pixel i Satisfy: 6°≤S i <25°, the natural renewal capability value of the i-th pixel is modified to R i =R i ′+1; If the slope information S of the i-th pixel i Satisfy: 25°≤S i or S i <6°, then the natural renewal capability value of the i-th pixel is modified to R i =R i ′; Step 523: According to the slope information of the i-th pixel, the natural update capacity value R i ″To make correction: If the slope information A of the i-th pixel i Satisfy: 112.5°≤A i <247.5°, then the natural renewal capability value of the i-th pixel is modified to R i ″′=R i ″+1; If 0≤A i <112.5°or 247.5°≤A i <360°, the natural renewal capability value of the i-th pixel is modified to R i ″′=R i ″; Step 53: R i ″′ is taken as the corrected natural update capability value of the i-th pixel.

7. A method for determining the priority of artificial restoration of burnt areas in cold temperate coniferous forests according to claim 6, characterized in that: The specific process of step six is ​​as follows: Step 6. Obtain the median R of the corrected natural update capability value of each pixel median ; If R i ″′>R median , then set the type of the i-th pixel to A; If R i ″′<R median , then set the type of the i-th pixel to B; If R i ″′=R median , then the type of the i-th pixel is C; Step 6.2: Merge adjacent pixels of the same type into one patch, and then calculate the area of ​​each patch separately. j , j = 1, 2,…, J, J represents the number of patches after fusion.

8. The method for determining the priority of artificial restoration of burnt areas in cold-temperate coniferous forests according to claim 7, characterized in that: The specific process of step seven is as follows: (1) For plaque type B: If the area of ​​the jth patch of type B is Area j ≥1ha, ha means hectare, then the priority of artificial restoration of the jth type B patch is high; If the area of ​​the jth type B patch satisfies: 0.1ha≤Area j <1ha, then the priority of artificial restoration of the jth type B patch is low; If the area of ​​the jth patch of type B is Area j <0.1ha, it is determined that the jth type B patch does not need artificial restoration; (2) For plaque type C: If the area of ​​the jth type of patch is C Then the priority of artificial restoration of the jth type C plaque is medium, where n represents the number of type C plaques; If the area of ​​the jth type of patch is C Then it is determined that the jth type C plaque does not need to be manually restored; (3) For plaque type A: All type A plaques do not require artificial restoration.

Images