CN118521916A

CN118521916A - Forest disturbance detection method and system

Info

Publication number: CN118521916A
Application number: CN202410969278.3A
Authority: CN
Inventors: 洪亮; 何丽; 刘仁
Original assignee: Yunnan Normal University
Current assignee: Yunnan Normal University
Priority date: 2024-07-19
Filing date: 2024-07-19
Publication date: 2024-08-20
Anticipated expiration: 2044-07-19
Also published as: CN118521916B

Abstract

The invention provides a forest disturbance detection method and system, and belongs to the field of forest protection and monitoring. Firstly, obtaining Landsat series data and related data of preset months in a detection year in a detection area; identifying a forest area according to the related data; in a forest area, calculating NBR indexes of pixels corresponding to each Landsat image based on Landsat series data, and obtaining a complete first NBR time sequence track through a maximum value method and an interpolation method; generating a second NBR time sequence track by pixel fitting through an SG filter with constraint; then LANDTRENDR algorithm is utilized to process the second NBR time sequence track to detect the pixels and years of forest disturbance; and analyzing the difference between the two NBR time sequence tracks before and after the forest disturbance year, and correcting the year. The reconstruction of high-density time sequence data is carried out by utilizing the data filling method, and the NBR time sequence track of the pixels is fitted by utilizing the SG filtering algorithm with constraint, so that the detection accuracy of forest disturbance events is improved.

Description

Forest disturbance detection method and system

Technical Field

The invention belongs to the field of forest protection and monitoring, and particularly relates to a forest disturbance detection method and system.

Background

Forest disturbance refers to the phenomenon of forest canopy coverage and obvious reduction or disappearance of forest biomass caused by nature (drought, storm, fire, plant diseases and insect pests, etc.) or artificial driving factors (such as forest deforestation, townization, forest disfigurement, etc.). Forest disturbance is the main driving force of forest ecological dynamic change, disturbance history influences the growth state of a forest stand, and the type and the intensity of disturbance change the species composition and the structure of the forest stand. Accurate, timely and continuous acquisition of forest disturbance information is important to forest sustainable management and regional carbon sink estimation.

Currently, landsat time series data is generally used for forest disturbance detection. The Landsat time series data are medium resolution (30 m) Landsat series satellite data which are continuously observed in long time sequence (1972 to date) and downloaded in free opening, and can provide continuous and important data sources for forest disturbance detection. Forest disturbance detection based on Landsat time series data generally adopts a remote sensing change detection method, for example: track change detection algorithms (Trajectory-based Change Detection Algorithm, TBCD), additional season and region interruption algorithms (Breaks For Additive Season and Trend, BFAT), forest disturbance continuous detection algorithms (The Continuous Monitoring of Forest Disturbance Algorithm, CMFDA), continuous change detection and classification (Continuous Change Detection and Classification, CCDC), landsat-based disturbance and recovery trend detection (Landsat-based Detection of Trends in Disturbance and Recovery, LANDTRENDR) algorithms and the like, which all utilize statistical program set thresholds to extract disturbance signals from spectrum time sequence tracks, and have differences in time frequency of remote sensing images, indexes of detection change, real-time detection and the like, so that forest disturbance is expanded from forest stand substitution to detection of gradual-type forest disturbance events caused by insect pests, meta-cuts, forest degradation and the like.

In the prior art, LANDTRENDR algorithm detects gradual change and abrupt change forest interference events by utilizing spectrum time sequence tracks formed by cloud-free images in the same season each year, the integration of multiple-stage Landsat images to select cloud-free pixels is flexible, and the method is realized on a GEE platform, provides a new opportunity for rapidly and efficiently detecting large-scale forest disturbance, and becomes one of the most widely applied algorithms for forest disturbance detection at present. However, the algorithm captures the forest disturbance occurrence time in the spectrum time sequence track by using a time segmentation method, on the premise of not losing important detail characteristics, the p value of the time sequence data F test is utilized to iteratively simplify the original spectrum time sequence track of each pixel to a small number of straight line segments for recording when related signals change, and finally forest disturbance detection is carried out by using the simplified time sequence track and vegetation coverage change criterion, so that the influence of the quality of long time sequence Landsat data and the complexity of disturbance events on large-area-scale forest disturbance detection is ignored; in addition, LANDTRENDR algorithm reduces the influence of the over-fitting and under-fitting phenomena of the time division algorithm on the remote sensing detection of the disturbance event of the complex forest by setting a series of control parameters and filtering processes. However, in a complex terrain area, forest stand species are complex in composition, various in forest disturbance types and large in disturbance intensity space heterogeneity, high-precision detection of mutant type and gradual type forest disturbance events is difficult to synchronously realize by a single control parameter, and delay or advance phenomenon occurs to the detection time of some forest disturbance events.

Disclosure of Invention

In view of the above-mentioned drawbacks or shortcomings in the prior art, the present invention aims to provide a forest disturbance detection method and system, based on an improved LANDTRENDR algorithm, the LANDTRENDR algorithm is improved from two aspects of quality of Landsat time series data and difference of spectrum time series track, reconstruction of high-density time series data is performed by using a data filling method, real change information of the data is reserved, and the spectrum time series track of pixels is fitted by using an SG filtering algorithm with constraint, so that the spectrum information of important time nodes of a mutant forest disturbance event is reserved to the greatest extent, and detection precision of the complex forest disturbance event is improved.

In order to achieve the above purpose, the embodiment of the present invention adopts the following technical scheme:

In a first aspect, an embodiment of the present invention provides a method for detecting forest disturbance, where the method includes the following steps:

step S1, landsat series data of preset months in a detection year in a detection area are obtained, and Sentinel-1/2 data of the last year in the detection year and 30m resolution digital elevation model SRTM DEM data acquired by a space plane radar topography mapping task are obtained;

S2, dividing a detection area into a forest area and a non-forest area according to Sentinel-1/2 data and SRTM DEM data of the last year in the detection year;

Step S3, preprocessing Landsat series data based on the divided forest areas, taking a normalized burning rate NBR as a forest disturbance detection index, calculating an NBR index of a pixel corresponding to each Landsat image in the forest areas, synthesizing the NBR maximum value of the pixel corresponding to the current pixel in each year to obtain cloud-removed NBR time sequence data of the detection year, filling the NBR data of the pixel in the missing year by using an interpolation method, reconstructing to obtain the NBR time sequence data of the complete detection year of the pixel in the forest areas, and generating an NBR time sequence track observed by a satellite as a first NBR time sequence track;

S4, eliminating noise of the NBR time sequence data by utilizing a constraint SG filter and reserving spectrum information of important time nodes of a forest disturbance event, and fitting pixel by pixel to generate a continuous NBR time sequence track with the constraint SG filter as a second NBR time sequence track;

S5, detecting pixels and years of forest disturbance by utilizing LANDTRENDR algorithm based on NBR time sequence track with constraint SG filtering in the range of the divided forest area, and assigning the years to the forest disturbance pixels;

S6, analyzing differences between the first NBR time sequence track and the second NBR time sequence track before and after the forest disturbance year, and correcting the year of the forest disturbance pixel; the year of the forest disturbance pixel is corrected, which comprises the following steps:

step S61, comparing the first NBR time sequence track and the second NBR time sequence track before and after the forest disturbance year, and dividing the comparison result into eight types:

Comparing two NBR time sequence tracks from the previous year of the forest disturbance year to the forest disturbance year, wherein the comparison result is divided into 4 types, and ① is the trend that the first NBR time sequence track and the second NBR time sequence track are both declining; ② is that the first NBR timing trace is up and the second NBR timing trace is down; ③ th is that the first NBR time sequence track is lowered, and the second NBR time sequence track is raised; ④ is that the first NBR time sequence track and the second NBR time sequence track are in an ascending trend;

Comparing two NBR timing traces from the year of forest disturbance to the year after the disturbance year, there are also four types: ⑤ is that the first NBR time sequence track and the second NBR time sequence track are in an ascending trend; ⑥ is that the first NBR timing trace is up and the second NBR timing trace is down; ⑦ th is that the first NBR timing trace is down and the second NBR timing trace is up; ⑧ is that the first NBR time sequence track and the second NBR time sequence track are in a descending trend;

Step S62, correcting the years of the detected forest disturbance pixels according to eight different types of comparison results, wherein the correction process is as follows:

When the scene ①、⑤ is met, the fact that two time sequence tracks are consistent is indicated, and the forest disturbance time detected by the LANDTRENDR algorithm accords with the actual situation; when the scene ②、④、⑥ is satisfied, the year in which the disturbance actually occurs is corrected to be LANDTRENDR, which is the year before the year in which the algorithm detects the year; when the scene ③、⑦、⑧ is satisfied, the year correction of the disturbance is LANDTRENDR, and the calculation formula is shown in formula (5):

（5）

In the formula (5), the amino acid sequence of the compound, Is the corrected forest disturbance year,Is the year of forest disturbance detected by LANDTRENDR algorithm,Is the comparison result type of two NBR time sequence tracks;

and S7, evaluating the precision of the forest disturbance detection result.

As a preferred embodiment of the present invention, the dividing process of the forest area and the non-forest area in step S2 is as follows:

step S21, dividing a detection area into four areas by using a grid segmentation method, randomly generating M sample points in each area, taking the sample points which belong to forest types in a corresponding database of the sample points and have a tree height of more than 5M as forest samples, and taking the rest as non-forest samples;

and S22, carrying out forest/non-forest classification mapping by using a random forest algorithm based on the polarization and spectrum characteristics of the median synthesized image of the Sentinel-1 and the cloud-free Sentinel-2 and the altitude, gradient and slope characteristics calculated by the SRTM DEM data.

As a preferred embodiment of the invention, the Landsat series of data comprises Landsat Collection Tier1 surface reflectance images available on the Google Earth Engine GEE platform for a predetermined month each year of the year.

As a preferred embodiment of the present invention, the Landsat Collection t 2t ir 1 surface reflectance image includes: the image collected by three sensors of the thematic drawing instrument TM, the enhanced thematic drawing instrument ETM+, and the land imager OLI.

As a preferred embodiment of the present invention, preprocessing the Landsat series data in step S3 includes:

step S31, performing cloud/shadow/snow masking by utilizing CFmask algorithm;

step S32, performing normalization processing on all Landsat wave bands collected by the OLI by using a correction coefficient to reduce the influence of differences among TM, ETM+ and OLI sensors on forest disturbance detection;

Step S33, calculating NBR indexes of pixels corresponding to each Landsat image, synthesizing NBR time sequence data in detection years in a forest area by using the maximum value of NBR in a preset month each year, and generating an NBR time sequence track observed by a satellite based on the NBR time sequence data as a first NBR time sequence track.

As a preferred embodiment of the present invention, the NBR has the following formula:

（1）

In the formula (1), the components are as follows, AndRespectively the earth surface reflectivities of near infrared and short wave infrared wave bands in Landsat images;

Year of detection NBR timing trace for satellite observationsThe method comprises the following steps:

（2）

