A forest land landslide identification method based on multi-source time-series remote sensing data

Abstract

A forest landslide identification method based on multi-source time sequence remote sensing data belongs to the technical field of remote sensing image processing and geological disaster monitoring. The forest landslide identification method comprises the following steps: acquiring and preprocessing time sequence Landsat images, land cover data, forest fire data, SPEI and DEM data of a research area; constructing a forest land mask based on early land cover data, and extracting pure forest land time sequence images; adopting a continuous change detection and classification algorithm to perform continuous change detection on NDVI time sequence, and extracting surface disturbance events; then, constructing multi-level non-landslide disturbance removal rules including human factor discrimination, meteorological factor discrimination, topographic factor constraint and NDVI recovery characteristic physical constraint, filtering non-landslide forest land disturbance events, and reserving disturbance events meeting landslide characteristics; and finally outputting a forest landslide spatial distribution map and a time evolution map. The forest landslide identification method can realize rapid and accurate identification of the position and time of forest landslide occurrence in a large range of regions.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing image processing and geological disaster monitoring technology, and in particular relates to a method for identifying forest landslides based on multi-source time-series remote sensing data. Background Technology

[0002] Landslides are a common and serious type of natural disaster in forested mountainous areas. In recent years, affected by global warming, the frequency and intensity of extreme precipitation events have continued to rise, and the risk of landslides in forest areas has also been increasing. Against this backdrop, the rapid and accurate acquisition of the time and spatial location of forest landslides, and the construction of a complete landslide event inventory based on this information, is of great and urgent practical significance for improving regional geological disaster prevention capabilities and protecting the lives and property of the people.

[0003] Traditional landslide identification methods primarily rely on field surveys. While these methods offer high accuracy in localized areas, they generally suffer from significant drawbacks, including high manpower requirements, limited coverage, and difficulty in accessing complex terrain, making them unsuitable for large-scale, high-frequency landslide identification. Remote sensing technology, with its inherent advantages of non-contact observation, wide coverage, and short revisit cycles, has gradually become a crucial technical means for landslide identification and monitoring. Compared to traditional field surveys, remote sensing monitoring can simultaneously capture surface changes before and after landslides on both spatial and temporal scales, significantly improving monitoring efficiency and the comprehensiveness of information acquisition.

[0004] After a landslide, surface vegetation is typically destroyed over a large area, exposing soil and rock. This process is manifested in multi-temporal remote sensing imagery as a sharp decline in the Normalized Difference Vegetation Index (NDVI). Landslide identification methods based on remote sensing technology utilize these anomalous vegetation changes to identify and extract landslide areas. However, the decline in NDVI is not unique to landslides: human factors such as urban expansion and agricultural encroachment on forest land, as well as natural factors such as forest fires and extreme droughts, can also cause significant changes in NDVI. This can easily lead to a high false detection rate in practical applications, limiting the reliability of landslide identification.

[0005] Currently, conventional landslide remote sensing identification methods mainly fall into the following categories: First, relying on high-resolution imagery for manual visual interpretation. While this method offers high accuracy, it is severely limited by manpower and operational efficiency, making it unsuitable for large-scale surveys. Second, feature extraction methods based on single-temporal images. These methods lack the ability to capture dynamic surface processes and struggle to track the evolution of landslides. Third, change detection methods based on dual-temporal or multi-temporal images. Although these methods can reflect surface changes to some extent, their application requires prior knowledge of the approximate time and location of landslides, and the accuracy of the identified landslide occurrence time is directly constrained by the time interval of the selected images. All of these methods have limitations to varying degrees when facing automated landslide identification tasks involving long time series and large areas.

[0006] With the full availability of Landsat time-series remote sensing imagery, monitoring land surface changes using dense time-series remote sensing imagery has become a routine technique. However, its application in forest landslide identification and monitoring remains relatively limited. How to fully leverage the rich spatiotemporal information contained in dense time-series remote sensing imagery to achieve high-precision automated identification of the location and timing of forest landslides is a key technical challenge that urgently needs to be addressed in the field of geological disaster remote sensing monitoring. Summary of the Invention

[0007] To address the shortcomings of existing technologies in the automated identification of long-term, large-scale forest landslides, this invention provides a forest landslide identification method based on multi-source time-series remote sensing data. This method fully integrates the complementary advantages of multi-source time-series remote sensing data, including long-term Landsat satellite imagery, global land cover products, forest fire records, standardized precipitation evapotranspiration index, and digital elevation models. Through dense time-series continuous change detection technology and a multi-level non-landslide disturbance coupling discrimination strategy, it achieves rapid and accurate identification of the location and timing of forest landslides over large areas.

[0008] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0009] A method for identifying forest landslides based on multi-source time-series remote sensing data includes the following steps:

[0010] S1: Acquire temporal Landsat images, land cover data, forest fire data, standardized precipitation evapotranspiration index (SPEI), and DEM data of the study area, and preprocess all data to achieve the unification of spatial resolution and coordinate system of multi-source data.

[0011] S2: Based on the early land cover data of the study area, a forest mask is constructed. The mask is extracted from the preprocessed cloudless optical image to obtain a pure forest time series image. Then, the Continuous Change Detection and Classification (CCDC) algorithm is used to perform continuous change detection on the NDVI time series of the forest image, extract surface disturbance events, and simultaneously obtain the spatial location and time of the disturbance.

[0012] S3: Construct multi-level non-landslide disturbance removal rules by combining multi-source auxiliary data, and sequentially perform coupled screening of four levels: human factor discrimination, meteorological factor discrimination, topographic factor constraint, and NDVI recovery feature physical constraint.

[0013] The physical constraints of the NDVI recovery features are obtained by extracting the NDVI spectral recovery slope after disturbance, and by taking advantage of the characteristics of high slope and slow vegetation recovery in landslide areas, filtering out non-landslide forest disturbance events and retaining forest disturbance events that conform to landslide characteristics.

[0014] S4: The forest disturbance events retained after all levels of screening are identified as landslide candidate points, and the spatial distribution map and temporal evolution map of forest landslides are output.

