A method for updating the forest stand age of forest subcompartments based on multi-source remote sensing data

Abstract

The present invention relates to a method for updating the forest stand age of forest sub-compartments based on multi-source remote sensing data. Based on time-series NDVI, random layout segmentation is used to find forest change points (time, location, quantity) to determine the stand age distribution; then, the sub-compartments are segmented based on multi-source high-resolution images, and finally, the stand age extracted from the time series is updated to the re-segmented sub-compartments. This method saves the labor and time costs of stand age surveys and provides certain support for forest stand age updating.

Description

Technical Field

[0001] The present invention belongs to the research field of forest remote sensing, and particularly relates to a method for updating the forest stand age of forest sub-compartments based on multi-source remote sensing data. Background Art

[0002] China has established a three-level forest resource inventory system: the National Forest Continuously Inventory (NFCI), the Forest Management Planning Inventory (FMPI), and the Forest Operational Design Inventory (FODI). Among them, FODI is a detailed survey of the smallest forest management unit (FSC) in China and is the only spatial data with forest stand age obtained based on the survey. However, the survey is conducted every five years, so the stand age information is relatively lagging.

[0003] Currently, there are mainly two ways to extract the stand age based on remote sensing data. One is to construct a stand age estimation model based on single-phase or multi-phase remote sensing data (such as tree height) and combine ground survey, meteorological and other data to extract the single-phase stand age spatial distribution at the national or global scale. The other is to extract the forest disturbance years through change detection based on time series remote sensing data, and then calculate the forest age. These remote sensing-based stand age products can, to a certain extent, serve well the estimation of forest stand parameters in the forest growth model.

[0004] However, due to the limitations of remote sensing data itself and algorithms, such as cloud and snow, noise, spatio-temporal resolution, breakpoint detection accuracy and other factors, there are certain errors and uncertainties in the obtained stand age information. For example, global segmentation is a commonly used method for detecting multiple change points in time series. It searches for the first change point from the entire dataset. After detecting the first breakpoint, the dataset is divided into two sub-segments according to the detected breakpoint. Then, a similar process is performed on each sub-segment. This binary segmentation strategy is widely applied to time series remote sensing data for change detection in different ecosystems.

[0005] However, according to the analysis of existing research results, it is found that the previous methods for obtaining stand age information have the following problems: 1) In the research of estimating stand age with tree height, the stand age is easily affected by the extraction accuracy of tree height, and there is a lack of other auxiliary information for reference; 2) In the research of extracting stand age by segmenting time series, there are problems such as misjudgment of change points and inability to effectively predict the starting and ending points; 3) The extracted stand age information is not further updated in combination with the sub-compartment boundary. Summary of the Invention

[0006] The purpose of the present invention is to provide a method for updating the forest stand age of forest sub-compartments based on multi-source remote sensing data, which can effectively improve the efficiency of updating the forest stand age of forest sub-compartments and reduce the cost of obtaining the sub-compartment stand age and rotation period. This method can be applied to planted forests such as eucalyptus, Chinese fir, and rubber.

[0007] To achieve the above object, the technical solution of the present invention is: A method for updating the forest stand age of forest sub-compartments based on multi-source remote sensing data, comprising the following steps:

[0008] Step 1: Prepare all available Landsat NDVI data from 2000 to 2021, GF-6 images and Sentinel-2 images in 2021, forest field survey data and forest sub-compartment data in 2021;

[0009] Step 2: Divide it into two parts: 1) Time series part: Randomly divide the NDVI time series of each pixel into M segments, and determine the maximum cumulative sum CUSUM of each segment, as well as the corresponding breakpoint quantity, position and order through sSIC calculation; 2) Multi-source high-resolution data part: Preprocess the GF-6 images and Sentinel-2 images respectively, and then perform data fusion;

[0010] Step 3: Divide it into two parts: 1) Time series part: Based on the change points and the characteristics of each segment, use a random forest classifier to classify the land use in different periods, obtain the spatial distribution map of eucalyptus, and determine the time, decade, rotation period and forest age of forest harvesting and stable recovery in 2021; 2) Multi-source high-resolution data part: Vectorize the sub-compartments in the forest survey data and perform multi-scale segmentation based on the fused high-resolution images to obtain the latest eucalyptus sub-compartments;