In the formula (2), the amino acid sequence of the compound, Is the 1 st year of the detection yearThe maximum value of NBR of the picture element,Is the 2 nd year in the detection yearThe maximum value of NBR of the picture element,Is the last year of the detection yearsThe maximum value of NBR of the picture element.

As a preferred embodiment of the present invention, in step S3, the interpolation method is used to fill NBR data of pixels in missing years, and NBR time sequence data of complete detection years of pixels in a forest area is obtained by reconstruction, specifically comprising:

Based on two data adjacent to the missing year, performing numerical estimation by using a linear relation, wherein a calculation formula is shown in a formula (3):

（3）

in the formula (3), the amino acid sequence of the compound, Is of the missing yearThe value of the speculation is calculated,Representing the year of the absence of the drug,AndRespectively, the upper limit of the interpolation intervalAnd lower limit ofYear corresponds toAnd (5) observing values.

As a preferred embodiment of the present invention, step S4 of NBR timing trace fitting with constrained SG filtering comprises:

step S41, fitting the NBR time sequence data observed by the satellite by using an SG filtering algorithm to obtain an SG filtered NBR time sequence track of the detection year;

And step S42, correcting the SG-filtered NBR time sequence track by using the constraint condition to obtain the NBR time sequence track with the constraint SG filter.

As a preferred embodiment of the invention, the SG filtering algorithm with constraint is adopted to correct the primary NBR time sequence track, which specifically comprises the following steps:

when the absolute change rate of the pixel NBR of the satellite-observed NBR time sequence track and the SG-filtered NBR time sequence track is larger than 0.2, the satellite-observed NBR is used for replacing the NBR at the time point, otherwise, the SG-filtered NBR is reserved, and the calculation formula is shown in the formula (4):

（4）

In the formula (4), the amino acid sequence of the compound, AndNBR observed by satellites and NBR estimated by SG filters are respectively integrated with all pixels to obtain NBR time sequence tracks with constraint SG filters.

In a second aspect, an embodiment of the present invention further provides a forest disturbance detection system, where the system includes: the system comprises a data acquisition module, a forest area identification module, a first NBR time sequence track generation module, a second NBR time sequence track generation module, a LANDTRENDR algorithm detection module, a year correction module and a precision evaluation module; wherein,

The data acquisition module is used for acquiring Landsat series data of preset months in the detection year in the detection area, sentinel-1/2 data of the last year in the detection year and 30m resolution digital elevation model SRTM DEM data acquired by a space plane radar topography mapping task;

The forest area identification module is used for dividing the detection area into a forest area and a non-forest area according to the Sentinel-1/2 data and the SRTM DEM data of the last year in the detection year;

The first NBR time sequence track generation module is used for preprocessing Landsat series data based on the divided forest area, taking a normalized combustion rate NBR as a forest disturbance detection index, calculating an NBR index of each Landsat image corresponding to a pixel in the forest area, synthesizing the cloud removed NBR time sequence data of the detection year by utilizing the NBR maximum value of the corresponding current pixel in each year, filling the NBR data of the pixel in the missing year by utilizing an interpolation method, reconstructing to obtain the NBR time sequence data of the complete detection year of the pixel in the forest area, and generating an NBR time sequence track observed by a satellite as a first NBR time sequence track;

The second NBR time sequence track generation module is used for eliminating noise of the NBR time sequence data by utilizing the SG filter with constraint and keeping spectrum information of important time nodes of forest disturbance events, and generating a continuous NBR time sequence track with constraint SG filter by pixel fitting as a second NBR time sequence track;

The LANDTRENDR algorithm detection module is used for detecting pixels and years of forest disturbance occurrence by utilizing a LANDTRENDR algorithm based on NBR time sequence tracks with constraint SG filtering in the divided forest area range, and assigning the years to the forest disturbance pixels;

The year correction module is used for analyzing differences between the first NBR time sequence track and the second NBR time sequence track before and after the forest disturbance year and correcting the year of the forest disturbance pixel; when correcting the year of the forest disturbance pixel, the following steps are further executed:

（5）

The precision evaluation module is used for evaluating the precision of the forest disturbance detection result.

The technical scheme provided by the embodiment of the invention has the following beneficial effects:

The forest disturbance detection method and system provided by the embodiment of the invention effectively reserve the real spectrum information of important time nodes of various forest disturbance events, and compared with the original LANDTRENDR algorithm, the forest disturbance detection method and system are more accurate in detecting the forest disturbance events. By adopting the forest disturbance detection method to detect forest disturbance, the main concentrated year of occurrence of forest disturbance event can be accurately obtained, and by adopting and actively implementing related forest protection policies, the forest disturbance area is obviously reduced, and the forest is effectively protected.

Of course, it is not necessary for any one product or method of practicing the invention to achieve all of the advantages set forth above at the same time.

Drawings

In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for the description of the embodiments will be briefly described below, and it is apparent that the drawings in the following description are only some embodiments of the present invention, and other drawings may be obtained according to these drawings without inventive effort for a person of ordinary skill in the art.

FIG. 1 is a flow chart of a forest disturbance detection method according to an embodiment of the present invention;

FIG. 2 is a diagram of the original SG filtering result in the prior art;

FIG. 3 is a graph showing the results of the improved SG filtering in an embodiment of the invention;

FIG. 4 is a graph showing the difference between the satellite observations and the SG filtered spectrum timing traces in an embodiment of the present invention;

Fig. 5 is a visual interpretation of forest disturbance in an embodiment of the invention: (a) NBR timing traces, (b 1-b 8) are Landsat images of important time nodes (false color synthesis), respectively;

