A method for large-scale forest disturbance detection and disturbance attribution analysis

Abstract

The present invention belongs to the field of forest disturbance monitoring, and provides a large-scale forest disturbance detection and disturbance attribution analysis method, including: obtaining a multi-year forest cover composite map, optimizing using a pixel compensation model, generating a cover composite map, improving the basic CCDC algorithm, extracting disturbances, calculating disturbance data, calculating feature data, data fusion, calculating associated data, and determining the target attribution result. By optimizing and predicting the multi-year forest cover composite map, the present invention improves the quality and quantity of images used for analysis, and improves the classification accuracy and reliability; by improving the basic CCDC algorithm, the improvement of the disturbance extraction effect is realized; by calculating feature data and associated data, intuitive data is provided for disturbance attribution analysis, and the efficiency of disturbance attribution analysis is improved.

Description

Technical Field

[0001] The present invention relates to the field of forest disturbance monitoring, and particularly to a method for large-scale forest disturbance detection and disturbance attribution analysis. Background Art

[0002] With the increase of global climate change and human activities, the frequency and intensity of forest disturbances have increased significantly, seriously affecting the stability of forest ecosystems. Existing forest disturbance monitoring methods mostly rely on ground observations and traditional remote sensing data, suffering from problems such as poor data timeliness, limited spatial coverage, and inability to effectively capture large-scale disturbance dynamics.

[0003] The existing technologies mainly include using medium-resolution MODIS and AVHRR data for coarse-grained disturbance monitoring, and using Landsat data for fine-grained change detection. Although these methods have solved the problem of spatial resolution to a certain extent, there is still much room for improvement in data processing capabilities, algorithm accuracy, and monitoring timeliness, and there is also a lack of comprehensive analysis of multiple factors. Summary of the Invention

[0004] In order to overcome the deficiencies of the prior art, the purpose of the present invention is to provide a method for large-scale forest disturbance detection and disturbance attribution analysis, which improves the quality and quantity of images for analysis, as well as the classification accuracy and reliability, by optimizing and predicting multi-year forest cover composite maps; improves the disturbance extraction effect by improving the basic CCDC algorithm; and provides intuitive data for disturbance attribution analysis and improves the efficiency of disturbance attribution analysis through the calculation of feature data and associated data.

[0005] To achieve the above purpose, the present invention provides the following solutions:

[0006] A method for large-scale forest disturbance detection and disturbance attribution analysis, comprising:

[0007] Collect and fuse long-term Landsat images and Google Maps historical images to obtain a multi-year forest cover composite map;

[0008] Input the multi-year forest cover composite map into a pre-trained pixel compensation model for optimization to obtain a high-precision forest cover composite map;

[0009] Input the high-precision forest cover composite map into a pre-trained neural network model for prediction to obtain a generated cover composite map;

[0010] Add the disturbance maximum amplitude band and ID mapping value to the basic CCDC algorithm to obtain an improved CCDC algorithm;

[0011] The improved CCDC algorithm is used to extract the disturbance year, the maximum disturbance amplitude, and the number of disturbances from the high-precision forest cover composite map and the generated cover composite map, respectively, to obtain a forest disturbance map; the forest disturbance map includes: a year data map, a maximum amplitude data map, and a number data map;

[0012] Extracting the disturbance spatial distribution, disturbance density, and change trend of the forest disturbance map to obtain disturbance data;

[0013] Performing discriminant analysis on the disturbance data using preset factors to obtain forest fire attribution maps, human attribution maps, and other attribution maps; the human attribution maps include: logging attribution maps, farming attribution maps, and construction activity attribution maps; the other attribution maps include: climate attribution maps and disaster attribution maps;

[0014] Fusing the forest fire attribution map, the human factor attribution map, and the other attribution maps to obtain a fused attribution map;

[0015] The global Moran index, local Moran index and Getis-Ord Gi* index are used to analyze the fusion attribution map to obtain the target attribution result.

[0016] Preferably, the calculation formula of the change trend includes:

[0017]

[0018] Wherein, θ is the linear change slope of the trend; FDD i is the forest disturbance density in year i; n is the length of time.

[0019] Preferably, the training process of the pixel compensation model includes:

[0020] collecting a basic forest cover composite map and a fine forest cover composite map; the fine forest cover composite map having a higher precision than the basic forest cover composite map;