[0011] Step 4: Based on the results of Step 3, use the latest eucalyptus sub-compartments as input data, define the area, capture the forest stand age raster, and extract the pixel values of the forest stand age raster;

[0012] Step 5: Based on the results of Step 4, assign the extracted forest stand age pixel values to the sub-compartment raster;

[0013] Step 6: Convert the raster data after assignment into vector sub-compartments, that is, complete the update of the sub-compartment forest stand age.

[0014] In an embodiment of the present invention, in Step 1, select Landsat NDVI data with cloud cover below 80% in the study area from 2000 to 2021 on USGS, including TM, ETM+ SLC-off, ETM+ SLC-on and OLI; Remove the noise including clouds, shadows and stripes in each image according to the quality file information;

[0015] Collect land cover type samples through field survey records and visual interpretation on Google Earth; Set the classification system according to project requirements, and divide the land cover into 6 types: eucalyptus, non-eucalyptus forest, farmland, orchard, vegetation-free and logging slash;

[0016] The main information of the forest sub-compartments includes the average breast diameter, average tree height, age, survey date, and number of trees per hectare of the sub-compartments;

[0017] The multi-source high-resolution data includes GF-6 image data and Sentinel-2 images.

[0018] In an embodiment of the present invention, step 2 is specifically implemented as follows:

[0019] 1) Time series part: First, divide the NDVI time series into M segments; secondly, determine the corresponding breakpoint quantity, position and order;

[0020] a. Randomly extract subsequences, that is, the vector (NDVI s , NDVI s+1 ,..., NDVI e ), where s and e are the start and end observations of a subsequence, and calculate the CUSUM statistic for each subsequence; maximize each CUSUM, select the largest CUSUM among the M CUSUMs, and use it as the first candidate breakpoint to be tested according to the information criterion or a predefined threshold; once a breakpoint is determined to be significant, perform the same procedure on its left and right time series;

[0021]

[0022]

[0023] s < b < e, n = e - s + 1, b is the observation time in the middle of s and e, e m , s m , n m are the end position, start position and number of observations of the m-th segment, and NDVI(m, t) is the time series value of a pixel in the m-th sub-segment at time t; is the CUSUM statistic of b in the time series from s to e; f(m, b) is the sub-segment with the largest and its corresponding change point; is to take the m value and b value corresponding to the largest CUSUM;

[0024] b. Use sSIC to determine the position and quantity of breakpoints:

[0025]

[0026] α = 1 represents the penalty term of the standard Schwarz criterion; K is the number of candidate breakpoints; is the final variance of the maximum likelihood estimate of the corresponding residual; N is the number of observations in the time series data;

[0027] 2) Multi-source high-resolution data part: First, perform image preprocessing on GF-6 and Sentinel-2 respectively; then use the GS fusion algorithm to complete image fusion;

[0028] Radiometric correction and atmospheric correction are performed on GF-6; radiometric calibration and atmospheric correction are performed on Sentinel-2, and the multispectral bands are uniformly resampled to a resolution of 10 m.

[0029] Based on the GF-6 panchromatic image as a reference, geometric registration is performed on the Sentinel-2 multispectral image; the GS fusion algorithm is used to perform a GS transformation on the resampled multispectral image and the simulated panchromatic image, and the first band of the image after the GS transformation is replaced with the real panchromatic band, and finally a GS inverse transformation is performed to obtain a multispectral image with a resolution of 2 m.

[0030] In an embodiment of the present invention, the specific implementation of step 3 is as follows:

[0031] 1) Time series part: First, based on the change points and the characteristics of each segment, combined with the 2021 field survey data, random forest classification is used to obtain the spatial distribution of eucalyptus in 2021; on this basis, the rotation cycle and forest age of eucalyptus in 2021 are determined.

[0032] a. Based on the characteristics of the subsequence NDVI and terrain, the land cover types of each period are classified, and its characteristics include harmonic terms, segmented average NDVI, change trend, segmented length, change amplitude, elevation, slope, and aspect; combined with the 2021 field survey data of the land cover types collected, a random forest classification model for the latest time period is established, and the established model is used to classify all segments to obtain the spatio-temporal distribution results of eucalyptus in 2021.