FIG. 6 is a graph of comparing forest disturbance mapping accuracy with statistics for different experimental schemes in an embodiment of the invention;

fig. 7 is an area ratio diagram of forest disturbance in yunnan province in an embodiment of the present invention: (a-d) are statistical graphs of forest disturbance area occupation ratios of different years, different elevation areas, annual average temperature and annual average precipitation areas respectively.

Detailed Description

The present inventors have conducted intensive studies on the existing forest disturbance detection method after finding the above-mentioned problems. The quality of Landsat time sequence data is found to be one of key elements for accurately detecting forest disturbance by LANDTRENDR algorithm. The cloud-free Landsat image with better quality in the season growing every year can be used for achieving a forest disturbance detection result with higher precision. However, high quality Landsat time series in cloudy rain areas (e.g., tropical and subtropical) are substantially sparse and irregularly distributed, and the lack of images during critical seasons causes spectral timing traces to be over-fitted, thereby affecting the accuracy of change detection. At present, filling and fitting of time sequence data are main methods for solving the problems, namely, firstly filling spectrum values of missing years by using a moving average method, and then removing data noise by using SG filtering (Savitzky-Golay filter) to obtain long-time sequence Landsat data. The SG filtering utilizes a least square method to calculate a smoothing coefficient corresponding to a higher-order polynomial function between continuous spectrum values in a sliding window, and utilizes the smoothing coefficient and discrete convolution to calculate a smoothing value of a window center point, and can select different window widths at any position of a spectrum change curve so as to meet the requirements of fitting different forest disturbance type spectrum time sequence tracks. However, since the spectrum timing trace of the abrupt forest disturbance event (such as forest fire, forest felling, etc.) does not conform to the assumption of polynomial fitting (i.e., the assumption data is generated by a polynomial function with respect to some unknown parameters), the SG filter may misjudge the abrupt signal of the disturbance time node as data noise, thereby causing pseudo-variation information in the fitted spectrum timing trace. Therefore, the uncertainty of LANDTRENDR algorithm on mutant forest disturbance event detection can be reduced by improving the SG filtering algorithm to keep the spectrum information of important time nodes.

In addition, in general, the spectrum value of the year in which the forest disturbance occurs is necessarily smaller than that of the previous year or the next year, and the satellite observation and the spectrum time sequence track after SG filtering between the previous year and the next year of the forest disturbance should be consistent. Therefore, based on the forest disturbance result detected by LANDTRENDR algorithm, the remote sensing detection precision of the complex forest disturbance event is improved by correcting the detection result by utilizing the satellite observation before and after the disturbance year and the difference of the spectrum track after SG filtering.

It should be noted that the above drawbacks and solutions of the prior art solutions are all results obtained by the inventor after practice and careful study, and thus the discovery process of the above problems and the solutions presented below by the embodiments of the present invention for the above problems should be all contributions of the inventor to the present invention during the process of the present invention.

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. It will be apparent that the described embodiments are only some, but not all, embodiments of the invention. The components of the embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations. It should be noted that, in the case of no conflict, the embodiments of the present invention and features in the embodiments may also be combined with each other.

It should be noted that: like reference numerals and letters denote like items in the following figures, and thus once an item is defined in one figure, no further definition or explanation thereof is necessary in the following figures. In the description of the present invention, the terms "first," "second," "third," "fourth," and the like are used merely to distinguish between descriptions and are not to be construed as indicating or implying relative importance.

Based on the above analysis, the embodiment of the invention provides a forest disturbance detection method and system, because LANDTRENDR algorithm essentially uses the intrinsic drive of data to divide a time sequence track into a series of straight line segments, uses the time and spectrum value of the dividing point to perform forest disturbance detection, ignores the influence of the quality of long time sequence Landsat data and the complexity of disturbance event on large-scale forest disturbance detection, improves LANDTRENDR algorithm from the two aspects of the quality of Landsat time sequence data and the difference of NBR time sequence track, and improves the detection precision of complex forest disturbance event, and specifically comprises: firstly, reconstructing long-time-sequence Landsat data by using a linear interpolation method based on seasonal synthetic data with frequent regional forest disturbance events and higher data availability; secondly, eliminating data noise by utilizing an SG filter with constraint, and reserving spectrum information of important time nodes of a mutation type forest disturbance event, and generating a continuous and reliable NBR time sequence track pixel by pixel; finally, carrying out forest disturbance mapping by utilizing LANDTRENDR algorithm, and further correcting the detection result of the complex forest disturbance event based on satellite observation and the difference of NBR time sequence track after SG filtering; and the following objectives are achieved: (1) Based on Landsat data in seasons with frequent regional forest disturbance events and higher data availability, reconstructing high-density time sequence data by using a data filling method and reserving real change information of the data; (2) Fitting a spectrum time sequence track of the pixel by using a constrained SG filtering algorithm, so that spectrum information of important time nodes of the mutant forest disturbance event is reserved to the greatest extent; (3) Carrying out forest disturbance detection by utilizing LANDTRENDR algorithm, and further correcting the detection result based on satellite observation and the difference of the NBR time sequence track after SG filtering; (4) The effectiveness of the improved LANDTRENDR algorithm on forest disturbance mapping is quantitatively evaluated.

As shown in fig. 1, the forest disturbance detection method provided by the embodiment of the invention includes the following steps:

Step S1, landsat series data of preset months in the detection year in the detection area, sentinel-1/2 data of the last year in the detection year and 30m resolution digital elevation model (Shuttle RADAR TERRAIN Mission DigitalElevationModel, SRTM DEM) data acquired by a space plane radar topography mapping task are acquired.