[0015] Preferably, the long-term Landsat imagery acquired in step S1 is surface reflectance data from all available Landsat series satellites covering the study period, with a spatial resolution of 30m; the land cover data uses a global land cover dataset with an annual update frequency, with a spatial resolution of 30m; the forest fire data uses a global forest fire dataset with an annual update frequency, with a spatial resolution of 30m; the standardized precipitation evapotranspiration index dataset (SPEI) uses a 24-month scale; and the spatial resolution of the DEM data is no less than 30m.

[0016] Preferably, in step S1, the Landsat series satellites include at least one of Landsat5, Landsat7, and Landsat8, with a time span from 2000 to 2022; the global land cover dataset is the GLC_FCS30D dataset; the global forest fire dataset is the UMD_Fire dataset; and the DEM data is SRTM DEM data.

[0017] Preferably, the data preprocessing in step S1 includes a categorized processing procedure, which removes clouds, cloud shadows and oversaturated invalid pixels from Landsat images using quality assessment bands, and resamples the standardized precipitation evapotranspiration index dataset to the same spatial resolution as the Landsat images using projection transformation and nearest neighbor sampling, ensuring that the multi-source data are completely matched in spatial scale.

[0018] Preferably, when constructing the forest mask in step S2, all forest category pixels are extracted from the land cover data of the observation start year to form a mask layer. The forest categories include evergreen broad-leaved forest, evergreen broad-leaved forest-open, deciduous broad-leaved forest, deciduous broad-leaved forest-open, evergreen coniferous forest, evergreen coniferous forest-open, deciduous coniferous forest, deciduous coniferous forest-open, mixed forest, and mixed forest-open. Only cloudless optical images within the forest area are retained for subsequent disturbance detection.

[0019] Preferably, the CCDC algorithm in step S2 achieves change detection by establishing a time series regression model containing trend term and seasonal harmonic term for each pixel. When the deviation between the actual observed value and the model predicted value exceeds 3 times the root mean square error, it is marked as an abnormal observation. After accumulating 6 or more consecutive abnormal observations, it is confirmed as a surface disturbance event. The algorithm uses the LASSO regression method to estimate the model parameters and compresses the unimportant model coefficients to zero through the L1 regularization penalty term to suppress overfitting and improve the generalization ability of time series change detection in response to the uneven spatial distribution and sparsity of Landsat observation data.

[0020] Preferably, the human factor discrimination rule in step S3 is as follows: extract the land cover type from the land cover dataset that is available for the first year after the disturbance occurs and is no more than 3 years after the disturbance occurs. If the type changes to impermeable surface or cultivated land, it is determined to be a disturbance caused by human activity, and the disturbance event is directly removed.