[0033] b. Determine the eucalyptus planting history, rotation cycle, and forest age: The cutting point and stable point are two main types of changes in the eucalyptus NDVI time series; the stable point is the transition from the rapid growth stage to the undisturbed stable stage; the rules are as follows: (1) If the land cover changes from logging to eucalyptus, mark the start time and end time of this section as the cutting time and stable time; the two sections are divided into one rotation cycle; (2) The change point from non-eucalyptus to eucalyptus is marked as the cutting time / stable time, and the second section is marked as the first planting period; (3) If the land cover changes from eucalyptus to eucalyptus and the age doubles, mark the change point as the cutting time and stable time; (4) If there is only one stable point in the NDVI time series and non-eucalyptus was planted before, the forest age increases by six months; (5) Post-process the classification results: When a section is classified as a logging site and the length of the section lasts for more than three years, re-mark it as eucalyptus. If the duration is less than three years and it cannot be classified due to lack of feature data, mark the section as logging; After determining the eucalyptus planting history in 2021, combined with the eucalyptus spatio-temporal distribution results in step a in 2021, finally obtain the spatial distribution of eucalyptus forest age with a spatial resolution of 30m×30m in 2021;

[0034] 2) Multi-source high-resolution data part: For eucalyptus sub-compartments, adopt the strategy of first vector segmentation and then multi-scale segmentation to re-divide the sub-compartment boundaries;

[0035] Based on eCognition Developer, conduct checkerboard vector segmentation on eucalyptus sub-compartments, overlay the vector segmentation results on the fused 2m image, and perform multi-scale segmentation to obtain the re-segmented eucalyptus sub-compartments.

[0036] In an embodiment of the present invention, the specific implementation of step 4 is as follows:

[0037] Update the eucalyptus forest age information to the re-segmented eucalyptus sub-compartments in step 3. Its essence is to update the pixel values of the forest age raster data to the sub-compartment vector data according to the corresponding regions;

[0038] Using the feature set as the input data, before defining the pixel value extraction range, the feature set needs to be converted to raster data. Automatically convert the vector sub-compartments to raster data based on python+gdal, set the geographic environment, first read the projection information of the raster forest age, and set the projection coordinate system parameters of the raster sub-compartments to be the same as those of the forest age; secondly, set the output pixel size of the raster sub-compartments to 30m×30m, and output raster sub-compartments with the same resolution and projection coordinates;

[0039] Through the pixel center method, capture the pixels of each forest age data that fall within the range of each pixel of the sub-compartment raster, read the pixel value of each pixel of the forest age data, and establish an attribute table.

[0040] In one embodiment of the present invention, in step 5, based on the result of step 4, the table is associated with the small-class raster data, and the pixel value is assigned to the small-class raster.

[0041] In one embodiment of the present invention, in step 6, based on the result of step 5, the small-class raster has extracted the stand age information, and the small-class raster is converted back to vector data. Thus, the stand age update of the eucalyptus small class is completed.

[0042] Compared with the prior art, the present invention has the following beneficial effects: The method of the present invention updates the stand age information on the basis of the re-segmented small class, presents more detailed features in space, and the patches are distinguished by stand age, which is convenient for forestry workers to quickly locate the eucalyptus stand age suitable for operation. The present invention can effectively improve the efficiency of forest stand age update of forest small classes and reduce the cost of obtaining the stand age and rotation period of small classes. This method can be applied to plantations such as eucalyptus, Chinese fir, and rubber. Brief Description of the Drawings

[0043] Figure 1 It is the flow chart of stand age update in the present invention.

[0044] Figure 2 It is the spatial distribution map of eucalyptus age in 2021. Detailed Embodiment

[0045] The technical solution of the present invention will be specifically described below with reference to the drawings.

[0046] As Figure 1 shown, the present invention provides a method for updating the stand age of forest small classes based on multi-source remote sensing data, including the following steps:

[0047] Step 1: Prepare all available Landsat NDVI data from 2000 to 2021 (or the most recent year required by the project), GF-6 image and Sentinel-2 image in 2021, forest field survey data and forest small class data in 2021.