In this step, landsat series data includes Landsat Collection Tier1 surface reflectance images available from Google Earth Engine (Google EARTH ENGINE, GEE) platforms for a predetermined month each year in the year, e.g., GEE platform available from 1 month 1 day to 4 months 1 day each year in 1986-2023. The Landsat Collection Tier1 surface reflectance image includes: images collected by three sensors, a thematic mapping instrument (THEMATIC MAPPER, TM), an enhanced thematic mapping instrument (ENHANCED THEMATIC MAPPER Plus, etm+), and a terrestrial imager (Operational LAND IMAGER, OLI).

And S2, dividing the detection area into a forest area and a non-forest area according to the Sentinel-1/2 data and the SRTM DEM data of the last year in the detection year.

In this step, the dividing process of the forest area and the non-forest area is specifically as follows:

S21, dividing a detection area into four areas by using a grid segmentation method, and randomly generating M sample points in each area, wherein M is preferably more than or equal to 5000; and mapping these sample points to a database, for example, a database of ESAWorldCover m 2021 v200 products, wherein sample points belonging to forest types and having a tree height of >5m are used as forest samples, and the rest are used as non-forest samples.

And S3, preprocessing the Landsat series data based on the divided forest areas, taking a normalized burning rate (Normalized Burn Ratio, NBR) as a forest disturbance detection index, calculating an NBR index of each Landsat image corresponding to a pixel in the forest areas, synthesizing to obtain cloud-removed NBR time sequence data of the detection year by using the NBR maximum value of each year corresponding to the current pixel, filling the NBR data of the pixel in the missing year by using an interpolation method, reconstructing to obtain the NBR time sequence data of the complete detection year of the pixel in the forest areas, and generating an NBR time sequence track observed by satellites as a first NBR time sequence track.

In this step, in order to ensure the quality of the NBR time series data, preprocessing the Landsat series data includes:

step S31, performing cloud/shadow/snow masking by utilizing CFmask algorithm;

Step S33, calculating NBR indexes of pixels corresponding to each Landsat image, synthesizing NBR time sequence data in detection years in a forest area by using the maximum value of NBR in a preset month each year, and generating an NBR time sequence track observed by a satellite based on the NBR time sequence data as a first NBR time sequence track. Taking forest disturbance detection in 1986-2023 of Yunnan province as an example, taking the NBR corresponding to healthy and dense forests as a forest disturbance detection index, the embodiment selects the NBR as the forest disturbance detection index in consideration of relatively high NBR and high sensitivity to forest disturbance and recovery; the NBR indexes of pixels corresponding to each Landsat image are calculated respectively, landsat series data of 1 month 1 day to 4 months 1 day each year in 1986-2023 are taken as an example, NBR maximum values of 1 month to 4 months each year are utilized to synthesize NBR time sequence data of pixels in a forest region of a detection region in 1986-2023 years, and an NBR time sequence track observed by a satellite is generated based on the NBR time sequence data to be used as a first NBR time sequence track.

Wherein, the calculation formula of NBR is:

（1）

In the formula (1), the components are as follows, AndThe surface reflectivities of near infrared and short wave infrared bands in the Landsat image, respectively.

Year of detectionNBR timing trace for satellite observationsThe method comprises the following steps:

（2）

For example, NBR timing data from 1986-2023 are:。

The reconstruction is carried out to obtain NBR time sequence data of complete detection years of pixels in a forest area, the NBR time sequence data are comprehensively influenced by complex terrains and climate conditions, high-quality cloud-free Landsat images in the forest area are distributed unevenly in space, landsat data accumulation years in different areas are different, and data loss conditions exist in partial areas in certain years. In order to generate complete and continuous NBR time sequence data in the detection year, the step fills the NBR in the missing year by using the NBR in the existing year and a linear interpolation method, namely, based on two adjacent data before and after the missing year, numerical value estimation is carried out by using a linear relation, and the calculation formula is shown as a formula (3):

（3）

And S4, eliminating noise of the NBR time sequence data by using the SG filter with the constraint and reserving spectrum information of important time nodes of the forest disturbance event, and fitting pixel by pixel to generate a continuous NBR time sequence track with the constraint SG filter as a second NBR time sequence track.

In this step, NBR timing track fitting with constrained SG filtering includes:

step S41, fitting the NBR time sequence data observed by the satellite by utilizing an SG filtering algorithm to obtain an SG filtered NBR time sequence track of the detection year;

In this step, the SG filter is widely used for fitting of spectrum timing tracks because it can reduce the influence of data noise caused by cloud detection errors, sensor degradation and other factors in long-time-series data to a large extent and effectively retain spectrum variation information. The core idea of the filtering algorithm based on the SG filter is to perform weighted filtering on the data in the sliding window, wherein the weight is obtained by performing least square fitting on a high-order polynomial of NBR data in the window. However, since the forest disturbance event belongs to the mutation data, the mutant NBR time sequence track does not satisfy the assumption of polynomial linear fitting, so that the SG filter often fails to fit the NBR time sequence track of the sudden forest disturbance event (such as forest fire, forest felling, etc.), as shown in fig. 2. Based on the above, the step adopts the SG filtering algorithm with constraint to correct the SG filtered NBR time sequence track, and specifically comprises the following steps:

（4）

Fig. 3 shows SG filtering results with constraints. As can be seen from fig. 2 and 3, the NBR timing trace fitted with the constrained SG filtering algorithm is more closely related to the actual spectral trace of the forest disturbance event.