[0021] Preferably, the meteorological factor discrimination rule in step S3 adopts a sequential strategy of fire first and drought later. Specifically, the fire records of the disturbance event in the intermittent year are extracted from the forest fire dataset, and the number of times the pixel where the event is located is marked as a high-confidence fire point is counted (count(UMD_Fire). If count(UMD_Fire)≥1, it is determined to be a disturbance caused by fire and is removed. For the disturbance event that passes the fire screening, the monthly SPEI value in the intermittent year is further extracted from the SPEI dataset, and the frequency of SPEI value less than -1.5 is counted (count(SPEI<-1.5). If count(SPEI<-1.5)≥1, it is determined to be a severe drought disturbance and is removed.

[0022] The topographic factor constraint rule in step S3 is as follows: extract the topographic slope of the study area based on DEM data, set a slope threshold, construct a landslide-prone area mask with a slope greater than the threshold, and eliminate disturbance events in the gentle slope area that do not meet the landslide topographic conditions; the slope threshold is between 10° and 20°, preferably 15°.

[0023] Preferably, the specific implementation of the physical constraint for NDVI recovery features in step S3 is as follows: extract the NDVI time series observations within 1 to 3 years after the disturbance occurs, and calculate the NDVI recovery slope k using least squares linear regression, as shown in equation (2):

[0024] In the formula, t i For the i-th observation time, NDVI i For the corresponding NDVI value, n is the number of valid observations and n≥6. and t i and NDVI i The arithmetic mean of the slope is used; an NDVI recovery slope threshold is set. When the slope is greater than the slope threshold and the NDVI recovery slope k is not greater than the NDVI recovery slope threshold, it is determined to be a landslide disturbance event and retained; otherwise, it is removed. The NDVI recovery slope threshold is between 0.03 / year and 0.08 / year, preferably 0.05 / year.

[0025] In this invention, the determination of the critical threshold is based on the following principles and experimental verification:

[0026] (1) Basis for determining the slope threshold of 15°: Based on the statistical analysis of the slopes of known landslide points in the study area, the slopes of 137 historical landslide points in South China from 2000 to 2022 were statistically analyzed. The slope distribution of landslide points showed a normal distribution characteristic, with a mean of 14.8° and a standard deviation of 3.2°. Selecting 15°, which is close to the mean, as the threshold can cover more than 90% of the landslide points and effectively exclude disturbances in gentle slope areas with a slope of less than 10°. Sensitivity analysis showed that when the slope threshold varied within the range of 10° to 20°, the landslide identification F1 index did not change by more than 5%, indicating that the threshold has good robustness.

[0027] (2) Basis for determining the NDVI recovery slope threshold of 0.05 / year: Based on the time-series comparative analysis of NDVI in typical landslide areas (n=45) and forest logging areas (n=52), the mean NDVI recovery slope in the landslide area was 0.042 / year (standard deviation 0.011) within 3 years after the disturbance, while the mean recovery slope in the forest logging area was 0.087 / year (standard deviation 0.015), and the difference between the two was significant (p<0.01). Selecting 0.05 / year as the threshold can maximize the distinction between the two types of disturbances and achieve the optimal Youden index on the validation set.

[0028] (3) Basis for determining the number of consecutive abnormal observations to 6: Based on the observation density analysis of Landsat time series, the average number of effective observations per year in the study area at 30m resolution is 8-12. Setting the cumulative threshold of consecutive anomalies to 6 can effectively filter out single noise disturbances such as cloud residue and sensor anomalies, while ensuring that the real disturbance events can be confirmed within 3-4 months, avoiding missed detections due to an excessively high threshold.

[0029] Preferably, the method further includes a step of evaluating the accuracy of the recognition results output in step S4. The accuracy evaluation uses user accuracy (UA), producer accuracy (PA), and the F1 index to quantitatively assess the recognition results, wherein:

[0030]

[0031]

[0032]

[0033] In the formula, TP represents a true positive, indicating a correctly identified landslide pixel; FP represents a false positive, indicating that a non-landslide was incorrectly identified as a landslide; and FN represents a false negative, indicating that a landslide was missed and was not identified as a non-landslide.

[0034] Compared with the prior art, the technical solution provided by the present invention has the following beneficial effects:

[0035] First, this invention utilizes dense time-series Landsat images to construct a pixel-by-pixel NDVI time-series model, and uses the CCDC algorithm to achieve continuous automatic detection of surface disturbance events. It does not require pre-setting the time window or spatial search range for landslide occurrence, fundamentally overcoming the inherent defect in traditional dual-temporal or multi-temporal change detection methods where the temporal accuracy is constrained by the image acquisition interval.

[0036] Second, this invention constructs a four-level non-slope disturbance elimination framework that covers human factors, meteorological factors, topographic factors, and NDVI recovery characteristics. It systematically couples and discriminates information on land cover type change, fire records, drought index, topographic slope, and vegetation recovery rate, effectively eliminating major confounding factors such as urban expansion, agricultural encroachment, forest fires, extreme drought, and deforestation, and significantly reducing the false detection rate.

[0037] Third, this invention introduces the NDVI spectral recovery slope after disturbance as a physical constraint condition, and further refines the screening by utilizing the unique characteristics of "high slope and slow vegetation recovery" in landslide areas, effectively distinguishing landslides from disturbance types with similar morphologies such as deforestation.

[0038] Fourth, all data sources in this invention are publicly available remote sensing products, and the overall methodology is highly automated, making it directly applicable to large-scale, long-term forest landslide surveys. In practical verification, this method achieved a user accuracy (UA) of 78.3%, a producer accuracy (PA) of 81.6%, and an F1 index of 79.9% in mountainous forest areas of South China, indicating that the method has high accuracy and reliability in forest landslide identification tasks. Attached Figure Description

[0040] Figure 1 This is a schematic diagram of the overall process of a forest landslide identification method based on multi-source time-series remote sensing data provided by the present invention.

[0041] Figure 2 This is a detailed sub-flowchart of the multi-level non-slope disturbance removal strategy in this invention.

[0042] Figure 3 This is a schematic diagram illustrating the principle of CCDC algorithm time series model fitting and disturbance event detection in this invention.

[0043] Figure 4 This is a comparison chart of NDVI time series between typical landslide areas and non-landslide disturbance areas. Detailed Implementation

[0045] To facilitate understanding and implementation of the present invention by those skilled in the art, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0046] See Figure 1 It can be known that: Figure 1 This diagram illustrates the overall process of a forest landslide identification method based on multi-source time-series remote sensing data provided by the present invention. The diagram, presented as a top-down flowchart, shows the four main steps of the method and their logical relationships.

[0047] The top rectangle in the diagram represents step S1, "Multi-source Remote Sensing Data Acquisition and Preprocessing." The five data sources used in this invention are: Landsat imagery, GLC_FCS30D land cover dataset, UMD_Fire forest fire dataset, SPEI standardized precipitation evapotranspiration index, and SRTM DEM digital elevation model data. This step completes the acquisition of multi-source data and the unification of spatial resolution and coordinate system.

[0048] In spatial resampling and projection transformation processing, the present invention adopts the following unified standards: (1) Coordinate system unification: all data are unified to the WGS84 UTM projection coordinate system, and the corresponding UTM zone is selected according to the longitude zone of the study area. The original Landsat image, GLC_FCS30D, UMD_Fire and SRTM DEM data are all uniformly resampled to this coordinate system. (2) Resampling method selection: different resampling strategies are adopted for different types of remote sensing data: for discrete or categorical data such as Landsat image, land cover classification data, and forest fire data, the nearest neighbor sampling method is adopted to maintain the integrity of the original classification labels and avoid introducing non-existent land types; for DEM data, the bilinear interpolation method is adopted to maintain the continuity of the terrain surface; for the SPEI continuous meteorological index, the nearest neighbor sampling method is adopted to avoid the introduction of false spatial variation by interpolation, while retaining the numerical meaning of the original grid data. (3) Explanation of the rationality of SPEI spatial resampling: As a meteorological drought index, the spatial variation of SPEI is controlled by climatic conditions, and there is no significant spatial heterogeneity at the 30m scale. Resampling the 0.5° grid data to 30m does not introduce new spatial information. Its significance lies in achieving spatial dimensional alignment with Landsat imagery, which facilitates subsequent extraction of the drought index per pixel. The SPEI value of each 30m pixel is inherited from its corresponding 0.5° grid point, ensuring the authenticity of drought intensity information.

[0049] The downward arrow leads to step S2, "Forest Mask Construction and Time Series Change Detection." The box explains the two core operations of this step: first, constructing a forest mask using the forest categories in the GLC_FCS30D dataset to obtain pure forest time series images; second, using the CCDC algorithm to perform continuous change detection on the NDVI time series and extract the set of surface disturbance events.

[0050] The downward arrow leads to step S3, "Multi-level Non-slope Disturbance Removal." This step is represented by a large outer frame containing four sequentially connected sub-frames, corresponding to four filtering levels: The first level is human-induced removal (land cover type change discrimination), with impervious surfaces and farmland marked on the right; the second level is meteorological-induced removal (fire + drought cascade discrimination), with disturbances caused by fire and drought marked on the right; the third level is topographic-related constraints (slope threshold filtering), with disturbances in gentle slope areas marked on the right; and the fourth level is NDVI recovery feature physical constraints (recovery slope filtering), with rapid recovery disturbances marked on the right. The levels are connected sequentially by downward arrows, indicating a progressive filtering process, with the output of each level serving as the input for the next.

[0051] After all four levels of filtering, the results are connected to step S4, "Landslide Identification Result Output," via a down arrow. The box indicates that the output includes a landslide spatial distribution map and a time evolution map. The bottom rectangle is for "Accuracy Evaluation (UA / PA / F1)," which indicates that the final identification result is quantitatively evaluated using user accuracy, producer accuracy, and the F1 index.

[0052] See Figure 2 It can be known that: Figure 2 This is a detailed sub-flowchart of the multi-level non-slope disturbance removal strategy in step S3 of this invention, using standard flowchart symbols to represent the judgment logic. In the diagram, rounded rectangles represent the start and end nodes of the data set, diamonds represent condition judgment nodes, and rectangles represent intermediate data sets or operation results.

[0053] The process begins with the “Forest Disturbance Event Set (from S2)” at the top, which is a collection of all forest disturbance events obtained after detection by the CCDC algorithm in step S2.

[0054] The first diamond-shaped judgment node is "After disturbance, the land cover type changes to impermeable surface or arable land?", corresponding to the first level of human factor discrimination. If the judgment result is "yes", then follow the right arrow to enter the "removal" operation and remove the disturbance event from the set; if the judgment result is "no", then follow the down arrow to enter the intermediate data set "disturbance_data1".

[0055] The second diamond-shaped decision node is "count(UMD_Fire)≥1?", corresponding to the fire factor judgment in the second level. If "yes", the right side outputs "Remove (Fire)"; if "no", it proceeds to the third diamond-shaped decision node "count(SPEI<-1.5)≥1?", corresponding to the drought factor judgment. If "yes", the right side outputs "Remove (Drought)"; if "no", it proceeds to "disturbance_data2". The fire and drought judgment adopts a sequential strategy, first excluding the fire factor, and then assessing the drought impact.

[0056] The fourth diamond-shaped decision node is "slope ≤ 15°?", corresponding to the third-level terrain factor constraint. If "yes", the right side outputs "remove (gentle slope)"; if "no", it proceeds to "disturbance_data3".

[0057] The fifth diamond-shaped judgment node is "NDVI recovery slope k>0.05 / year?", corresponding to the fourth-level NDVI recovery feature physical constraint. If "yes", the right side outputs "reject (fast recovery)", indicating that the vegetation recovery speed of the disturbance event is too fast and does not meet the landslide characteristics; if "no", that is, NDVI recovery slope k≤0.05 / year, it goes down to the bottom rounded rectangle termination node "landslide candidate point (final identification result)", indicating that the disturbance event has passed all four levels of screening and has been identified as a landslide candidate point with high confidence.

[0058] See Figure 3 It can be known that: Figure 3 This is a schematic diagram illustrating the principle of the CCDC (Continuous Change Detection and Classification) algorithm in this invention for time series model fitting and disturbance event detection. The horizontal axis in the figure represents the time axis (in Julian days), and the vertical axis represents the NDVI value, which ranges from 0.2 to 0.8.

[0059] The figure contains the following key graphic elements: Solid dots represent normal observations before the disturbance occurred. These observations are distributed near the fitted curve of the time-series regression model, with NDVI values ​​remaining at a relatively high level between 0.7 and 0.8, reflecting the normal spectral state of healthy vegetation in the forest area. The solid curve represents the time-series regression model fitted by the CCDC algorithm based on historical observation data. This model includes a trend term and a seasonal harmonic term, which can capture the annual seasonal fluctuation pattern of NDVI. Above and below the model curve are short dashed lines forming a 3×RMSE confidence interval, used to define the fluctuation range of normal observations.

[0060] The vertical dashed line in the figure marks the occurrence time of "disturbance events (discontinuity nodes)" (approximately 2015). After the discontinuity nodes, hollow circles represent anomalous observations where the NDVI drops significantly to between 0.3 and 0.5, deviating from the model prediction by more than 3 times the RMSE threshold. When the cumulative number of consecutive anomalous observations reaches 6 or more, the algorithm confirms that a surface disturbance event has occurred.

[0061] After the perturbation event was confirmed, the CCDC algorithm rebuilt a new time-series regression model, represented by a long dashed line. The new model started fitting from the lower NDVI level after the perturbation, reflecting the slow recovery process of vegetation after the perturbation. The legend is located below the figure, labeling four graphic elements: normal observation (solid circle), abnormal observation (hollow circle, >3×RMSE), time-series regression model (solid line), and the new model after the perturbation (long dashed line).

[0062] See Figure 4 It can be known that: Figure 4This is a time-series comparison of NDVI (Non-Landslide Vibration Index) between a typical landslide area and a non-landslide disturbed area (forest logging, post-fire recovery). It visually demonstrates the differences in vegetation recovery rates after different disturbance types, serving as the theoretical basis for the physical constraints of the fourth-level NDVI recovery characteristics in this invention. The horizontal axis represents the time following the disturbance (in years), from before the disturbance to four years after; the vertical axis represents the NDVI value, ranging from 0.0 to 0.8. The moment of disturbance (0 on the horizontal axis) is marked by a vertical dashed line.

[0063] The figure contains three curves with different line types, representing the time series variation characteristics of NDVI under three different disturbance types:

[0064] (1) The solid line represents the NDVI change curve in the landslide area. Before the disturbance, the NDVI remained at a relatively high level of about 0.75. After the disturbance, the NDVI dropped sharply to an extremely low level of about 0.1. Subsequently, within a time window of 1 to 4 years after the disturbance, the NDVI only slowly recovered to about 0.2, with a recovery slope of k≤0.05 / year. The figure is marked with a dashed box "Landslide: Slow vegetation recovery interval, k≤0.05 / year", indicating that the vegetation recovery process in the landslide area is extremely slow due to the severe damage to the surface soil structure and the large-scale stripping of the topsoil.

[0065] (2) The long dashed line represents the NDVI change curve in the forest logging area. Before the disturbance, the NDVI also remained at a high level (about 0.75). After the disturbance, the NDVI also dropped sharply, but within 1 to 2 years after the disturbance, it rebounded significantly and rapidly. By the fourth year, the NDVI had recovered to a high level of about 0.7, close to the pre-disturbance state. The figure is marked with a dashed box "Forest Logging: Rapid Recovery", indicating that the soil structure in the logging area is well preserved, which is conducive to the rapid reconstruction of vegetation.

[0066] (3) The dotted line represents the NDVI change curve of the recovery area after the fire. The NDVI before the disturbance was about 0.68. After the disturbance, the decrease was similar to the previous two. The recovery rate was between that of landslide and deforestation. By the fourth year, the NDVI recovered to about 0.6.

[0067] A comparison of the three curves clearly shows that the NDVI recovery rate of landslide areas after disturbance is significantly lower than that of other disturbance types such as deforestation and post-fire recovery. This invention utilizes this unique characteristic of "high slope and slow vegetation recovery" by setting an NDVI recovery slope threshold (k≤0.05 / year) and combining it with a slope threshold (>15°) to effectively distinguish landslides from other similar forest disturbance types. The legend is located below the figure, labeling four graphic elements: landslide (solid line), deforestation (long dashed line), post-fire recovery (dotted line), and disturbance time (short dashed line).

[0068] See Figures 1-4 As shown, the present invention provides a forest landslide identification method based on multi-source time-series remote sensing data, which can be automated using computer software technology, and mainly includes the following steps:

[0069] Step S1: Acquisition and Preprocessing of Multi-Source Remote Sensing Data

[0070] The core objective of this step is to acquire multi-source remote sensing data covering a sufficient time span within the study area, and to perform standardization preprocessing on data from diverse sources and in different formats to achieve strict alignment in spatial resolution and coordinate reference system, thus establishing a unified data foundation for subsequent analysis. The datasets acquired in this embodiment specifically include the following five categories:

[0071] The first category comprises all available Landsat satellite surface reflectance imagery data from 2000 to 2022, covering Landsat 5, Landsat 7, and Landsat 8 satellites, with a spatial resolution of 30m. Since its launch in 1972, the Landsat series satellites have accumulated decades of continuous Earth observation data. The 16-day revisit cycle provides intensive time-series observations of the study area, ensuring ample data for capturing gradual and abrupt changes in surface vegetation. This data is available through the U.S. Geological Survey (USGS) Earth Explorer platform.

[0072] The second category is the temporal global land cover dataset GLC_FCS30D, spanning from 2000 to 2022, with a spatial resolution of 30m and annual updates. This dataset provides detailed land cover classification information, covering various land cover types such as forests, cultivated land, and impervious surfaces, playing a crucial role in forest boundary delineation and land cover change assessment. This dataset is available through the Zenodo data sharing platform.

[0073] The third category is the time-series forest fire dataset UMD_Fire, spanning from 2000 to 2022, with a spatial resolution of 30m and annual updates, released by the University of Maryland. This dataset records forest fire events globally and their confidence levels, and can be used to identify and rule out forest vegetation changes caused by fires. This dataset is available through the GLAD Labs.

[0074] The fourth category is the Standardized Precipitation-Evapotranspiration Index (SPEI) dataset, spanning from 2000 to 2022, with a spatial resolution of 0.5° and multiple time scales available. This invention employs the 24-month SPEI index, which effectively reflects the cumulative impact of medium- to long-term drought conditions on vegetation growth, thus identifying vegetation decline caused by severe drought events. This dataset can be obtained through the SPEI global drought monitoring system.

[0075] The fifth category is SRTM (Digital Elevation Model) data, with a spatial resolution of 30m. DEM data is the fundamental data source for extracting terrain slope information and plays a crucial role in subsequent slope threshold selection. This data can be obtained through the USGS Earth Explorer platform.

[0076] In the data preprocessing stage, this invention adopts differentiated processing strategies for different types of data. For Landsat imagery data, pixels contaminated by clouds and cloud shadows are identified and removed using the Quality Assessment (QA) bands inherent in each image. Simultaneously, oversaturated and invalid pixels are filtered out based on reasonable ranges of surface reflectance for each band, ensuring that all observations used in subsequent analysis are high-quality clear-sky data. For the SPEI dataset, since its original spatial resolution is 0.5° (approximately 55km), which differs significantly from the 30m resolution of the Landsat imagery, projection transformation and spatial resampling are required. Specifically, the nearest neighbor sampling method is used to resample the SPEI dataset to a 30m resolution, making it perfectly match the Landsat imagery in spatial scale. The nearest neighbor sampling method is chosen because SPEI is a continuous meteorological index characterizing drought severity; this sampling method ensures that the original values ​​are not altered by the interpolation algorithm, thus guaranteeing the accuracy of drought intensity information.

[0077] Step S2: Woodland Mask Construction and Time Series Change Detection

[0078] This step includes two core components: first, constructing a forest mask for the study area to strictly limit the analysis scope to the forest region; second, using a continuous change detection algorithm to perform automated change detection on the NDVI time series of the forest region, extracting all surface disturbance events and their corresponding spatiotemporal information.

[0079] In the forest mask construction stage, this invention extracts all forest category pixels from the global land cover dataset GLC_FCS30D, which is based on the observation start year (2000 in this embodiment). The classification codes corresponding to the forest categories in the GLC_FCS30D dataset include: 51 (evergreen broad-leaved forest), 52 (evergreen broad-leaved forest - open), 61 (deciduous broad-leaved forest), 62 (deciduous broad-leaved forest - open), 71 (evergreen coniferous forest), 72 (evergreen coniferous forest - open), 81 (deciduous coniferous forest), 82 (deciduous coniferous forest - open), 91 (mixed forest), and 92 (mixed forest - open). After merging the pixels corresponding to the above classification codes, a unified forest mask layer is formed. Subsequently, this mask is used to spatially crop the preprocessed temporal Landsat image from step S1, retaining only cloudless optical image pixels located within the forest area. The reason for using land cover data from the starting year of the observation to construct the mask is to ensure that it reflects the original state of forest cover before the disturbance event, and to prevent potential landslide areas from being missed due to changes in the forest area later.

[0080] In the time-series change detection stage, this invention selects NDVI as the core spectral index characterizing vegetation status. NDVI is one of the most widely used vegetation indices in the field of remote sensing, with values ​​typically ranging from -1 to 1. Healthy and dense vegetation areas correspond to higher NDVI values, while bare soil, rocks, and other non-vegetated areas correspond to lower NDVI values. When a landslide occurs in a forest area, the vegetation suffers large-scale damage and is accompanied by exposed ground, resulting in a sudden and continuous decrease in NDVI values. This significant spectral response characteristic provides a reliable physical basis for detecting landslide events using time-series methods. Figure 3 As shown, the CCDC algorithm can accurately locate the time node when NDVI undergoes a sudden change by fitting a time series model to a historical observation period, thereby capturing the precise occurrence time of the disturbance event.

[0081] This invention employs the CCDC algorithm to continuously detect changes in the NDVI time series of each forest pixel. CCDC is a classic dense time series change detection method. Its basic principle is to build a time series regression model for each pixel based on historical observation data, which includes a trend term and a seasonal harmonic term, and use this model to predict subsequent observations. When the deviation between the subsequent actual observation value and the model prediction value exceeds a set threshold, a surface disturbance event is determined to have occurred, and the algorithm then reconstructs a new time series model after the event. The mathematical expression of the CCDC change detection model is shown in equation (1):

[0082] (1)

[0083] In equation (1), Let be the model prediction value for the i-th band of the t-th Julian Day; Let be the regression intercept of the i-th band; The regression slope of the i-th band reflects the long-term trend of NDVI. and , respectively, are the cosine and sine coefficients of the k-th harmonic term in the i-th band, used to fit the annual seasonal fluctuation pattern of NDVI; n is the order of the harmonic term; T is the number of days in a year, usually 365. In this embodiment, the harmonic order n is set to 3, that is, the third harmonic is used to fully characterize the seasonal variation of vegetation.

[0084] In determining disturbance events, CCDC employs a strategy of accumulating consecutive outlier observations: when the difference between a new actual observation and the model prediction exceeds three times the root mean square error (RMSE), the observation is marked as an outlier; if six or more consecutive outlier observations occur, a surface disturbance event is confirmed, and the time corresponding to the first outlier observation is recorded as the disturbance occurrence time. The cumulative number of consecutive outliers is set to six, a default threshold widely used by the CCDC algorithm in Landsat time series scenarios. This value maintains sufficient sensitivity to real disturbance events while avoiding false detections due to random factors such as single-observation noise or cloud remnants.

[0085] Due to the spatial unevenness of Landsat observation data—the number of effective observations varies significantly across different regions, and the distribution of effective observations over time is sparse and discontinuous due to factors such as cloud cover—this invention employs the LASSO (Least Absolute Shrinkage and Selection Operator) regression method for model parameter estimation to improve the model's generalization ability and effectively suppress overfitting. LASSO regression, by introducing an L1 regularization penalty term, automatically compresses unimportant model coefficients to zero during parameter estimation, achieving robust parameter estimation results even with insufficient observation data, thereby enhancing the generalization ability of time-series change detection.

[0086] In this invention, the key parameters of the CCDC algorithm adopt publicly available default values. Those skilled in the art can directly repeat the implementation based on these parameters without creative effort. Specifically, `minObservations` represents the minimum number of effective observations required to start the time series model, with a default value of 6; `chiSquareProbability` is the chi-square probability threshold for controlling change detection sensitivity, with a default value of 0.99; and `minNumOfYearScaler` is the minimum time span scaling factor required for modeling, with a default value of 0.5. These parameters can be adaptively adjusted according to the cloud cover and observation density of the study area: when there are sufficient effective observations in the study area, `minObservations` can be increased to 8 to enhance model stability; when it is necessary to increase change detection sensitivity, `chiSquareProbability` can be decreased to 0.95; in long-term, densely observed scenarios, `minNumOfYearScaler` can be adjusted to 1.0 to ensure the robustness of the time series model.

[0087] After processing by the CCDC algorithm, the NDVI time series of each forest pixel is divided into several time periods, and the discontinuities between adjacent segments are the detected surface disturbance events. This invention fully records the spatial coordinates and occurrence time of each disturbance event, thereby forming a preliminary set of forest disturbance events for subsequent in-depth discrimination and analysis.

[0088] Step S3: Multi-level non-landslide disturbance removal

[0089] The forest disturbance event set obtained in step S2 contains, in addition to landslide events, a large number of disturbances caused by non-landslide factors such as human activities, meteorological disasters, and topographical constraints. This step constructs a systematic, multi-level set of non-landslide disturbance removal rules, such as... Figure 2 As shown, by filtering layer by layer, the aforementioned confounding factors are eliminated sequentially, ultimately retaining landslide candidate points with high credibility. The multi-level elimination rule designed in this invention includes four levels: human factor elimination, meteorological factor elimination, topographic factor constraint, and NDVI recovery feature physical constraint. The specific implementation methods of each level are as follows:

[0090] Level 1: Human Factor Identification and Removal. For each forest disturbance event obtained in step S2, the land cover type of the disturbance event in the first available year after the disturbance occurred (i.e., the year following the year of the disturbance; if the data for the following year is unavailable, it is extended to the most recent available year, but no more than 3 years after the disturbance occurred) is extracted from the time-series land cover dataset GLC_FCS30D. If the land cover type of the pixel changes from forest to impervious surface (GLC_FCS30D classification code 190) or cultivated land (classification codes 10, 11, 12, 20), the disturbance event is determined to be caused by urban expansion or agricultural reclamation activities and is removed. After this level of filtering, the set of retained forest disturbance events is denoted as disturbance_data1. The basis for this distinction is that after a landslide, the surface is generally covered with bare soil or gravel, which will not be classified as impermeable surface or arable land in land cover classification; while the reduction of forests due to urbanization and agricultural encroachment will clearly transform the land cover into construction land or arable land. Therefore, it can be distinguished by the direction of change in land cover type after disturbance.

[0091] Level 2: Meteorological Factor Identification and Removal. This level performs a two-stage concatenated discrimination for each disturbance event in disturbance_data1, using forest fire data and drought index data respectively. First, fire records for the disturbance event within the year of the disturbance (i.e., the year of the discontinuity node detected by CCDC, hereinafter referred to as the discontinuity year) are extracted from the UMD_Fire dataset, and the number of times the pixel where the event is located is marked as a high-confidence fire point is counted, count(UMD_Fire). If count(UMD_Fire)≥1, that is, there is at least one high-confidence fire record for the disturbance event within the discontinuity year, the disturbance event is determined to be caused by a forest fire and is removed. For disturbance events that pass the fire screening (i.e., count(UMD_Fire)<1), the drought discrimination stage begins: the monthly SPEI values ​​for the disturbance event within the discontinuity year are extracted from the 24-month scale SPEI dataset, and the number of months with SPEI values ​​less than -1.5 is counted, denoted as count(SPEI<-1.5). When the SPEI value is below -1.5, the corresponding area is considered to be in a state of severe drought. If count(SPEI<-1.5)≥1, the disturbance event is determined to be caused by vegetation decline due to extreme drought stress and is therefore excluded. After this level of screening, the set of forest disturbance events retained is denoted as disturbance_data2. This level adopts a cascaded discrimination strategy of fire first and then drought because fire discrimination has higher confidence and discrimination efficiency. Prioritizing the exclusion of fire factors before assessing the impact of drought is beneficial to improving the overall screening efficiency and reducing the risk of misjudgment.

[0092] The third level: Topographic constraints. This level utilizes DEM data to extract topographic slope information of the study area and constructs a landslide-prone area mask using a slope threshold of 15°. Specifically, areas with slopes greater than 15° are retained as landslide-prone areas, while disturbance events in disturbance_data2 located in gentle slope areas with slopes less than or equal to 15° are removed. This threshold setting is based on the basic understanding of geological hazards: a certain topographic slope is a necessary condition for landslides to occur. Overly gentle terrain is unlikely to accumulate sufficient gravitational potential energy to drive the sliding of soil and rock masses; therefore, vegetation disturbance events detected in gentle slope areas are highly unlikely to be caused by landslides. After this level of screening, the set of retained forest disturbance events is denoted as disturbance_data3.

[0093] Level 4: Physical Constraints Based on NDVI Recovery Characteristics. This level further incorporates the spectral recovery slope of NDVI after disturbance to construct physical constraints, allowing for a more refined distinction between landslides and other similar forest disturbance types. The principle is that after a landslide, the surface soil structure is severely damaged, with large-scale stripping and accumulation of the topsoil, resulting in extremely slow vegetation recovery. The recovery rate of NDVI after disturbance is significantly lower than other disturbance types such as deforestation and post-fire recovery. For example... Figure 4 As shown, the NDVI in the landslide area remained at a low level for several years after the disturbance, while the NDVI in the logging area showed a significant rebound within 1 to 2 years. For the disturbance events in disturbance_data3 that were filtered by terrain, the NDVI time-series observations within 1 to 3 years after the disturbance were extracted, and the NDVI recovery slope k was calculated using least squares linear regression. The calculation formula is shown in equation (2):

[0094] (2)

[0095] In equation (2), For the i-th observation time (in Julian days), Let be the NDVI observation value corresponding to that moment, and n be the number of valid observations within a time window of 1 to 3 years after the perturbation, with n ≥ 6 required to ensure the statistical reliability of the regression results. and They are respectively and The arithmetic mean of the values ​​is used. An NDVI recovery slope threshold of 0.05 / year is set. When the slope is greater than 15° and the NDVI recovery slope k ≤ 0.05 / year, it is identified as a landslide disturbance event and retained; otherwise, it is identified as a non-landslide disturbance such as logging or rapid post-fire recovery and is removed. This invention utilizes the unique characteristics of "high slope and slow vegetation recovery" in landslide areas to further filter out non-landslide forest disturbances, improving the reliability of identification.

[0096] After the above four levels of screening, the forest disturbance events that are finally retained are the landslide candidate sites with high credibility.

[0097] Step S4: Landslide identification results output and accuracy evaluation

[0098] The landslide candidate points selected in step S3 are used as the final identification results and visualized. A spatial distribution map of forest landslides is generated based on the spatial coordinate information of each landslide candidate point, which intuitively shows the geographical distribution pattern of landslide events in the study area. At the same time, a temporal evolution map of forest landslides is generated based on the occurrence time of each disturbance event recorded by the CCDC algorithm, revealing the occurrence pattern and evolution trend of landslide events in the time dimension.

[0099] To objectively evaluate the recognition performance of the method of this invention, three classic classification accuracy evaluation indicators are used to quantify the results: User Accuracy (UA), Producer Accuracy (PA), and F1 index. The calculation formulas for the three are shown in equations (3), (4), and (5), respectively:

[0100] (3)

[0101] (4)

[0102] (5)

[0103] In the above formula, TP (True Positive) represents the number of pixels correctly identified as landslides; FP (False Positive) represents the number of non-landslide pixels incorrectly identified as landslides; and FN (False Negative) represents the number of landslide pixels missed. User accuracy (UA) measures the proportion of actual landslides among all pixels identified as landslides, reflecting the precision of the identification results; producer accuracy (PA) measures the proportion of successfully identified actual landslide pixels, reflecting the recall of the identification results; and the F1 index is the harmonic mean of UA and PA, comprehensively considering both precision and recall to fully evaluate the overall quality of landslide identification results.

[0104] In practical accuracy verification, existing landslide survey lists or concurrent high-resolution remote sensing images can be used as reference ground truth data. Spatial overlay analysis is used to compare the identification results with the reference data one by one, and the values ​​of TP, FP, and FN are statistically analyzed and then substituted into the above formulas to calculate each accuracy index. In this embodiment, a typical mountainous forest area in South China is used as the verification area, and a manually interpreted high-resolution Google Earth image landslide list is used as reference data. Statistical calculations show that the user accuracy (UA) of this method is 78.3%, the producer accuracy (PA) is 81.6%, and the F1 index is 79.9%, indicating that this method has good overall performance in distinguishing landslides from other forest disturbances.

[0105] It should be understood that any parts not described in detail in this specification belong to the prior art. The above description of the preferred embodiments is quite detailed, but it should not be considered as a limitation on the scope of protection of this invention. Those skilled in the art, under the guidance of this invention, can make substitutions or modifications without departing from the scope of protection of the claims of this invention, and all such substitutions or modifications fall within the scope of protection of this invention. The scope of protection of this invention should be determined by the appended claims.

Claims

1. A method for identifying forest landslides based on multi-source time-series remote sensing data, characterized in that, Includes the following steps: S1: Acquire temporal Landsat images, land cover data, forest fire data, standardized precipitation evapotranspiration index (SPEI), and DEM data of the study area, and preprocess all data to achieve the unification of spatial resolution and coordinate system of multi-source data; S2: Based on the early land cover data of the study area, a forest mask is constructed. The mask is extracted from the preprocessed cloudless optical image to obtain a pure forest time series image. Then, the CCDC algorithm is used to continuously detect changes in the NDVI time series of the forest image, extract surface disturbance events, and simultaneously obtain the spatial location and time of the disturbance. S3: Construct multi-level non-landslide disturbance removal rules by combining multi-source auxiliary data, and sequentially perform coupled screening of four levels: human factor discrimination, meteorological factor discrimination, topographic factor constraint, and physical constraint based on NDVI recovery features, to retain forest disturbance events that meet the characteristics of landslides; S4: The forest disturbance events retained after all levels of screening are identified as landslide candidate points, and the spatial distribution map and temporal evolution map of forest landslides are output.

2. The forest landslide identification method based on multi-source time-series remote sensing data according to claim 1, characterized in that: The long-term Landsat imagery acquired in step S1 is surface reflectance data from all available Landsat series satellites covering the study period, with a spatial resolution of 30m; the land cover data uses the global land cover dataset with an annual update frequency, with a spatial resolution of 30m; the forest fire data uses the global forest fire dataset with an annual update frequency, with a spatial resolution of 30m; the standardized precipitation evapotranspiration index dataset (SPEI) uses a 24-month scale; and the spatial resolution of the DEM data is no less than 30m.

3. The forest landslide identification method based on multi-source time-series remote sensing data according to claim 2, characterized in that: In step S1, the Landsat series satellites include at least one of Landsat 5, Landsat 7, and Landsat 8, with a time span from 2000 to 2022; the global land cover dataset is the GLC_FCS30D dataset; the global forest fire dataset is the UMD_Fire dataset; and the DEM data is SRTM DEM data.

4. The forest landslide identification method based on multi-source time-series remote sensing data according to claim 1 or 2, characterized in that: Step S1 involves data preprocessing, which includes a categorized processing flow. For Landsat images, cloud, cloud shadows, and oversaturated invalid pixels are removed using quality assessment bands. For the standardized precipitation evapotranspiration index dataset, projection transformation and nearest neighbor sampling are used to resample it to the same spatial resolution as the Landsat images, ensuring that the multi-source data are completely matched in spatial scale.

5. The forest landslide identification method based on multi-source time-series remote sensing data according to claim 1, characterized in that: In step S2, when constructing the forest mask, all forest category pixels are extracted from the land cover data of the observation start year to form a mask layer. The forest categories include evergreen broad-leaved forest, evergreen broad-leaved forest-open, deciduous broad-leaved forest, deciduous broad-leaved forest-open, evergreen coniferous forest, evergreen coniferous forest-open, deciduous coniferous forest, deciduous coniferous forest-open, mixed forest, and mixed forest-open. Only cloudless optical images within the forest area are retained for subsequent disturbance detection.

6. The forest landslide identification method based on multi-source time-series remote sensing data according to claim 5, characterized in that: In step S2, the CCDC algorithm detects changes by establishing a time series regression model for each pixel that includes a trend term and a seasonal harmonic term. When the deviation between the actual observed value and the model prediction value exceeds 3 times the root mean square error, it is marked as an abnormal observation. After accumulating 6 or more consecutive abnormal observations, it is confirmed as a surface disturbance event. The algorithm uses the LASSO regression method to estimate the model parameters and compresses unimportant model coefficients to zero through the L1 regularization penalty term to suppress overfitting and improve the generalization ability of time series change detection, taking into account the uneven spatial distribution and sparsity of Landsat observation data.

7. The forest landslide identification method based on multi-source time-series remote sensing data according to claim 1, characterized in that: The human factor discrimination rule in step S3 is as follows: extract the land cover type from the land cover dataset that is available for the first year after the disturbance and is no more than 3 years after the disturbance. If the type changes to impermeable surface or cultivated land, it is determined to be a disturbance caused by human activity and the disturbance event is directly removed.

8. The forest landslide identification method based on multi-source time-series remote sensing data according to claim 7, characterized in that: The meteorological factor discrimination rule in step S3 adopts a sequential strategy of first fire and then drought. Specifically, the fire records of the disturbance event in the intermittent year are extracted from the forest fire dataset. The number of times the pixel where the event is located is marked as a high-confidence fire point is counted (UMD_Fire). If count(UMD_Fire)≥1, it is determined to be a disturbance caused by fire and is removed. For disturbance events that pass the fire screening, the monthly SPEI values ​​for the intermittent years are further extracted from the SPEI dataset, and the frequency of SPEI values ​​less than -1.5 is counted (SPEI<-1.5). If count(SPEI<-1.5)≥1, it is judged as a severe drought disturbance and is removed.

9. The forest landslide identification method based on multi-source time-series remote sensing data according to claim 8, characterized in that: The topographic factor constraint rule in step S3 is as follows: extract the topographic slope of the study area based on DEM data, set a slope threshold, construct a landslide-prone area mask with a slope greater than the threshold, and eliminate disturbance events in the gentle slope area that do not meet the landslide topographic conditions; the slope threshold is between 10° and 20°.

10. The forest landslide identification method based on multi-source time-series remote sensing data according to claim 1 or 9, characterized in that: The specific implementation method of the physical constraint for NDVI recovery features in step S3 is as follows: extract the NDVI time series observations within 1 to 3 years after the disturbance occurs, and calculate the NDVI recovery slope k using least squares linear regression: In the formula, t i For the i-th observation time, NDVI i For the corresponding NDVI value, n is the number of valid observations and n≥6. and t i and NDVI i The arithmetic mean of the slope is used; an NDVI recovery slope threshold is set. When the slope is greater than the slope threshold and the NDVI recovery slope k is not greater than the NDVI recovery slope threshold, it is determined to be a landslide disturbance event and retained; otherwise, it is removed. The NDVI recovery slope threshold is between 0.03 / year and 0.08 / year.

Images