[0048] Step 2: This step can be divided into two parts. 1) Time series part. Randomly divide the NDVI time series of each pixel into M segments, and determine the maximum cumulative sum (CUSUM) of each segment, as well as the corresponding breakpoint quantity, position and order through sSIC calculation; 2) Multi-source high-resolution data part. Preprocess the GF-6 image and Sentinel-2 image respectively, and then perform data fusion.

[0049] Step 3: This step can be divided into two parts. 1) Time series part. Based on the change points and the characteristics of each segment, use a random forest classifier to classify the land use in different periods to obtain the spatial distribution map of eucalyptus. Determine the time, decade, rotation cycle, and forest age of the forest logging and stable recovery in 2021. 2) Multi-source high-resolution data part. Vectorize the small forest classes in the forest survey data and then perform multi-scale segmentation based on the fused high-resolution images to obtain the latest eucalyptus small forest classes.

[0050] Step 4: Based on the results of Step 3, use the latest forest small classes as input data, define the area, capture the forest age raster, and extract the pixel values of the forest age raster.

[0051] Step 5: Based on the results of Step 4, assign the extracted forest age pixel values to the small class raster.

[0052] Step 6: Convert the raster data after assignment into vector small classes. That is, complete the update of the small class forest age.

[0053] The said Step 1: Select Landsat NDVI data (including TM, ETM+ SLC-off, ETM+ SLC-on, and OLI) with cloud cover below 80% in the study area from 2000 to 2021 on USGS. Remove noises such as clouds, shadows, and stripes in each image according to the quality file information.

[0054] Collect sufficient land cover type samples through on-site investigation records and visual interpretation on Google Earth. The classification system can be set according to project requirements. In this example, the focus is on the forest age of eucalyptus plantations, and the land cover is divided into 6 types: eucalyptus, non-eucalyptus forest, farmland, orchard, vegetation-free, and logging sites.

[0055] The main information of the forest small classes includes the average breast diameter, average tree height, age, survey date, number of trees per hectare, etc.

[0056] The multi-source high-resolution data includes GF-6 image data and Sentinel-2 images. The data acquisition time is February 2021.

[0057] The said Step 2: This step includes two parts of work.

[0058] 1) Time series part. First, divide the NDVI time series into M segments; secondly, determine the corresponding breakpoint quantity, location, and order.

[0059] a. Randomly extract some subsequences (M), that is, vectors (NDVI s , NDVI s+1 ,..., NDVI e) where s and e are the start and end observations of a subsequence, and the CUSUM statistic is calculated for each subsequence. Maximize each CUSUM, select the largest CUSUM among the M CUSUMs, and use it as the first candidate break point to be tested according to an information criterion or a predefined threshold. Once a break point is determined to be significant, then the same procedure is applied to the time series on its left and right sides.

[0060]

[0061]

[0062] S < b < e, n = e - s + 1, the s and e points are the start and end points of the time series, b is the observation time in the middle of the s and e points, and NDVI(m, t) is the time series value of a pixel at time t in the m-th sub-segment; is the CUSUM statistic of b in the time series from s to e; f(m, b) is the one with the largest sub-segment and its corresponding change point.

[0063] b. Determine the position and number of break points using sSIC.

[0064]

[0065] α = 1 represents the penalty term of the standard Schwarz criterion (SIC). Here, α = 1.01 is selected to ensure that the results are close to those obtained by SIC; K is the number of candidate break points; is the maximum likelihood estimate of the final variance of the corresponding residuals; N is the number of observations in the time series data.

[0066] 2) Multi-source high-resolution data part. First, perform image preprocessing on GF-6 and Sentinel-2 respectively; then use the GS fusion algorithm to complete image fusion.

[0067] Perform radiometric correction and atmospheric correction on GF-6; perform radiometric calibration and atmospheric correction on Sentinel-2, and uniformly resample the multi-spectral bands to a resolution of 10m.

[0068] Based on the GF-6 panchromatic image, perform geometric registration on the Sentinel-2 multi-spectral image. Use the GS fusion algorithm to perform GS transformation on the resampled multi-spectral image and the simulated panchromatic image, replace the first band of the image after the GS transformation in the previous step with the real panchromatic band, and finally perform GS inverse transformation to obtain a multi-spectral image with a resolution of 2m.