[0021] Extracting pixel information of the fine forest cover composite image to obtain an original pixel dataset;

[0022] Extracting high-precision pixel information of the multi-year forest cover composite map to obtain a target pixel dataset;

[0023] Build a knn classification model;

[0024] The original pixel data set is used as input data, the target pixel data set is used as verification data, the knn classification model is trained, and the pixel compensation model is obtained.

[0025] Preferably, the types of the maximum amplitudes of the perturbations include: low amplitude, medium-low amplitude, medium amplitude, medium-high amplitude, and high amplitude; the range of the low amplitude is: 0 to 0.09; the range of the medium-low amplitude is: 0.09 to 0.19; the range of the medium amplitude is: 0.19 to 0.30; the range of the medium-high amplitude is: 0.30 to 0.45; the range of the high amplitude is: greater than 0.45.

[0026] Preferably, the preset factors include: climate factors, human factors, and other factors; the other factors include: insect outbreaks and river course changes.

[0027] Preferably, the calculation formula of the global Moran's I index is:

[0028]

[0029] where I is the global Moran's I index; n is the number of disturbed patches; x i and x j are the attribute values of the i-th and j-th disturbed patches respectively; w i,j is the spatial weight between the i-th and j-th disturbed patches; is the average of the attribute values.

[0030] Preferably, the calculation formula of the local Moran's I index includes:

[0031]

[0032] where I i is the local Moran's I index.

[0033] Preferably, the calculation formula of the Getis-Ord Gi* index includes:

[0034]

[0035] where is the i-th Getis-Ord Gi* index.

[0036] The present invention discloses the following technical effects:

[0037] The present invention provides a method for large-scale forest disturbance detection and disturbance attribution analysis. By optimizing and predicting multi-year forest cover composite maps, the problems of low accuracy and small quantity of image sources are solved, and the improvement of the quality and quantity of images for analysis is achieved; by improving the basic CCDC algorithm, the problem of poor disturbance extraction effect of the existing CCDC algorithm is solved, and the improvement of the disturbance extraction effect is realized; by calculating feature data and correlation data, the problem of complex conventional disturbance attribution analysis process is solved, and intuitive feature data is provided for disturbance attribution analysis. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0039] Figure 1 Schematic diagram of the large-scale forest disturbance detection and disturbance attribution analysis process provided by the embodiment of the present invention;

[0040] Figure 2 Schematic diagram of the pixel compensation model training process provided by the embodiment of the present invention;

[0041] Figure 3 Schematic diagram of the k-nearest neighbor principle provided by the embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0042] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, rather than all embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts belong to the scope of protection of the present invention.

[0043] The object of the present invention is to provide a method for large-scale forest disturbance detection and disturbance attribution analysis. By optimizing and predicting multi-year forest cover composite maps, the quality and quantity of images for analysis are improved, and the classification accuracy and reliability are enhanced; by improving the basic CCDC algorithm, the improvement of the disturbance extraction effect is realized; by calculating feature data and correlation data, intuitive data is provided for disturbance attribution analysis, and the efficiency of disturbance attribution analysis is increased.

[0044] In order to make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below with reference to the drawings and specific embodiments.

[0045] Figure 1 Schematic diagram of the large-scale forest disturbance detection and disturbance attribution analysis process provided by the embodiments of the present invention, as Figure 1 shown, the present invention provides a large-scale forest disturbance detection and disturbance attribution analysis method, including:

[0046] Step 100: Collect and fuse long-term Landsat images and Google Maps historical images to obtain a multi-year forest cover composite map;

[0047] Step 200: Input the multi-year forest cover composite map into a pre-trained pixel compensation model for optimization to obtain a high-precision forest cover composite map;

[0048] Step 300: Input the high-precision forest cover composite map into a pre-trained neural network model for prediction to obtain a generated cover composite map;

[0049] Step 400: Add the disturbance maximum amplitude band and ID mapping values to the basic CCDC algorithm to obtain an improved CCDC algorithm;

[0050] Step 500: Use the improved CCDC algorithm to extract the disturbance year, disturbance maximum amplitude, and disturbance times from the high-precision forest cover composite map and the generated cover composite map respectively to obtain a forest disturbance map; the forest disturbance map includes: a year data map, a maximum amplitude data map, and a times data map;