And S5, detecting pixels and years of forest disturbance by utilizing LANDTRENDR algorithm based on NBR time sequence track with constraint SG filtering in the divided forest area, and assigning years to the forest disturbance pixels.

In this step, when the disturbance and recovery trend detection (Landsat-based detection of trends in disturbance and recovery, LANDTRENDR) algorithm based on Landsat is used to detect the disturbance year t of the forest, when the disturbance of the forest occurs, the NBR time sequence track will display a numerical value drop, and when the forest recovers, the NBR index gradually rises, and the LANDTRENDR algorithm will determine the year corresponding to the node where the NBR "drop-rise" change occurs as the year where the disturbance of the forest occurs.

And S6, analyzing differences between the first NBR time sequence track and the second NBR time sequence track before and after the forest disturbance year, and correcting the year of the forest disturbance pixel.

The method for correcting the year of the forest disturbance pixel comprises the following steps:

Step S61, comparing the first NBR time series track and the second NBR time series track before and after the forest disturbance year, and classifying the comparison result into eight types, as shown in fig. 4:

Comparing two NBR timing traces from the year of forest disturbance to the year after the disturbance year, there are also four types: ⑤ is that the first NBR time sequence track and the second NBR time sequence track are in an ascending trend; ⑥ is that the first NBR timing trace is up and the second NBR timing trace is down; ⑦ th is that the first NBR timing trace is down and the second NBR timing trace is up; ⑧ is that both the first NBR timing trace and the second NBR timing trace have a decreasing trend.

（5）

In the formula (5), the amino acid sequence of the compound, Is the corrected forest disturbance year,Is the year of forest disturbance detected by LANDTRENDR algorithm,Is the comparison result type of the two NBR timing tracks.

And S7, evaluating the precision of the forest disturbance detection result.

In the step, the effectiveness of forest disturbance event detection is evaluated by using a LANDTRENDR algorithm which is improved by quantitative evaluation of three verification samples from different sources, namely, the obtained corrected year forest disturbance pixels and years are evaluated. Quantitative evaluation indexes include overall accuracy (Overallaccuracy, OA), omission factor (Omission error, OE), error (CE), and F1-Score, which are calculated by the following formulas, respectively:

（6）

（7）

（8）

（9）

In the formula (6-9), The total number of pixels with the detection year consistent with the forest disturbance year is the number of pixels with the correct detection of the forest disturbance year; the total number of pixels of which the detection year is inconsistent with the forest disturbance year is the total number of pixels of error detection; The total number of samples without forest disturbance is the total number of pixels without disturbance when the forest disturbance is actually generated but the model detection result is that the disturbance is not generated; Respectively verifying the total number of pixels in the sample, which are subjected to forest disturbance and are not subjected to forest disturbance; Is the total number of verification samples, 。

Based on the same thought, the embodiment of the invention also provides a forest disturbance detection system, which comprises a data acquisition module, a forest area identification module, a first NBR time sequence track generation module, a second NBR time sequence track generation module, LANDTRENDR algorithm detection module, a year correction module and an accuracy evaluation module; wherein,

The data acquisition module is used for acquiring Landsat series data of preset months in the detection year in the detection area, and Sentinel-1/2 data and SRTM DEM data of the last year in the detection year;

the forest area identification module is used for dividing the detection area into a forest area and a non-forest area according to the Sentinel-1/2 data of the last year in the detection year and the SRTM DEM data;

the first NBR time sequence track generation module is used for preprocessing Landsat series data based on a divided forest area, taking a normalized burning rate NBR as a forest disturbance detection index, calculating an NBR index of each Landsat image corresponding to a pixel in the forest area, synthesizing the NBR maximum value of each year corresponding to the current pixel to obtain cloud removed NBR time sequence data of a detection year in the detection area, filling the NBR data of the pixel in the missing year by using an interpolation method, reconstructing to obtain complete NBR time sequence data of all pixels, and generating an NBR time sequence track observed by a satellite as a first NBR time sequence track;

The year correction module is used for analyzing differences between the first NBR time sequence track and the second NBR time sequence track before and after the forest disturbance year and correcting the year of the forest disturbance pixel;

Taking Yunnan province of China as an example, an interpretation example of the forest cutting event shown in FIG. 5 is that the time of forest disturbance is primarily determined by using NBR time sequence tracks in 1986-2023; and then selecting at least two clear images before and after the time of the suspected last occurrence of forest disturbance to further confirm the occurrence time of the forest disturbance. Through the method, 618 forest disturbance verification samples are selected in Yunnan province in total in actual operation, and the spatial distribution of the forest disturbance verification samples is shown in fig. 6.

The method for detecting forest disturbance based on the improved LANDTRENDR algorithm provided by the embodiment is adopted to execute the steps S1-S6 by taking Yunnan province in southwest China as a detection area, and the forest disturbance detection in 1986-2023 in the detection area is carried out.

The elevation of the joint of China, southeast Asia and south Asia in Yunnan province is declined from north (6740 m) to south (76.4 m), and the terrain types are complex and various, including terrains such as basin, hills, mountain land, plateau and the like. In addition, the research area belongs to subtropical plateau monsoon climate, but because of the complex terrain condition, the space and vertical difference of local climate are larger, so that the space heterogeneity of the temperature and precipitation condition of the area is obvious, complex and various forest landscapes are formed, and abundant forest resource endowments are included. From the 2021 government work report, the forest of Yunnan province covers about 2392.65 ten thousand hectares, covers about 65.04 percent, and has an accumulation of about 20.67 hundred million cubic meters. However, due to the comprehensive influence of factors such as climate, resources, geographical location and the like, the forest fire risk and the potential pest risk level in Yunnan province are high, and particularly forest disturbance events frequently occur in winter and spring (11 months to 5 months in the next year) with drought, little rain, sufficient illumination and mild and dry climate.