[0069] Step 3 described above: This step includes two parts of work.

[0070] 1) Time series part. First, based on the change points and the characteristics of each segment, combined with the field survey data in 2021, using random forest classification, the spatial distribution of eucalyptus in 2021 was obtained; on this basis, the rotation cycle and forest age of eucalyptus in 2021 were determined.

[0071] a. Based on the characteristics of the subsequence NDVI and terrain, the types of land cover for each period were classified (such as eucalyptus, non-eucalyptus forest, farmland, orchard, vegetation-free, and logging sites). Its characteristics include harmonic terms, segmented average NDVI, change trend, segmented length, change amplitude, elevation, slope, and aspect. Combining the field survey data of land cover types collected in 2021, a random forest classification model for the latest time period was established. The parameters of mtry and ntree were set to their default values. mtry is the number of variables used for the binary tree in the specified node, and the default value is the square root of the number of dataset variables. ntree is the number of decision trees included in the specified random forest, and the default is 500. Using the established model to classify all segments, the spatio-temporal distribution results of eucalyptus in 2021 were obtained.

[0072] b. Determine the planting history, rotation cycle, and forest age of eucalyptus. Cutting points and stable points are the two main types of changes in the NDVI time series of eucalyptus. The stable point is the transition from the rapid growth stage to the stable stage without interference. The rules are as follows: (1) If the land cover changes from logging to eucalyptus, mark the start time and end time of this segment as the cutting time and stable time; two segments are divided into one rotation cycle. (2) The change point from non-eucalyptus to eucalyptus is marked as the cutting time / stable time. The second segment is marked as the first planting period. (3) If the land cover range changes from eucalyptus to eucalyptus and the age doubles, mark the change point as the cutting time and stable time. (4) If there is only one stable point in the NDVI time series and non-eucalyptus was planted before, the forest age increases by 6 months. (5) Post-processing was performed on the classification results: when these segments are classified as logging sites and the length of these segments lasts for more than three years, they are re-marked as eucalyptus. If the duration is less than three years and they are not classified due to lack of feature data, this segment is marked as logging. After determining the planting history of eucalyptus in 2021, combined with the spatio-temporal distribution results of eucalyptus in 2021 in step a, the spatial distribution of the forest age of eucalyptus with a spatial resolution of 30m×30m in 2021 was finally obtained, as Figure 2 shown.

[0073] 2) Multi-source high-resolution data part. For the eucalyptus sub-compartments, a strategy of first vector segmentation and then multi-scale segmentation was adopted to re-divide the sub-compartment boundaries.

[0074] Based on eCognition Developer, perform a checkerboard vector segmentation on the eucalyptus small compartments. Superimpose the vector segmentation results on the fused 2m image and perform multi-scale segmentation to obtain the re-segmented eucalyptus small compartments.

[0075] Step 4 mentioned above: The core work is to update the eucalyptus forest age information into the eucalyptus small compartments re-segmented in Step 3. Its essence is to update the pixel values of the forest age raster data into the small compartment vector data according to the corresponding regions.

[0076] Taking the feature set as the input data, before defining the range of pixel values to be extracted, the feature set needs to be converted into raster data. Automatically convert the vector small compartments into raster data based on python+gdal. Set the geographical environment. First, read the projection information of the raster forest age, and set the projection coordinate system parameters of the raster small compartments to be the same as those of the forest age; secondly, set the output pixel size of the raster small compartments to be 30m×30m. Output raster small compartments with consistent output resolution and projection coordinates.

[0077] Through the pixel center method, capture the pixels of each forest age data that fall within the range of each pixel of the small compartment raster. Read the pixel value of each pixel of the forest age data and establish an attribute table.

[0078] Step 5 mentioned above: Based on the results of Step 4, associate the table with the small compartment raster data and assign the pixel values to the small compartment raster.

[0079] Step 6 mentioned above: Based on the results of Step 5, the forest age information has been extracted from the small compartment raster. Convert the small compartment raster back into vector data. Thus, the forest age update of the eucalyptus small compartments is completed.