[0051] Step 600: Extract the disturbance spatial distribution, disturbance density, and change trend from the forest disturbance map respectively to obtain disturbance data;

[0052] Step 700: Use preset factors to perform discriminant analysis on the disturbance data respectively to obtain a forest fire attribution map, a human factor attribution map, and other attribution maps; the human factor attribution map includes: a logging attribution map, an agricultural cultivation attribution map, and a construction activity attribution map; other attribution maps include: a climate attribution map and a disaster attribution map;

[0053] Step 800: Fuse the forest fire attribution map, the human factor attribution map, and other attribution maps to obtain a fused attribution map;

[0054] Step 900: Analyze the fused attribution map using the global Moran's I index, local Moran's I index, and Getis-Ord Gi* index to obtain the target attribution result.

[0055] Specifically, the CCDC algorithm (Change Detection and Classification) is an algorithm used in the fields of remote sensing and image processing, mainly for detecting and classifying changes in images or time series data. This algorithm is usually used to monitor environmental changes, land use changes, urban development, deforestation, etc. The core steps of the CCDC algorithm generally include: Data preprocessing: Correct, denoise, standardize, etc. the input remote sensing data to ensure the quality and usability of the data. Time series analysis: Use time series remote sensing data to analyze the changes between different time points. By introducing temporal information, the changing areas can be identified more effectively. Change detection: Use statistical methods or machine learning techniques to detect the areas that have changed in the time series. Methods such as threshold method, classification method, regression method, etc. can be used. Change classification: Classify the detected changes, such as classifying the changes into types such as urban construction, farmland changes, forest reduction, etc. Result verification and analysis: Verify and analyze the change detection and classification results by comparing with existing ground truth data or more detailed remote sensing data to improve the accuracy and reliability. The advantage of the CCDC algorithm is that it can process large-scale remote sensing data sets and can extract effective information from multi-temporal data, which is suitable for the monitoring and evaluation of dynamic changes. It has important application value for analysis and decision support (such as resource management, environmental monitoring, etc.). In this embodiment, the CCDC algorithm can segment the time series data of pixels, divide the entire time series into several sub-time series, and each sub-time series corresponds to a time series model. Based on these time series, the number and length (in years) of line segments, the slope, and the year corresponding to the vertex can be obtained.

[0056] Furthermore, the data set source of this embodiment is referred to Table 1.

[0057] Table 1

[0058]

[0059]

[0060] Preferably, the calculation formula for the change trend includes:

[0061]

[0062] where θ is the slope of the trend linear change; FDD i is the forest disturbance density in the i-th year; n is the time length.

[0063] Reference Figure 2 , the training process of the pixel compensation model includes:

[0064] Step 201: Collect the basic forest cover composite map and the detailed forest cover composite map; the accuracy of the detailed forest cover composite map is higher than that of the basic forest cover composite map.

[0065] Step 202: Extract the pixel information of the detailed forest cover composite map to obtain the original pixel dataset.

[0066] Step 203: Extract the pixel information of the high-precision multi-year forest cover composite map to obtain the target pixel dataset.

[0067] Step 204: Construct a knn classification model.

[0068] Step 205: Use the original pixel dataset as the input data and the target pixel dataset as the validation data to train the knn classification model to obtain a pixel compensation model.

[0069] Preferably, the types of maximum perturbation amplitudes include: low amplitude, medium-low amplitude, medium amplitude, medium-high amplitude, and high amplitude; the range of low amplitude is: 0 to 0.09; the range of medium-low amplitude is: 0.09 to 0.19; the range of medium amplitude is: 0.19 to 0.30; the range of medium-high amplitude is: 0.30 to 0.45; the range of high amplitude is: greater than 0.45.

[0070] Optionally, the preset factors include: climate factors, human factors, and other factors; other factors include: insect outbreaks and river course changes.

[0071] Specifically, the calculation formula of the global Moran's I is:

[0072]

[0073] where I is the global Moran's I; n is the number of disturbance patches; x i 、x j are the attribute values of the i-th and j-th disturbance patches respectively; w i,j is the spatial weight between the i-th and j-th disturbance patches; is the average of the attribute values.