In the 1986-2023 Yunnan forest disturbance detection result manufactured by the method of the embodiment, in the division of the forest area and the non-forest area, 5096 sample points which meet the requirements are counted, wherein 3126 forest samples and 1970 non-forest samples are obtained; 1529 forest/non-forest samples (911 forests, 618 non-forests) are randomly selected by using a high-resolution Google Earth image in 2023 to verify drawing precision, and the verification result shows that the overall precision of forest/non-forest classification in Yunnan province is 88.07% and the kappa coefficient is 0.7407.

In order to quantitatively evaluate the effectiveness of the improved LANDTRENDR algorithm on the remote sensing detection of the disturbance event of the complex forest, 6 experimental schemes are set, as shown in table 1; and quantitatively evaluating the detection precision of the algorithm on the complex forest disturbance event by using the verification samples from three aspects of time sequence data reconstruction, spectrum time sequence track fitting and forest disturbance detection respectively, as shown in table 2. The schemes I and II are used for evaluating influences of Landsat time sequence data in different time periods on forest disturbance detection, the schemes II and III are used for evaluating influences of different data filling methods on results, the schemes IV and V are used for evaluating contributions of a constraint SG filtering algorithm on the results, and the schemes V and VI are used for evaluating contributions of a correction method of the LANDTRENDR algorithm on forest disturbance detection. In addition, all experiments in this example were based on the GoogleEarthEngine platform, using the LANDTRENDR algorithm segmentation threshold in the Du et al (2022) study for forest disturbance detection.