[0080] The above are the preferred embodiments of the present invention. All changes made according to the technical solution of the present invention, when the functions and effects produced do not exceed the scope of the technical solution of the present invention, fall within the protection scope of the present invention.

Claims

1. A method for updating the forest stand age of forest sub-compartments based on multi-source remote sensing data, characterized in that, It includes the following steps: Step 1: Prepare all available Landsat NDVI data from 2000 to 2021, GF-6 images and Sentinel-2 images in 2021, forest field survey data and forest sub-compartment data in 2021; Step 2: Divide it into two parts: 1) Time series part: Randomly divide the NDVI time series of each pixel into M segments, and determine the maximum cumulative sum CUSUM, as well as the corresponding breakpoint quantity, position and order of each segment through sSIC calculation; 2) Multi-source high-resolution data part: Preprocess the GF-6 images and Sentinel-2 images respectively, and then perform data fusion; Step 3: Divide it into two parts: 1) Time series part: Based on the change points and the characteristics of each segment, use the random forest classifier to classify the land use in different periods, obtain the spatial distribution map of eucalyptus, and determine the time, decade, rotation cycle and forest age of forest logging and stable recovery in 2021; 2) Multi-source high-resolution data part: Vectorize the sub-compartments in the forest survey data and perform multi-scale segmentation based on the fused high-resolution images to obtain the latest eucalyptus sub-compartments; Step 4: Based on the results of Step 3, use the latest eucalyptus sub-compartments as input data, define the area, capture the forest age raster, and extract the pixel values of the forest age raster; Step 5: Based on the results of Step 4, assign the extracted forest age pixel values to the sub-compartment raster; Step 6: Convert the raster data after assignment into vector sub-compartments, that is, complete the update of the sub-compartment forest age; The specific implementation of Step 2 is as follows: 1) Time series part: First, divide the NDVI time series into M segments; Secondly, determine the corresponding breakpoint quantity, position and order; a. Randomly extract subsequences, namely vectors (NDVI s , NDVI s+1 ,..., NDVI e ), where s and e are the start and end observations of a subsequence, and calculate the CUSUM statistic for each subsequence; maximize each CUSUM, select the largest CUSUM among the M CUSUMs, and use it as the first candidate breakpoint to be tested according to an information criterion or a predefined threshold; once a breakpoint is determined to be significant, perform the same procedure on its left and right time series; s < b < e, n = e - s + 1, where b is the observation time between s and e, and e m 、s m 、n m are the end position, start position, and number of observation values of the m-th segment, and NDVI(m, t) is the time series value of a pixel at time t in the m-th sub-segment; is the CUSUM statistic of b in the time series from s to e; f(m, b) is the sub-segment with the maximum and its corresponding change point; is to take the m value and b value corresponding to the maximum CUSUM; b. Use sSIC to determine the position and quantity of breakpoints: α = 1 represents the standard Schwarz criterion penalty term; K is the number of candidate breakpoints; is the final variance of the maximum likelihood estimate of the corresponding residuals; N is the number of observations in the time series data; 2) Multi-source high-resolution data part: First, perform image preprocessing on GF-6 and Sentinel-2 respectively; then use the GS fusion algorithm to complete image fusion; Perform radiometric correction and atmospheric correction on GF-6; perform radiometric calibration and atmospheric correction on Sentinel-2, and uniformly resample the multi-spectral bands to a resolution of 10m; Based on the GF-6 panchromatic image as a reference, perform geometric registration on the Sentinel-2 multi-spectral image; use the GS fusion algorithm to perform GS transformation on the resampled multi-spectral image and the simulated panchromatic image, replace the first band of the image after GS transformation with the real panchromatic band, and finally perform GS inverse transformation to obtain a multi-spectral image with a resolution of 2m; The specific implementation of Step 3 is as follows: 1) Time series part: First, based on the change points and the characteristics of each segment, combined with the field survey data in 2021, use random forest classification to obtain the spatial distribution of eucalyptus in 2021; On this basis, determine the rotation cycle and forest age of eucalyptus in 2021; a. Classify the types of land cover in each period based on the characteristics of the subsequence NDVI and terrain Its characteristics include harmonic terms, segmented average NDVI, change trends, segmented lengths, change amplitudes, elevations, slopes, and aspects; combining the field survey data of land cover types collected in 2021, a random forest classification model for the latest time period is established, and all segments are classified using the established model to obtain the spatio-temporal distribution results of eucalyptus in 2021; b. Determine the eucalyptus planting history, rotation cycle, and forest age: Cutting points and stable points are the two main types of changes in the eucalyptus NDVI time series; stable points are transitions from the rapid growth stage to the undisturbed stable stage; the rules are as follows: (1) If the land cover changes from logging to eucalyptus, mark the start time and end time of this segment as the cutting time and stable time; two segments are divided into one rotation cycle; (2) The change point from non-eucalyptus to eucalyptus is marked as the cutting time / stable time, and the second segment is marked as the first planting period; (3) If the land cover range changes from eucalyptus to eucalyptus and the age doubles, mark the change point as the cutting time and stable time; (4) If there is only one stable point in the NDVI time series and non-eucalyptus was previously planted, the forest age increases by 6 months; (5) Post-process the classification results: When a segment is classified as a logging site and the length of the segment lasts for more than three years, it is re-marked as eucalyptus. If the duration is less than three years and it is not classified due to lack of feature data, the segment is marked as logging; After determining the eucalyptus planting history in 2021, combining with the spatio-temporal distribution results of eucalyptus in 2021 in step a, the spatial distribution of eucalyptus forest age with a spatial resolution of 30m×30m in 2021 is finally obtained; 2) Multi-source high-resolution data part: For eucalyptus sub-compartments, adopt the strategy of first vector segmentation and then multi-scale segmentation to re-divide the sub-compartment boundaries; Based on eCognition Developer, perform checkerboard vector segmentation on eucalyptus sub-compartments, overlay the vector segmentation results on the fused 2m image, and perform multi-scale segmentation to obtain re-segmented eucalyptus sub-compartments.