[0074] Optionally, the Global Moran's I is a statistic used in spatial data analysis, especially for detecting spatial autocorrelation. It was proposed by statistician Patrick A. Moran in 1950. The Global Moran's I is used to evaluate the distribution pattern of a variable in space to determine whether the variable shows characteristics of clustering, dispersion, or random distribution.

[0075] Preferably, the calculation formula of the local Moran's I includes:

[0076]

[0077] Among them, I i is the local Moran's I index.

[0078] Optionally, the Local Moran's I is a statistic used for spatial data analysis. It is an extension of the global Moran's I and is used to detect the spatial autocorrelation of a specific geographical location. The Local Moran's I can help identify the clustering or distribution characteristics at specific locations in spatial data, especially suitable for identifying hotspots (areas with high-value clusters) and cold spots (areas with low-value clusters).

[0079] Specifically, the calculation formula of the Getis-Ord Gi* index includes:

[0080]

[0081] Among them, is the i-th Getis-Ord Gi* index.

[0082] Optionally, the Getis-Ord Gi* index (also known as Getis-Ord hot spot analysis) is a statistic used for spatial data analysis, aiming to detect local clustering patterns in spatial data. This index can identify "hot spots" (areas with high-value clusters) and "cold spots" (areas with low-value clusters) and evaluate whether these clustering phenomena are significant.

[0083] Preferably, the obtained multi-year forest cover composite map is input into the pixel compensation model. Each sub-pixel point corresponds to a set of color weights, and the top four color ratios with the highest proportion are selected for combination as the optimized color of the sub-pixel point.

[0084] Specifically, kNN (k-Nearest Neighbor) is a basic and simple classification algorithm. As a supervised learning method, the KNN model requires labeled training data. The class of a new sample is determined by the k training sample points closest to the new sample according to the classification decision rule. The k-nearest neighbor (kNN) method is a basic classification and regression method; it is a model based on labeled training data; it is a supervised learning algorithm. The three key points of the basic approach are: First, determine the distance metric; second, select the value of k (find the k instance points in the training set that are closest to the point to be estimated); third, the classification decision rule. In classification tasks, the "voting method" can be used, that is, select the label category that appears most frequently among these k instances as the prediction result; in regression tasks, the "average method" can be used, that is, take the average of the real-value outputs of these k instances as the prediction result; weighted average or weighted voting can also be performed based on the distance, and the closer the instance, the greater the weight. The k-nearest neighbor method does not have an explicit learning process. It is a famous representative of lazy learning. Such learning techniques only save the samples during the training stage, with zero training time overhead, and then process them after receiving the test samples.

[0085] Reference Figure 3 , the main methods of k-nearest neighbor clustering analysis are as follows: Initially, randomly assign K (hyperparameter) cluster centers, and these K cluster centers are randomly selected from the sample points. Divide the sample points to be classified into each cluster according to the nearest neighbor principle. After the division, recalculate the centroid of each cluster by the average method to determine the new cluster center. Keep iterating until the moving distance of the cluster center is less than a given value. The evaluation criterion of the algorithm is the sum of squared errors criterion. When the sum of squared errors reaches the optimal (small) value, it can make the intra-class as compact as possible and the clusters as separated as possible. The loss function is a non-convex function, so there are many minimum values, and it may not fall into the global optimum every time, and may fall into the local optimum. Use the loss function for global optimization during iteration. The remaining data is grouped into clusters according to the silhouette coefficient; the calculation formula of the silhouette coefficient is:

[0086]

[0087] where s is the silhouette coefficient; b is the sum of the distances between the data to be assigned and all points in the next nearest cluster; a is the sum of the distances between the data to be assigned and other points in the same cluster; max() is the maximum value function. The range of the silhouette coefficient is (-1, 1), where the value closer to 1 means that the sample is very similar to the samples in its own cluster and not similar to the samples in other clusters. When the sample point is more similar to the sample points outside the cluster, the silhouette coefficient is negative. When the silhouette coefficient is 0, it means that the two clusters have the same similarity and the two clusters should have been one cluster.