The acquisition of forest disturbance verification samples includes three sources: (1) Acquiring time and place of occurrence of historical forest fires by using statistical data of forest fires in the detection areas 1986-2023; (2) The method comprises the steps of (1) using samples marked as forests in FAST (https:// doi.org/10.1016/j.scib.2017.03.011) data in a full-ball land coverage fine resolution observation and detection (Finer Resolution Observation and Monitoring of Global Land Cover, FROM-GLC) project as forest disturbance verification samples, (3) randomly selecting the forest disturbance verification samples in a detection area by using a visual interpretation method based on 1986-2023 long time sequence GoogleEarth images.

Table 1 design of experimental protocol

As can be seen from fig. 6 and table 2: (1) The NBR time sequence data synthesized in 1-4 months, which has higher data availability and is frequent by forest disturbance events, can effectively improve the forest disturbance drawing precision in a cloudy and rainy region, and compared with NBR time sequence data in a growing season (6-9 months), F1-Score and OA are improved by 0.5080 and 25.85 percent, and CE and OE are reduced by 43.97 percent and 9.58 percent; (2) The complete and continuous NBR time sequence data in 1986-2023 can be rebuilt by linear interpolation to effectively improve the forest disturbance detection precision of the region with continuous missing Landsat observation values, and compared with the NBR filling the missing year by using a moving average method, the F1-Score and the OA are improved by 0.121 and 10.89 percent, and the CE and the OE are reduced by 15.07 percent and 3.94 percent. (3) The SG filter with constraint can better reserve spectrum information of important time nodes of sudden forest disturbance events, and error detection rate are reduced; the F1-Score and OA were improved by 0.133 and 15.48% and the CE and OE were reduced by 20.61% and 3.18% compared to the original SG filtering algorithm. (4) Compared with the forest disturbance mapping result of LANDTRENDR algorithm, the remote sensing detection result of the complex forest disturbance event can be effectively improved based on the difference between the satellite observation and the spectrum time sequence track after SG filtering fitting, the false detection rate is reduced, the F1-Score and the OA are improved by 0.028 and 2.74%, and the CE and the OE are reduced by 5.27% and 0.30%. In conclusion, compared with the original LANDTRENDR algorithm, the improved LANDTRENDR method has the advantages that the detection precision of the complex forest disturbance event is remarkably improved, F1-Score and OA are improved by 0.8577 and 62.44%, CE and OE are reduced by 92.21% and 22.33%, and the method can effectively reserve complete and reliable spectrum time sequence tracks of various forest disturbance events, reduce the omission ratio and the false detection ratio and realize accurate forest disturbance drafting of complex terrain areas.

Table 2 forest disturbance mapping accuracy evaluation of different experimental schemes

Based on Landsat time sequence data in 1986-2023, the improved LANDTRENDR algorithm is utilized to realize forest disturbance drawing in Yunnan province. From the experimental results, it can be seen that: (1) The forest disturbance in Yunnan province is more severe, the area of the forest disturbance in 1986-2023 is about 1,985,820.9km ², which occupies 49.69% of the total area of the area; wherein, the area of forest disturbance in 1989 is the largest, about 15,498.11 km ², accounting for 3.93% of the total area of the area; the forest disturbance area was minimal in 2022, approximately 315.66 km ², accounting for 0.08% of the total area of the area (fig. 7 a). (2) The forest disturbance in Yunnan province has obvious spatial heterogeneity under the comprehensive influence of factors such as climate, topography, human activity intensity and the like, for example, the forest disturbance year in northwest is mainly concentrated in the 90 th year of the 20 th century, and the change of the surface coverage is mainly changed from forest to cultivated land; forest disturbance in tropical areas such as Xishuangbanna, pu' er, german and the like in the south is mainly concentrated in a large-area planting period (2000-2015), and the surface coverage change mainly changes from a forest to a rubber forest; (3) Thanks to the implementation of national policies such as returning forest, construction of natural protection areas, forest protection utilization and the like, the forest disturbance area is gradually reduced from 2010 in Yunnan province.

From the area-scale forest disturbance area occupancy statistics (fig. 7): (1) The forest disturbance intensities of different elevation areas in Yunnan province have larger difference (figure 7 b), wherein the forest disturbance intensity of the elevation area of 1000-2000m is maximum, about 76.26 percent of the area occupied by the forest disturbance intensity is the largest, and the forest disturbance area occupation ratio of the elevation area of more than 3000m is the smallest, which is 6.44 percent; in addition, the forest disturbance occurrence time of different elevation areas is also greatly different, the disturbance area occupation ratio of the elevation area <1000m area 2007-2012 is highest, and the occupation area is about 17.07%; while forest disturbance in an area with an altitude of more than 1000m is mainly concentrated in 1987-2002, and the disturbance area of the forest disturbance is about 57.96% of the total forest disturbance area of the area. (2) As the annual average temperature of the area rises, the forest disturbance intensity is in a trend of rising and then falling (fig. 7 c), when the annual average temperature is between (6 and 11), the forest disturbance intensity of the area is the largest, and the disturbance area accounts for about 76.42% of the total area of the area; in addition, the occurrence frequency of forest disturbance events in different temperature-uniformization areas is mainly concentrated in 1987-2002, and the forest disturbance area accounts for about 56.63% of the total area of the area disturbance. (3) The forest disturbance intensity and the annual average precipitation amount of the area are in positive correlation (figure 7 d), and when the annual average precipitation amount is more than 122mm, the forest disturbance intensity is maximum and occupies 77.19% of the total area of the area; in addition, the forest disturbance intensities of the regions with precipitation amounts in different years all have a decreasing trend along with the time, wherein the forest disturbance intensity is the largest in 1987-2002 and accounts for about 54.82% of the total area of the forest disturbance of the region. In summary, the forest disturbance intensities in different regions in Yunnan province have obvious spatial heterogeneity, the forest disturbance intensity in regions with the altitude of (2000, 3000 m), the annual average temperature of (6, 11 ℃) and the annual average precipitation of >122mm is the largest, and the disturbance areas respectively occupy 76.26%, 76.42% and 77.19% of the total area of the regions; and the forest disturbance events mainly occur in a concentrated manner in 1987-2002, and the disturbance occurrence area is 56.47% of the total forest disturbance area in the area on average.

According to the technical scheme, the method for detecting the forest disturbance in the complex terrain area effectively reserves real spectrum information of important time nodes of various forest disturbance events, and has better effectiveness for detecting the complex forest disturbance events compared with an original LANDTRENDR algorithm. In the cloud of Yunnan province, the forest disturbance is severe in 1986-2023, and the forest disturbance area accounts for about 49.69% of the total area of the area; forest disturbance events mainly occur intensively in 1986-2002; thanks to the positive implementation of the related forest protection policy, the forest disturbance area in Yunnan province is remarkably reduced from 2002. The natural factors and the humane factors drive the occurrence of the forest disturbance event in Yunnan province together, but natural factors such as climate change, natural disasters, water and soil loss and the like are main driving factors for the frequent occurrence of the forest disturbance event in the region.

The above description is only of the preferred embodiments of the present invention and the description of the technical principles applied is not intended to limit the scope of the invention as claimed, but merely represents the preferred embodiments of the present invention. It will be appreciated by persons skilled in the art that the scope of the invention referred to in the present invention is not limited to the specific combinations of the technical features described above, but also covers other technical features formed by any combination of the technical features described above or their equivalents without departing from the inventive concept. All other embodiments, which can be made by a person skilled in the art without making any inventive effort, are intended to be within the scope of the present invention.

Claims

1. A method for forest disturbance detection, the method comprising the steps of:

（5）

and S7, evaluating the precision of the forest disturbance detection result.

2. The method of claim 1, wherein the dividing process of the forest area and the non-forest area in step S2 is as follows:

3. The method of claim 1, wherein the Landsat series data comprises Landsat Collection tie 1 surface reflectance images available on google earth engine GEE platform for a predetermined month each year of the year.

4. A method of forest disturbance detection according to claim 3, wherein the Landsat Collection t 2t ir 1 surface reflectance image includes: the image collected by three sensors of the thematic drawing instrument TM, the enhanced thematic drawing instrument ETM+, and the land imager OLI.

5. The method for detecting forest disturbance according to claim 4, wherein preprocessing the Landsat series data in step S3 includes:

step S31, performing cloud/shadow/snow masking by utilizing CFmask algorithm;

6. The method for detecting forest disturbance according to claim 4, wherein the calculation formula of the NBR is:

（1）

（2）

7. The forest disturbance detection method according to claim 4, wherein in step S3, the interpolation method is used to fill up the NBR data of the pixels in the missing year, and the NBR time sequence data of the complete detection year of the pixels in the forest area is obtained by reconstruction, specifically including:

（3）

8. The method of forest disturbance detection according to claim 4, wherein step S4 of NBR timing trajectory fitting with constrained SG filtering includes:

9. The method for detecting forest disturbance according to claim 8, wherein the primary NBR timing track is corrected by using an SG filtering algorithm with constraint, and specifically comprising:

（4）

10. A forest disturbance detection system, the system comprising: the system comprises a data acquisition module, a forest area identification module, a first NBR time sequence track generation module, a second NBR time sequence track generation module, a LANDTRENDR algorithm detection module, a year correction module and a precision evaluation module; wherein,

（5）