2. The method for updating the forest stand age of forest sub-compartments based on multi-source remote sensing data according to claim 1, wherein In step 1, select Landsat NDVI data with cloud cover below 80% in the study area from 2000 to 2021 on USGS, including TM, ETM+SLC-off, ETM+SLC-on, and OLI; remove the noise including clouds, shadows, and stripes in each image according to the quality file information; Collect land cover type samples through field survey records and visual interpretation on Google Earth; set the classification system according to project requirements, and divide the land cover into 6 types: eucalyptus, non-eucalyptus forest, farmland, orchard, vegetation-free, and logging site; The main information of forest sub-compartments includes the average breast diameter, average tree height, age, survey date, and number of trees per hectare of the sub-compartment; The multi-source high-resolution data includes GF-6 image data and Sentinel-2 images.

3. A method for updating the forest stand age of forest sub-compartments based on multi-source remote sensing data according to claim 1, characterized in that, The specific implementation of step 4 is as follows: Update the eucalyptus forest age information to the re-segmented eucalyptus sub-compartments in step 3, which essentially updates the pixel values of the forest age raster data to the sub-compartment vector data according to the corresponding regions; Taking the feature set as the input data, before defining the extraction range of pixel values, the feature set needs to be converted into raster data. Based on Python + GDAL, the vector sub-compartments are automatically converted into raster data. Set the geographic environment. First, read the projection information of the raster forest age, and set the projection coordinate system parameters of the raster sub-compartments to be the same as those of the forest age. Secondly, set the output pixel size of the raster sub-compartments to 30m × 30m, and output raster sub-compartments with the same output resolution and projection coordinates. Through the pixel center method, capture the pixels of each forest age data falling within the range of each pixel of the sub-compartment raster, read the pixel value of each pixel of the forest age data, and establish an attribute table.

4. A method for updating the forest stand age of forest sub-compartments based on multi-source remote sensing data according to claim 1, characterized in that In step 5 described above, based on the result of step 4, associate the table with the sub-compartment raster data, and assign the pixel value to the sub-compartment raster.

5. A method for updating the forest stand age of forest sub-compartments based on multi-source remote sensing data according to claim 1, characterized in that In step 6 described above, based on the result of step 5, the sub-compartment raster has extracted the forest age information, and convert the sub-compartment raster back to vector data. Thus, the forest age update of the eucalyptus sub-compartments is completed.

Images