[0088] First, the official CCDC code was used for disturbance monitoring, but there were some problems in extracting the year with the largest change. For example, in the entire research period, the NDVI of a certain pixel changed three times, and the change years were 1995, 2005, and 2018, with the change amplitudes being -0.0833, -0.3828, and -0.3237 respectively. According to the absolute value of the change rate, the year with the largest change was 2005, but the CCDC algorithm extracted 2018. After checking the underlying source code of CCDC, it was found that when the CCDC algorithm calculated the year of the largest change breakpoint, it did not use the "magnitude" band, did not use ID mapping to obtain values, but determined tBreak according to the year of the latest breakpoint. Therefore, in this study, the CCDC code was modified, and the disturbance maximum amplitude band and ID mapping value were added to the basic CCDC algorithm to obtain the improved CCDC algorithm.

[0089] Furthermore, several parameters of the CCDC algorithm need to be adjusted. The value range of the parameters was determined by visually inspecting the known disturbance sample points (Table 2 shows the set values of the value range in this embodiment). The fitting degree of the trajectories under different parameter combinations with the reality was visually evaluated to determine the best parameter combination across samples for each forest area. Based on the CCDC algorithm, a forest disturbance map was output. The disturbance map includes the year of the largest disturbance, the number of disturbances, and the maximum disturbance amplitude: the year with the largest disturbance amplitude during the detection period, that is, for a given date range and spectral band, the time when the largest detected change amplitude occurs is visualized; the number of disturbances, that is, the total number of detected disturbances within the specified time period; the maximum disturbance amplitude, that is, the largest change value detected within the specified time period and spectral band, which is measured by the difference between the end point and the start point of adjacent time periods (the slope of the fitting curve, and the larger the slope, the more severe the disturbance).

[0090] Table 2

[0091]

[0092]

[0093] Specifically, the maximum disturbance amplitude detected by the CCDC algorithm is the slope of the fitting curve. The larger the slope (mag) of the maximum disturbance amplitude, the more severe the disturbance. According to the natural breakpoints, mag is divided into 5 categories: low (0 < mag ≤ 0.09), medium - low (0.09 < mag ≤ 0.19), medium (0.19 < mag ≤ 0.30), medium - high (0.30 < mag ≤ 0.45), and high (mag > 0.45).

[0094] Specifically, the driving factors of forest disturbance mainly include logging, agriculture (deforestation caused by slash - and - burn agriculture and commercial agricultural expansion), forest fires, urbanization, and other disturbance factors (such as insect outbreaks, wind disasters, ice disasters, floods, or river course changes).

[0095] Preferably, for the discrimination of forest fire disturbances: overlay the GABAM product with the forest disturbance area, and conduct conditional discrimination pixel by pixel and year by year to extract forest fire disturbances. The specific discrimination method is as follows: if the disturbance occurrence year of a certain pixel is the same as the fire occurrence year, it is determined as a forest fire disturbance area; otherwise, it is other disturbances.

[0096] Furthermore, for the discrimination of human disturbances: assume that a certain pixel is disturbed in the t-th year, and the (t - 1)-th year in the CLCD is forest, but: if the (t + 1)-th and (t + 2)-th years are shrubs, grasslands or bare lands, it is determined as logging; if the (t + 1)-th and (t + 2)-th years are impervious, it is determined as construction; if the (t + 1)-th and (t + 2)-th years are cultivated lands, it is determined as farming.

[0097] Specifically, the value range of I is [-1, 1]. I > 0 indicates positive correlation, I = 0 indicates no correlation, and I < 0 indicates negative correlation. This tool evaluates and expresses clustering, discrete or random patterns. For example, z-score or p-value represents statistical significance. If I > 0, it is clustering; otherwise, it is discrete.

[0098] Furthermore, the natural break method is used to visually express the Getis-Ord Gi* index in different periods. The Gi_Bin field will identify the cold and hot spots of statistical significance. The pixels with confidence intervals of +3 to -3, +2 to -2, and +1 to -1 are the statistical significances with confidence levels of 99%, 95%, and 90% respectively; while the confidence interval of 0 has no statistical significance.

[0099] The beneficial effects of the present invention are as follows:

[0100] By optimizing and predicting the multi-year forest cover composite map, the present invention improves the quality and quantity of the images used for analysis, and improves the classification accuracy and reliability; by improving the basic CCDC algorithm, the improvement of the disturbance extraction effect is realized; through the calculation of feature data and associated data, intuitive data is provided for the disturbance attribution analysis, and the disturbance attribution analysis efficiency is improved.

[0101] In this specification, each embodiment is described in a progressive manner. The key point of each embodiment is to illustrate the differences from other embodiments. The same or similar parts among the embodiments can be referred to each other.

[0102] In this article, specific examples are used to elaborate on the principle and implementation manner of the present invention. The descriptions of the above embodiments are only used to help understand the method and its core idea of the present invention; at the same time, for those of ordinary skill in the art, based on the idea of the present invention, there will be changes in the specific implementation manner and application scope. In summary, the content of this specification should not be construed as a limitation to the present invention.

Claims

1. A method for large-scale forest disturbance detection and disturbance attribution analysis, characterized in that including: Collecting and fusing long-term Landsat images and Google Maps historical images to obtain a multi-year forest cover composite map; Inputting the multi-year forest cover composite map into a pre-trained pixel compensation model for optimization to obtain a high-precision forest cover composite map; Inputting the high-precision forest cover composite map into a pre-trained neural network model for prediction to obtain a generated cover composite map; Adding the maximum amplitude band of perturbation and ID mapping values to the basic CCDC algorithm to obtain an improved CCDC algorithm; Using the improved CCDC algorithm to extract the perturbation year, maximum perturbation amplitude, and number of perturbations from the high-precision forest cover composite map and the generated cover composite map respectively to obtain a forest perturbation map; the forest perturbation map includes: a year data map, a maximum amplitude data map, and a number data map; Extracting the perturbation spatial distribution, perturbation density, and change trend from the forest perturbation map respectively to obtain perturbation data; Performing discriminant analysis on the perturbation data using preset factors respectively to obtain a forest fire attribution map, a human factor attribution map, and other attribution maps; the human factor attribution map includes: a logging attribution map, an agricultural cultivation attribution map, and a construction activity attribution map; the other attribution maps include: a climate attribution map and a disaster attribution map; Fusing the forest fire attribution map, the human factor attribution map, and the other attribution maps to obtain a fused attribution map; Analyzing the fused attribution map using the global Moran's I index, local Moran's I index, and Getis-Ord Gi* index to obtain a target attribution result.

2. The large-scale forest disturbance detection and disturbance attribution analysis method according to claim 1, wherein The calculation formula of the change trend includes: where θ is the slope of the linear change of the trend; FDD i is the forest disturbance density in the i-th year; n is the time length.

3. A large-scale forest disturbance detection and disturbance attribution analysis method according to claim 1, characterized in that The training process of the pixel compensation model includes: Collecting a basic forest cover composite map and a fine forest cover composite map; the accuracy of the fine forest cover composite map is greater than that of the basic forest cover composite map; Extracting the pixel information of the fine forest cover composite map to obtain an original pixel data set; Extracting the pixel information of the high-precision multi-year forest cover composite map to obtain a target pixel data set; Constructing a knn classification model; Using the original pixel data set as input data and the target pixel data set as verification data to train the knn classification model to obtain the pixel compensation model.

4. A large-scale forest disturbance detection and disturbance attribution analysis method according to claim 1, characterized in that The types of the maximum perturbation amplitude include: low amplitude, medium-low amplitude, medium amplitude, medium-high amplitude, and high amplitude; the range of the low amplitude is: 0 to 0.09; the range of the medium-low amplitude is: 0.09 to 0.19; the range of the medium amplitude is: 0.19 to 0.30; the range of the medium-high amplitude is: 0.30 to 0.45; the range of the high amplitude is: greater than 0.

45.

5. A method for large-scale forest disturbance detection and disturbance attribution analysis according to claim 1, characterized in that, The preset factors include: climate factors, human factors, and other factors; the other factors include: insect outbreaks and river course changes.

6. The large-scale forest disturbance detection and disturbance attribution analysis method according to claim 2, characterized in that The calculation formula of the global Moran's I index is: where I is the global Moran's index; n is the number of disturbed patches; x i , x j are the attribute values of the i-th and j-th disturbed patches respectively; w i,j is the spatial weight between the i-th and j-th disturbed patches; is the average of the attribute values.

7. A method for large-scale forest disturbance detection and disturbance attribution analysis according to claim 6, characterized in that The calculation formula of the local Moran's I index includes: Among them, I i is the local Moran's index.

8. A method for large-scale forest disturbance detection and disturbance attribution analysis according to claim 7, characterized in that The calculation formula of the Getis-Ord Gi* index includes: Among them, is the i-th Getis-Ord Gi* index.

Images