A method for identifying damaged land

Abstract

The application relates to the technical field of environmental monitoring, and provides a damaged land identification method, which comprises the following steps: acquiring long-time sequence remote sensing image data according to the spatial range and the monitoring period of a region to be investigated; generating a vegetation index curve covering the whole-year time sequence based on the long-time sequence remote sensing image data; the vegetation index curve is used for representing the annual vegetation growth phenology characteristics in the latest state of the current annual land surface; performing traversal analysis on the vegetation index curve, and respectively determining an annual maximum value index and an annual minimum value index; the time period corresponding to the annual maximum value index represents the vegetation growth season, and the time period corresponding to the annual minimum value index represents the non-vegetation growth season; determining an annual amplitude vegetation index according to the annual maximum value index and the annual minimum value index; the annual amplitude vegetation index is used for quantifying the growth fluctuation amplitude of the current annual land surface vegetation in the latest state; detecting an interannual change trajectory, correcting the abnormal value of the interannual change trajectory; and identifying damaged land according to the corrected change trajectory.

Description

Technical Field

[0001] This application relates to the field of environmental monitoring technology, and in particular to a method for identifying damaged land. Background Technology

[0002] In recent years, with the continuous acceleration of economic construction and resource development activities, the issues of land resource protection and management have become increasingly prominent. Land is not only the basic carrier of agricultural production, but also an important component of the ecological environment system. A comprehensive understanding of the current state of land damage is of great significance for scientifically formulating land resource protection policies, implementing ecological restoration projects, and promoting sustainable use.

[0003] According to statistics, in 2019, the area of ​​land damaged by mining nationwide reached 3.6105 million hectares, of which approximately 1.3404 million hectares were damaged by mines under construction and approximately 2.2701 million hectares were damaged by abandoned mines. Currently, the problem of unclear land inventory and unknown figures still exists in my country. Existing historically damaged land has not yet been restored, and newly damaged land continues to increase, which has become a major long-term challenge for land reclamation work. There is an urgent need to establish an efficient and accurate technology for identifying damaged land in order to achieve rapid monitoring and dynamic updating of the scope and temporal changes of damage. Summary of the Invention

[0004] In a first aspect, this application provides a method for identifying damaged land; comprising: acquiring long-term remote sensing image data based on the spatial range and monitoring period of the area to be investigated, wherein the long-term remote sensing image data includes blue light, green light, red light, near-infrared, and short-wave infrared bands; generating a vegetation index curve covering the entire year's time series based on the long-term remote sensing image data; the vegetation index curve being used to characterize the annual vegetation growth phenological characteristics under the latest surface conditions of the current year; performing traversal analysis on the vegetation index curve to determine the annual maximum value index and the annual minimum value index respectively; the period corresponding to the annual maximum value index represents the peak vegetation growth season, and the period corresponding to the annual minimum value index represents the non-growing vegetation season; determining the Intra-Annual Amplitude Vegetation Index based on the annual maximum value index and the annual minimum value index. IAAVI); the annual amplitude vegetation index is a relative quantity formed by the difference between the annual maximum value index and the annual minimum value index, used to quantify the growth fluctuation amplitude of the current year's surface vegetation in its latest state; detect the interannual IAAVI change trajectory, identify and correct outliers in the interannual IAAVI change trajectory; the interannual IAAVI change trajectory is a change trajectory of IAAVI values ​​including multiple years; identify damaged land based on the corrected interannual IAAVI change trajectory.

[0005] In some possible implementations, generating a vegetation index curve covering the entire year's time series based on the long-term remote sensing image data includes: selecting a vegetation index type; performing band calculations on each image based on the long-term remote sensing image data to obtain a vegetation index time series for the vegetation index type, wherein the vegetation index is used to indicate the land cover type and vegetation change characteristics; fitting the vegetation index time series of each pixel to obtain the temporal change trajectory of the corresponding pixel; performing breakpoint detection based on the temporal change trajectory; determining the annual fitting coefficient based on the breakpoint detection results; and generating a vegetation index curve covering the entire year's time series based on the temporal change trajectory and the annual fitting coefficient, wherein the vegetation index curve is used to characterize the annual vegetation growth phenological characteristics under the latest land surface conditions of the current year.

[0006] In some possible implementations, the step of performing band-wise calculations on the long-term remote sensing image data to obtain vegetation index values ​​includes: selecting vegetation index types based on the type, stage, and surface characteristics of land damage-reclamation in the area to be investigated; the vegetation index types include normalized vegetation index, enhanced vegetation index, comprehensive forest characteristic index, and / or normalized degradation fraction index; and performing band-wise calculations on the long-term remote sensing image data to determine the vegetation index values ​​corresponding to the vegetation index types.

[0007] In some possible implementations, the step of detecting breakpoints based on the temporal change trajectory and determining annual fitting coefficients based on the breakpoint detection results includes: identifying breakpoints in the temporal change trajectory, where the breakpoints indicate the transition of a pixel from one state to another in the time series; extracting the occurrence time of the breakpoint, the start and end times of the segments before and after the breakpoint, and the fitting coefficients; using historical years as the time unit, if there are no breakpoints within a year, the annual fitting coefficients are taken from the fitting coefficients of the unique segment corresponding to that year; if there are breakpoints within a year, the annual fitting coefficients are updated to the fitting coefficients of the last segment; the fitting coefficients are used to characterize the latest state of the land surface after the abrupt event at the breakpoint.

[0008] In some possible implementations, the vegetation index curve includes breakpoints, and the step of performing traversal analysis on the vegetation index curve to determine the annual maximum and minimum values ​​respectively further includes: correcting the annual maximum and minimum values ​​of the representative remote sensing index based on the updated annual fitting coefficient, wherein the annual maximum and minimum values ​​include the annual maximum and minimum values.

[0009] In some possible implementations, the step of identifying damaged land based on the annual amplitude vegetation index change trajectory includes: detecting the annual amplitude vegetation index change trajectory and identifying "pseudo-breakpoints"; the "pseudo-breakpoints" are outliers that are unrelated to actual surface changes; correcting the outliers to obtain the corrected interannual IAAVI change trajectory; and analyzing each pixel in the area to be investigated one by one based on the corrected interannual IAAVI change trajectory to identify damaged land.

[0010] In some possible implementations, identifying damaged land based on the annual amplitude vegetation index change trajectory includes: using the IAAVI change trajectory of all pixels in the area to be investigated during the monitoring period as input samples, employing the Otsu automatic thresholding method, and determining a damage determination threshold by automatically searching for the optimal segmentation point based on the principle of maximizing inter-class variance; the damage determination threshold is used to determine the surface damage status of each pixel in different years; and determining whether damage has occurred and the year of damage occurrence based on the damage determination threshold and the IAAVI change trajectory.

[0011] In some possible implementations, the method further includes: continuing to perform interannual IAAVI variation trajectory analysis on the damaged pixels in the area to be investigated to determine whether reclamation has been implemented and to determine the year of implementation of reclamation; if the determination result is that reclamation has been implemented, assessing the state after reclamation to determine whether the reclamation is sufficient; and generating a spatial distribution map of damaged land based on the land damage and reclamation status of each pixel.

[0012] In some possible implementations, the method further includes: classifying the land in the area to be investigated into types based on the harmonic fitting trajectory of the remote sensing index and the interannual IAAVI variation trajectory, combined with the damage and reclamation status and the adequacy of reclamation, including the following types: undamaged land; damaged land that has not been reclaimed; damaged land that has not been fully reclaimed; and damaged land that has been fully reclaimed.

[0013] The method provided in this application effectively solves problems such as interannual trend interference and low sensitivity to disturbance identification by using the "Intra-annual Amplitude Vegetation Index (IAAVI)" indicator.

[0014] The method provided in this application is based on the dynamic update of annual fitting coefficients driven by breakpoints, which corrects the annual extreme value index of the representative remote sensing index, ensuring that the remote sensing index can characterize the latest ecological state after the current year's surface disturbance, and effectively solves the "time lag" problem in the disturbance time point identification of the MVC method.

[0015] The method provided in this application is for outlier detection and correction of interannual IAAVI variation trajectory, which can systematically identify and correct "pseudo-breakpoints" generated during the CCDC model fitting process, thereby improving the accuracy of damaged land identification.

[0016] The "National Survey Implementation Plan for Land Damaged by Mining (Trial)" emphasizes the need to comprehensively ascertain the current status and baseline of land damaged by mining nationwide (including excavation, occupation, and subsidence), and to establish a dynamic data update and sharing mechanism. The method provided in this application's embodiments can achieve large-scale, automated, and time-series identification and monitoring of damaged land, quickly and accurately obtaining information such as land disturbance distribution, area changes, and correction progress, providing technical support for natural resources departments to establish a dynamic monitoring system that is "clearly identified, accurately calculated, and effectively managed."

[0017] The "National Green Mine Construction Evaluation Indicators" includes "establishing a dynamic monitoring system to dynamically monitor geological environment damage and restoration, land damage and reclamation, and ecosystem damage (degradation) and restoration" as one of the evaluation indicators. The method provided in this application's embodiments enables continuous monitoring and time-series analysis of the land surface. It can not only identify the time and extent of damage events but also track the effects of reclamation and restoration, providing quantitative evidence for green mine construction assessment and ecological restoration performance evaluation.

[0018] Besides mining areas, land disturbance issues are also prevalent in various scenarios such as urban expansion, transportation construction, water conservancy projects, and natural disasters. The method provided in this application is universal and can be widely applied to various tasks of identifying and monitoring changes in land surface disturbances. By continuously tracking changes in vegetation phenological characteristics throughout the year, it can identify abnormal fluctuations and potential degradation trends in the land surface, providing a scientific basis for land use planning, post-disaster assessment, and ecological risk early warning.

[0019] Secondly, this application provides a damaged land identification device, including a data acquisition module that acquires long-term remote sensing image data based on the spatial range and monitoring period of the area to be investigated. The long-term remote sensing image data includes blue, green, red, near-infrared, and short-wave infrared bands. A curve generation module generates a vegetation index curve covering the entire year's time series based on the long-term remote sensing image data; the vegetation index curve is used to characterize the annual vegetation growth phenological characteristics under the latest state of the land surface. A curve analysis module performs a comprehensive analysis of the vegetation index curve to determine the annual maximum and minimum values; the period corresponding to the annual maximum value represents the peak growth season of vegetation, and the period corresponding to the annual minimum value represents the non-growing season of vegetation. An IAAVI determination module determines the annual amplitude vegetation index based on the annual maximum and minimum values; the annual amplitude vegetation index is a relative quantity formed by the difference between the annual maximum and minimum values, used to quantify the growth fluctuation amplitude of the current year's land surface vegetation under the latest state. The damage identification module detects the interannual IAAVI change trajectory, identifies and corrects outliers in the interannual IAAVI change trajectory; the interannual IAAVI change trajectory is the trajectory of IAAVI value changes in multiple years; and the damaged land is identified based on the corrected interannual IAAVI change trajectory.

[0020] Thirdly, this application provides a computing device, comprising: at least one memory for storing a program; and at least one processor for executing the program stored in the memory, wherein when the program stored in the memory is executed, the processor is configured to execute the method as described in any one of the first aspects.

[0021] Fourthly, this application provides a computer-readable storage medium storing a computer program that, when run on a processor, causes the processor to perform the method as described in any one of the first aspects. Attached Figure Description

[0022] To more clearly illustrate the technical solutions of the various embodiments disclosed in this specification, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only a few embodiments disclosed in this specification. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] The accompanying drawings used in the description of the embodiments or prior art are briefly introduced below.

[0024] Figure 1 A flowchart illustrating the damaged land identification method provided in this application embodiment;

[0025] Figure 2This is a schematic diagram of the time-series variation trajectory after harmonic fitting provided in an embodiment of this application;

[0026] Figure 3 This diagram illustrates the time lag problem in damage monitoring using the MVC method.

[0027] Figure 4 Schematic diagram of annual amplitude vegetation index under different time segments provided in the embodiments of this application;

[0028] Figure 5 This is a schematic diagram illustrating the classification of damaged land types in the damaged land identification method provided in the embodiments of this application;

[0029] Figure 6 This is a schematic diagram illustrating the identification of the year of mine damage in the method provided in the embodiments of this application;

[0030] Figure 7 This is a schematic diagram illustrating the classification results of damaged land types in the mining area in the method provided in the embodiments of this application;

[0031] Figure 8 A damaged land identification device provided in this application embodiment;

[0032] Figure 9 A computing device provided in an embodiment of this application. Detailed Implementation

[0033] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions in the embodiments of this application will be described below with reference to the accompanying drawings.

[0034] In the description of the embodiments of this application, the words "exemplary," "for example," or "for instance" are used to indicate examples, illustrations, or explanations. Any embodiment or design described as "exemplary," "for example," or "for instance" in the embodiments of this application should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the words "exemplary," "for example," or "for instance" is intended to present the relevant concepts in a specific manner.

[0035] In the description of the embodiments in this application, the term "and / or" is merely a description of the association relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, B existing alone, and A and B existing simultaneously. Furthermore, unless otherwise stated, the term "multiple" means two or more. For example, multiple systems refer to two or more systems, and multiple terminals refer to two or more terminals.

[0036] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. The terms "comprising," "including," "having," and their variations all mean "including but not limited to," unless otherwise specifically emphasized.

[0037] In the description of the embodiments in this application, "some embodiments" are mentioned, which describe a subset of all possible embodiments. However, it is understood that "some embodiments" can be the same subset or different subsets of all possible embodiments, and can be combined with each other without conflict.

[0038] In the description of the embodiments of this application, the terms "first, second, third, etc." or module A, module B, module C, etc. are used only to distinguish similar objects and do not represent a specific ordering of objects. It is understood that, where permitted, a specific order or sequence can be interchanged so that the embodiments of this application described herein can be implemented in an order other than that illustrated or described herein.

[0039] In the description of the embodiments of this application, the reference numerals for the steps, such as S110, S120, etc., do not necessarily indicate that the steps will be executed in this manner. Where permissible, the order of the steps can be interchanged or executed simultaneously.

[0040] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. The terminology used herein is for the purpose of describing embodiments of this application only and is not intended to limit this application.

[0041] Traditional methods for identifying and monitoring damaged land typically employ a comparative approach after classifying "dual-temporal" or "multi-temporal" images. By comparing remote sensing images from different periods, changes in land use types, the direction of land type transitions, and spatial distribution characteristics can be identified. While these methods can reveal spatial changes in land cover, they are significantly generalized in their representation of the change process over time, making it difficult to continuously depict the dynamic evolution of land damage and restoration.

[0042] In recent years, with the development of high-resolution remote sensing data such as domestic high-resolution satellites, many scholars have begun to use high-resolution images from a limited number of periods combined with manual visual interpretation or deep learning to achieve high-precision identification of damaged land in local areas. This type of method has high resolution on a spatial scale, but it is difficult to achieve large-scale monitoring due to factors such as image acquisition frequency, cloud interference, and data cost. Moreover, the image coverage time is limited, and there is a lack of traceable data to support the damage and reclamation process in historical periods.

[0043] With the long-term accumulation and free access to medium- and high-resolution remote sensing data such as Landsat and Sentinel, time-series remote sensing analysis methods have been widely applied in the identification of damaged land. On the one hand, some scholars only select single high-quality images for analysis each year, or use filtering algorithms to smooth the selected images to construct land cover change trajectories for subsequent damaged land identification. On the other hand, with the continuous development of remote sensing cloud computing platforms, most studies utilize annual multi-imagery, generating annual images through the Maximum Value Composites (MVC) method using vegetation indices such as NDVI and EVI, and then combining filtering and curve fitting smoothing algorithms to construct land cover change trajectories, thereby identifying the land damage and restoration process. Based on the MVC method, time-series change detection algorithms for annual composite data have also been derived, such as the LandTrendr algorithm.

[0044] Such methods improve the continuity of time series to a certain extent, but have the following problems: (1) If the land damage occurs after the vegetation growth season of the year, that is, damage occurs in autumn or winter, the annual images of the MVC strategy directly applied cannot reflect the damage status of the year, resulting in misjudgment or lag in the identification of damage time; (2) Due to the differences between different sensors and the neglect of the temporal continuity between years in the data selection process, the change trajectory constructed based on single high-quality images or the MVC method often has significant temporal instability. Even if the land cover status is relatively stable and no significant disturbance occurs, the trajectory may still fluctuate abnormally, affecting the accurate expression of the land cover change trend.

[0045] This application provides a method for identifying damaged land, which improves upon the above-mentioned approach of directly applying the MVC strategy for annual damaged land identification.

[0046] This application provides a method for identifying damaged land. It utilizes a remote sensing cloud computing platform to acquire long-term remote sensing image data and optimizes the Continuous Change Detection and Classification (CCDC) algorithm to extract pixel-level time series features, thereby achieving accurate and rapid identification of damaged land. This effectively overcomes the problems of delayed identification of damage time and insufficient temporal stability in traditional methods during change detection.

[0047] Figure 1 A flowchart illustrating the damaged land identification method provided in this embodiment. Figure 1 As shown, it includes:

[0048] S11. Based on the spatial range and monitoring period of the area to be investigated, long-term remote sensing image data is acquired. The long-term remote sensing image data includes blue light, green light, red light, near-infrared, and short-wave infrared bands.

[0049] S12 generates a vegetation index curve covering the entire year's time series based on long-term remote sensing image data; the vegetation index curve is used to characterize the annual vegetation growth phenology characteristics under the latest surface conditions of the current year.

[0050] S13. Perform a traversal analysis on the vegetation index curve to determine the maximum and minimum values ​​for the year. The period corresponding to the maximum value represents the peak growth season for vegetation, while the period corresponding to the minimum value represents the non-growing season for vegetation.

[0051] S14. Determine the annual amplitude vegetation index (IAAVI) based on the annual maximum and annual minimum indicators. The annual amplitude vegetation index is a relative quantity formed by the difference between the annual maximum and annual minimum indicators, and is used to quantify the growth fluctuation of the current year's surface vegetation in its latest state.

[0052] S15 detects the interannual IAAVI variation trajectory, identifies and corrects outliers; the interannual IAAVI variation trajectory includes the variation trajectory of the annual amplitude vegetation index values ​​for multiple years.

[0053] S16, Identify and classify damaged land based on the revised interannual IAAVI variation trajectory.

[0054] The above steps will be further explained below with reference to specific implementation methods and accompanying drawings.

[0055] For S11, in some possible implementations, the area to be investigated and the monitoring period can be determined according to the actual situation and user needs, and long-term remote sensing image data can be obtained using a remote sensing cloud computing platform based on the spatial range of the area to be investigated and the monitoring period.

[0056] For example, in the Google Earth Engine (GEE) platform, all available Landsat TM / ETM+ / OLI / OLI-2 and Sentinel-2 MSI surface reflectance image data of the area to be investigated within the monitoring period can be obtained in chronological order of image capture time. Each image must include blue, green, red, near-infrared, shortwave infrared 1 and shortwave infrared 2 bands. The Landsat image processing level should be Collection 2 Level-2, and the Sentinel-2 image processing level should be Level 2A.

[0057] In some possible implementations, long-term remote sensing image data can be preprocessed, including cloud cover filtering, cloud and shadow removal, image cropping, and image mosaicking.

[0058] For example, the acquired Landsat TM / ETM+ / OLI / OLI-2 and Sentinel-2 MSI surface reflectance image data can be preprocessed in the GEE platform, including cloud cover filtering, cloud and shadow removal, image cropping and image mosaicking, to ensure the quality and consistency of remote sensing images in subsequent time series analysis.

[0059] For S12, in some possible implementations, step S12 includes the following steps S121-S124:

[0060] S121. Select vegetation index types based on the type, stage, and surface characteristics of land damage-reclamation in the area to be investigated. Vegetation index types include normalized vegetation index, enhanced vegetation index, comprehensive forest characteristic index, and / or normalized degradation fraction index.

[0061] S122, Perform band-wise calculations on the long-term remote sensing image data to determine the corresponding type of vegetation index time series. The long-term remote sensing image data includes preprocessed long-term remote sensing image data.

[0062] For example, the normalized difference vegetation index (NDVI) will be used as an example. The normalized difference vegetation index N... DVI It is an important indicator reflecting the surface vegetation cover and its productivity. It is one of the most widely used vegetation indices, with advantages such as high sensitivity in vegetation detection. Its calculation formula is as follows:

[0063] (1)

[0064] In equation (1), R NIR R RED These represent the reflectance values ​​in the near-infrared and red light bands, respectively.

[0065] S123, the vegetation index time series of each pixel is fitted to obtain the temporal change trajectory of the corresponding pixel.

[0066] In some possible implementations, to effectively address the temporal instability problem, a Continuous Change Detection and Classification (CCDC) algorithm can be employed in the GEE platform to perform harmonic function fitting on the observed vegetation index time series values ​​of each pixel within the survey area. The fitted harmonic function is then used as the formula to express the temporal change trajectory of the corresponding pixel. The harmonic function fitting formula is as follows:

[0067] (2)

[0068] In equation (2), i represents the pixel band, including the original bands such as blue light, green light, and red light, as well as the representative vegetation index band. Let be the predicted value of the temporal variation trajectory of the i-band of a pixel at time t, where T is the number of days in a year. The intercept is... This is the interannual variation coefficient. , This is the coefficient of variation within the year. , This represents the bimodal variation coefficient within the year. , This represents the coefficient of variation of the three peaks within the year.

[0069] For example, Figure 2 This is a schematic diagram of the time-series variation trajectory after harmonic fitting provided in an embodiment of this application. Taking pixel point I, with geographical coordinates of 112.401421° E, 39.047458° N, as an example... Figure 2 The NDVI time series observations and their harmonic fitting time series trajectory for this pixel from 2001 to 2024 are shown. From 2001 to 2024, the NDVI time series trajectory of this pixel is divided into 7 consecutive time periods, during which 6 abrupt events (A, B, C, D, E, and F) caused discontinuities.

[0070] MVC is one of the most widely used annual composite strategies in remote sensing vegetation time series analysis. This algorithm directly selects the maximum NDVI value of the year on a pixel-by-pixel basis to reflect the vegetation growth season status each year. However, when the damage occurs in the same year after the maximum NDVI value of the year, MVC will still select the maximum NDVI value of the year before the damage to represent the annual vegetation status. Although this result objectively reflects the highest growth level of vegetation in the spectral sense, it will create a "time lag" problem in the application scenario of damage time point identification, that is, the damage identification result lags behind the actual damage time.

[0071] For example, Figure 3 This diagram illustrates the time lag problem inherent in the MVC method during damage monitoring. Figure 3 As shown, the pixel was damaged in November 2021, and the NDVI fitting trajectory immediately dropped sharply. However, the MVC method still selected the NDVI peak value of the growing season (about 0.69) as the representative value for 2021, which led to it misjudging the year of damage as 2022.

[0072] To address the "time lag" problem in perturbation timing identification using the MVC method, the following step S124 can be performed in some possible implementations. S124 involves detecting breakpoints in the time-series change trajectory, and based on the breakpoint detection results, automatically determining and extracting annual fitting coefficients. These annual fitting coefficients include the intercept, slope, and three pairs of sine and cosine terms. The specific steps are as follows:

[0073] S1241, Based on the temporal change trajectory fitted by the CCDC algorithm in step S123, identify the breakpoints in the temporal change trajectory of each pixel. The breakpoints are used to indicate the transition of a pixel from one state to another in the time series.

[0074] S1242, extract the occurrence time, start and end times of the segments before and after the breakpoint, and fitting coefficients corresponding to each breakpoint.

[0075] S1243, using historical years as the time unit, for any given year of a pixel: if there are no discontinuities within the year, the annual fitting coefficients are directly taken from the fitting coefficients of the unique segment corresponding to that year. If there are discontinuities within the year, the annual fitting coefficients should be uniformly updated to the fitting coefficients of the last segment. This ensures the uniqueness of the fitting coefficients at the annual scale, which can be used to characterize the latest state of the Earth's surface after a sudden event.

[0076] For example, taking pixel I with geographic coordinates of 112.401421° E, 39.047458° N as an example, the annual fitting coefficient set for this pixel is shown in Table 1. This pixel experienced its first abrupt change in December 2007, so the fitting coefficients for 2007 should be uniformly updated to the fitting coefficients for the time period after the first breakpoint, i.e., the second segment; another abrupt change occurred in December 2009, so the fitting coefficients for 2009 should be uniformly updated to the fitting coefficients for the third segment; in this way, dynamic updating and stable extraction of the fitting coefficients for each year can be achieved.

[0077] Table 1. Annual Fitting Coefficients Set for the Temporal Variation Trajectory of NDVI at Pixel Point I

[0078]

[0079] S125 generates a vegetation index curve covering the entire year's time series based on the temporal change trajectory and annual fitting coefficient. The vegetation index curve is used to characterize the annual vegetation growth phenological characteristics under the latest surface conditions of the current year.

[0080] Specifically, after completing step S124, an annual traversal analysis can be performed on each pixel. The extracted annual fitting coefficients are substituted into the harmonic function expression shown in Equation 2 to generate a continuous representative vegetation index curve covering the entire year's time series.

[0081] For S13, in some possible implementations, the vegetation index curve can be ergonomically analyzed, and its annual maximum value index VI can be extracted respectively. max VI max The corresponding time period represents the peak season for vegetation growth and the annual minimum index VI. min VI min The corresponding time period represents the non-growing season of vegetation.

[0082] In some possible implementations, the annual extreme value index of the representative vegetation index curve can be corrected based on the updated annual fit coefficient. The annual extreme value index includes the annual maximum value index and the annual minimum value index.

[0083] At this point, if the maximum value VI is applied directly... max Trajectory analysis for damaged land identification effectively addresses the issues of temporal instability and damage time lag. However, the trajectory may still be affected by interannual trends, meaning that during periods when no pixel breaks, its annual VI... max It may show a continuous upward or downward trend, because the harmonic function fitted by the CCDC algorithm contains a slope term, and VI max The magnitude of change is related to the slope; larger interannual variation trends can cause significant interference and reduce the sensitivity of damage identification.

[0084] To effectively address the aforementioned issues of interannual trend interference and low identification sensitivity, this application implements step S14.

[0085] S14, based on the annual maximum value indicator VI max and the minimum value index VI within the year min The Intra-Annual Amplitude Vegetation Index (IAAVI) was determined, which is defined as the annual VI. max With VI min The difference is a relative quantity used to quantify the growth fluctuation of surface vegetation in the current year under its latest condition.

[0086] In some possible implementations, the modified annual maximum value index VI can be used as a reference. max and the revised annual minimum index VI min Confirm IAAVI.

[0087] The annual amplitude vegetation index can effectively remove the interannual trend effect introduced by the slope term, thereby improving the sensitivity of identifying damage-reclamation events. The specific calculation formula is as follows:

[0088] (3)

[0089] In equation (3), The vegetation index value represents the annual amplitude of a pixel in year y. This represents the maximum value of the remote sensing index for a pixel within the year y. This represents the minimum remote sensing index of a pixel within year y.

[0090] For example, Figure 4 This is a schematic diagram of the annual amplitude vegetation index (IAAVI) for different time segments provided in the embodiments of this application. Taking pixel point I with geographical coordinates of 112.401421° E, 39.047458° N as an example, Figure 4 The annual amplitude vegetation index (IAAVI) is shown at different time segments. A continuous fitting curve reflecting the latest surface state is constructed using this, and the intra-annual phenological difference is calculated. In the CCDC model, discontinuities are typically used to indicate the transition of pixels from one state to another in a time series. However, due to the combined effects of remote sensing noise, sensor differences, and seasonal disturbances, CCDC may introduce anomalous discontinuities unrelated to actual surface changes during fitting and segmented reconstruction. These anomalies often manifest as a state existing only briefly in the time series, such as lasting only 1-2 years, thus forming a non-realistic structure of multiple consecutive "pseudo-discontinuities." Failure to identify and remove these anomalies will significantly interfere with subsequent damage identification, reclamation monitoring, and type classification.

[0091] To effectively solve the above problems, this application embodiment executes S15 to detect outliers in the interannual IAAVI variation trajectory and identify and correct "pseudo-breakpoints" generated during the CCDC model fitting process.

[0092] S15 includes the following specific steps:

[0093] S151 detects the interannual IAAVI variation trajectory and performs outlier identification and correction for single years.

[0094] The sliding window size is set to 3, and the interannual IAAVI variation trajectory is traversed pixel by pixel (represented by s in the formula). When equation (4) is satisfied, the IAAVI value corresponding to year y is determined to be an outlier. The specific formula is as follows:

[0095] (4)

[0096] In equation (4), These correspond to the three years within the sliding window, respectively; s(y) is the value of the interannual IAAVI variation trajectory in year y. The preset threshold can be reasonably set according to the interannual fluctuation range of vegetation in the study area. For example, this application adopts... NDVI = 0.02; sign() indicates the direction of deviation; abs() indicates the absolute value.

[0097] After the traversal is complete, the year index array of the single-year outliers can be obtained as follows: .right Clustering is performed based on the difference between adjacent elements. When the difference between any two adjacent elements is 1, these consecutive years are grouped into the same single-year outlier group, thus obtaining... .

[0098] For each abnormal group According to equations (5) to (7), the replacement costs on both sides are constructed respectively, and the side with the smaller cost is selected to uniformly correct all outliers in the group. Let... A certain abnormal group (1≤m≤n≤j), then the specific formula is as follows:

[0099] (5)

[0100] (6)

[0101] (7)

[0102] In the formula, L and R represent the single-year abnormal groups, respectively. Non-abnormal neighbor values ​​on the left and right sides , These are the single-year abnormal groups. All outliers are replaced with the total cost of L and R, and the meanings of other variables are the same as in formula (4).

[0103] S152 detects the interannual IAAVI variation trajectory and performs short-term outlier identification and correction over two years.

[0104] The sliding window size is set to 4. The interannual IAAVI variation trajectory, corrected for anomalous conditions in each single year, is traversed pixel-by-pixel in the interannual IAAVI variation trajectory. In the formula... This means that when the condition of equation (8) is met, the IAAVI values ​​corresponding to year y and year y+1 are determined to be outliers. The specific formula is as follows:

[0105] (8)

[0106] In the formula, This represents the value of the interannual IAAVI variation trajectory in year y after correction for single-year anomalous conditions. and The preset threshold can be reasonably set according to the interannual fluctuation range of vegetation in the study area.

[0107] Exemplary, the embodiments of this application employ NDVI = 0.02, NDVI = 0.005; sign() indicates the direction of deviation; abs() indicates the absolute value.

[0108] After the traversal is complete, the year index array of the two outliers can be obtained as follows: .right Clustering is performed based on the differences between adjacent elements. When the difference between any two adjacent elements is 1 or 2, these consecutive years are grouped into the same single-year outlier group, thus obtaining... .

[0109] For each abnormal group According to equations (9) to (11), the replacement costs on both sides are constructed respectively, and the side with the smaller cost is selected to uniformly correct all abnormal points in the group.

[0110] set up A certain abnormal group (1≤m≤n≤j), then the specific formula is as follows:

[0111] (9)

[0112] (10)

[0113] (11)

[0114] In the formula, , They are respectively abnormal groups Non-abnormal neighbor values ​​on the left and right sides , The abnormal groups are respectively Replace all corresponding outliers with and The total cost, and the meanings of other variables are the same as in formula (8).

[0115] For step S16, each pixel in the area to be investigated can be analyzed one by one according to the corrected interannual IAAVI variation trajectory to determine whether it has been damaged and to determine the year in which the damage occurred.

[0116] In areas with stable vegetation cover, the intra-annual phenological variation of pixels is large, corresponding to higher inter-annual IAAVI values. When the surface is disturbed or vegetation is damaged, the intra-annual variation weakens, and the IAAVI value decreases significantly. Based on this characteristic, a damage assessment threshold T is set. IAAVI When the IAAVI value of a pixel in a certain year is lower than the threshold T IAAVI If the pixel is damaged, it is determined that the pixel is in a damaged state in that year; otherwise, it is considered that no damage has occurred.

[0117] Furthermore, the first time IAAVI falls below the threshold T IAAVI The year is determined as the year the pixel was damaged. If a pixel repeatedly shows damage exceeding T during the monitoring period... IAAVI Turn to below T IAAVI The change process indicates that the pixel has experienced multiple damage events; each year in which IAAVI changes from above the threshold to below the threshold is recorded as a new damage year.

[0118] Specifically, in order to objectively determine the damage assessment threshold T IAAVI The Otsu automatic thresholding method can be used, aiming to maximize inter-class variance. It automatically searches for the optimal segmentation point to achieve binary classification of continuous values. In practice, the interannual IAAVI variation trajectory of all pixels in the survey area during the monitoring period is used as input samples to construct a global IAAVI value distribution histogram. Based on the bimodal characteristics of the histogram, the Otsu algorithm determines a threshold globally that optimally distinguishes between "high IAAVI" (undamaged state) and "low IAAVI" (damaged state). The obtained threshold is used as T. IAAVI The value standard is used to determine the surface damage status of each pixel in different years.

[0119] In some possible implementations, the damaged land identification method provided in this application embodiment further includes determining whether reclamation has been carried out after the pixel is damaged, and obtaining reclamation year information.

[0120] Based on the determination results of step S16, the interannual IAAVI variation analysis is continued for the damaged pixels in the area to be investigated to determine whether reclamation has been implemented and to determine the year of reclamation implementation.

[0121] With the restoration of surface vegetation, the phenological characteristics of pixels will be significantly enhanced again, manifested as an increase in the annual phenological difference exceeding the damage assessment threshold T. IAAVI Therefore, when the IAAVI value of a pixel is higher than T in a certain year after the damage, IAAVI When that time comes, the pixel is determined to be subject to reclamation, and that year is recorded as the reclamation year.

[0122] To improve the reliability of the assessment, reclamation assessment should meet the following conditions: ① The reclamation year should be later than the most recent damage year of the pixel; that is, when a pixel has been damaged multiple times, the reclamation year must be later than the year corresponding to the last damage event; ② To avoid misjudgment caused by climate or abnormal observations, it can be required that IAAVI is higher than T for n consecutive years (e.g., 2–3 years). IAAVI ③ If the reclamation occurs in the last year of the monitoring period, it can also be considered reclamation, avoiding misjudgment due to the boundary of the monitoring period; ④ If the pixel's IAAVI does not exceed the threshold T during the monitoring period. IAAVI If so, it is determined to be in a state of "no reclamation has been implemented".

[0123] In some possible implementations, the damaged land identification method provided in this application embodiment further includes S17, determining whether the reclamation of the pixel is sufficient.

[0124] Based on the results of step S16, the annual IAAVI variation analysis is conducted on the damaged and reclaimed pixels in the area to be investigated to further assess their post-reclamation status and determine whether the reclamation is sufficient.

[0125] Reclamation adequacy refers to whether, after reclamation, the surface cover characteristics (especially vegetation growth) of a pixel have recovered to the level before damage, or are comparable to the productivity level of similar features in the surrounding area.

[0126] This step is based on the "Intra-Year Amplitude Vegetation Index (IAAVI)" and its constituent elements, specifically the "Intra-Year Maximum Value (VI)". max The dual-feature comprehensive judgment of the "indicator" is carried out as follows:

[0127] ① When the IAAVI value of a pixel recovers to or exceeds the pre-damage level after reclamation, i.e. If the pixel is fully reclaimed, it indicates that its seasonal vegetation change characteristics have returned to their original level.

[0128] ② When the IAAVI value of the reclaimed pixels does not fully recover to the pre-damage level, i.e. At that time, the maximum value within the year of its constituent elements (VI) should be further utilized. max The "VI" indicator is used for auxiliary judgment. If the reclamation results in VI... max When the value reaches or exceeds the stable level of similar land features before the damage or the surrounding undamaged area, it indicates that the pixel has functionally recovered to a new stable state and can be judged as fully reclaimed; otherwise, it is judged as insufficiently reclaimed.

[0129] The damaged land identification method provided in this application also determines the criteria for classifying damaged land types.

[0130] For example, Figure 5This diagram illustrates the classification of damaged land types using the damaged land identification method provided in this application embodiment. Figure 5 As shown, considering the harmonic fitting trajectory of typical remote sensing indices and the interannual IAAVI trajectory, and combining the damage and reclamation status and the adequacy of reclamation, the land in the survey area is classified into four types: ① Undamaged land; ② Damaged land that has not been reclaimed; ③ Damaged land that has not been fully reclaimed; ④ Damaged land that has been fully reclaimed.

[0131] Finally, the damaged land identification method provided in the embodiments of this application is used to verify the accuracy of damaged land identification.

[0132] Based on the damaged land classification standard described in step S18, and combined with the judgment results of steps S15 to S17, the damaged land in the area to be investigated is identified pixel by pixel, the land damage and reclamation status of each pixel is determined, and a spatial distribution map of damaged land is generated.

[0133] After identifying damaged land, accuracy verification is conducted using ground survey samples or high-resolution imagery data. Specifically, this includes: deploying verification samples within typical plots to ensure representativeness of all land types; constructing a confusion matrix based on the verification samples, and calculating indicators such as overall accuracy, user accuracy, producer accuracy, and Kappa coefficient to quantify identification accuracy; when accuracy indicators fail to meet preset thresholds, re-identification and verification are performed by adjusting classification thresholds or optimizing time-series trajectory fitting parameters until the accuracy meets requirements.

[0134] The method for identifying damaged land provided in this application has application scenarios including but not limited to:

[0135] (1) Damaged land identification and rapid survey. For regional-scale land resource monitoring tasks, this method can be used as an automated and rapid technical means to directly identify and map the distribution of damaged land in areas with strong surface disturbance, such as mining areas, transportation construction areas, and industrial land expansion areas, to assist manual surveys and on-site interpretation, and greatly improve the accuracy and timeliness of identification.

[0136] (2) Damaged land identification and results verification. In areas where damaged land surveys have been conducted but complete historical data or the latest monitoring data are lacking, the method provided in this application embodiment can verify and dynamically update existing damaged land results; by analyzing the interannual phenological change trajectory, the time of damage occurrence, the duration of damage and its ecological impact range can also be inverted, providing technical support for the identification of damaged land, supplementary surveys and database updates.

[0137] (3) Dynamic monitoring of the entire process of damage-reclamation. The method provided in this application embodiment can perform long-term dynamic tracking of damaged areas, continuously obtain information on changes in surface ecological conditions, quantify the vegetation restoration process and ecological function improvement effect after reclamation, and provide a scientific basis for reclamation supervision, performance evaluation and ecological restoration effectiveness assessment.

[0138] Example 1

[0139] Using a mining area as the area to be investigated, identify damaged land. The specific steps are as follows:

[0140] 1. Data Acquisition. The original remote sensing imagery data were all obtained from the U.S. Geological Survey, using the GEE platform to acquire Landsat TM / ETM data for the mining area from 1986 to 2024. + / OLI and Sentinel-2 MSI surface reflectance image data, Landsat image processing level is Collection 2 Level-2, Sentinel-2 image processing level is Level 2A.

[0141] 2. Data Preprocessing. The raw remote sensing image data is preprocessed using the GEE platform, including time and cloud cover filtering, cloud removal, image cropping, and mosaicking. The image filtering period is from January to December, and images with cloud cover greater than 40% are removed. Cloud removal is performed using the SimpleCloudScore algorithm on the GEE platform, and image cropping uses the Quality Mosaic algorithm.

[0142] 3. Select NDVI as the representative remote sensing index and calculate the NDVI value for each image using the GEE platform.

[0143] 4. Based on the CCDC algorithm of the GEE platform, harmonic function fitting is performed on the NDVI time-series observations of each pixel in the mining area. The time location of the breakpoint is identified according to the algorithm output, and representative coefficients are extracted year by year. Based on the interannual representative coefficients, an annual continuous fitting curve reflecting the latest state of the surface is constructed. On this basis, the corrected annual maximum and minimum NDVI indices are calculated, and the annual amplitude vegetation index (IAAVI) is further calculated.

[0144] 5. The damage assessment threshold T is selected using the Otsu automatic thresholding method. IAPD The value is 0.16, which determines information such as the year of damage, whether reclamation was carried out, and whether the reclamation was sufficient. The result of the damage year identification is as follows: Figure 6 The results of the classification of damaged land types are as follows: Figure 7 As shown.

[0145] The above is an introduction to the method for identifying damaged land provided in the embodiments of this application. It is understood that the sequence number of each step in the above embodiments does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application. Furthermore, in some possible implementations, each step in the above embodiments may be selectively executed according to actual circumstances; it may be partially or fully executed, without limitation here. In addition, all or part of any feature of any of the above embodiments can be freely and arbitrarily combined without contradiction; the combined technical solution is also within the scope of this application.

[0146] Next, based on the above, the damaged land identification device provided in the embodiments of this application will be described. For details regarding the concepts, formulas, etc., involved in the following content, please refer to the above text.

[0147] Figure 8 This application provides a damaged land identification device as an embodiment. For example... Figure 8 As shown, the damaged land identification device 80 includes a data acquisition module 81, a curve generation module 82, a curve analysis module 83, an IAAVI determination module 84, and a damage identification module 85.

[0148] In the device, the data acquisition module 81 acquires long-term remote sensing image data based on the spatial range of the area to be investigated and the monitoring period. The long-term remote sensing image data includes blue light, green light, red light, near-infrared, and short-wave infrared bands.

[0149] The curve generation module 82 generates vegetation index curves covering the entire year's time series based on long-term remote sensing image data; the vegetation index curves are used to characterize the annual vegetation growth phenology characteristics under the latest surface conditions of the current year.

[0150] The curve analysis module 83 performs a traversal analysis of the vegetation index curves to determine the maximum and minimum values ​​for the year. The period corresponding to the maximum value represents the peak growth season of vegetation, while the period corresponding to the minimum value represents the non-growth season of vegetation.

[0151] The IAAVI determination module 84 determines the annual amplitude vegetation index based on the annual maximum and annual minimum indicators. The annual amplitude vegetation index is a relative quantity formed by the difference between the annual maximum and annual minimum indicators, which is used to quantify the growth fluctuation of the current year's surface vegetation in its latest state.

[0152] The damage identification module 85 detects the interannual IAAVI change trajectory, identifies and corrects outliers in the interannual IAAVI change trajectory; the interannual IAAVI change trajectory is a trajectory of IAAVI value changes in multiple years; and identifies damaged land based on the corrected interannual IAAVI change trajectory.

[0153] As an example of a software functional unit, the data acquisition module 81 may include code running on a computing instance. The computing instance may include at least one of a physical host (computing device), a virtual machine, or a container. Furthermore, the aforementioned computing instance may be one or more. For example, the data acquisition module 81 may include code running on multiple hosts / virtual machines / containers. It should be noted that the multiple hosts / virtual machines / containers used to run the code may be distributed within the same region or in different regions. Further, the multiple hosts / virtual machines / containers used to run the code may be distributed within the same availability zone (AZ) or in different AZs, each AZ including one or more geographically proximate data centers. Typically, a region may include multiple AZs.

[0154] Similarly, multiple hosts / virtual machines / containers used to run this code can be distributed within the same Virtual Private Cloud (VPC) or across multiple VPCs. Typically, a VPC is set up within a region. Communication between two VPCs within the same region, as well as between VPCs in different regions, requires a communication gateway to be set up within each VPC to enable interconnection between VPCs.

[0155] This application also provides a computing device 90. For example... Figure 9 As shown, the computing device 90 includes a bus 92, a processor 94, a memory 96, and a communication interface 98. The processor 94, the memory 96, and the communication interface 98 communicate with each other via the bus 92. The computing device 90 can be a computing device or a terminal device. It should be understood that this application does not limit the number of processors and memories in the computing device 90.

[0156] Bus 92 can be a Peripheral Component Interconnect (PCI) bus or an Extended Industry Standard Architecture (EISA) bus, etc. Buses can be divided into address buses, data buses, control buses, etc. For ease of representation, Figure 9The bus 94 is represented by a single line, but this does not mean that there is only one bus or one type of bus. The bus 94 may include a path for transmitting information between various components of the computing device 90 (e.g., memory 96, processor 94, communication interface 98).

[0157] Processor 94 may include any one or more processors such as a central processing unit (CPU), a graphics processing unit (GPU), a microprocessor (MP), or a digital signal processor (DSP).

[0158] The memory 96 may include volatile memory, such as random access memory (RAM). The processor 94 may also include non-volatile memory, such as read-only memory (ROM), flash memory, hard disk drive (HDD), or solid state drive (SSD).

[0159] The memory 96 stores executable program code, which the processor 94 executes to implement all or part of the steps of the method in the above embodiments. That is, the memory 96 stores instructions for performing all or part of the steps of the method in the above embodiments.

[0160] The communication interface 98 uses transceiver modules, such as, but not limited to, network interface cards and transceivers, to enable communication between the computing device 90 and other devices or communication networks.

[0161] This application provides a computer-readable storage medium storing a computer program that, when run on a processor, causes the processor to perform the method as described in any one of the first aspects.

[0162] It is understood that the processor in the embodiments of this application may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, transistor logic devices, hardware components, or any combination thereof. A general-purpose processor may be a microprocessor or any conventional processor.

[0163] The method steps in the embodiments of this application can be implemented in hardware or by a processor executing software instructions. The software instructions can consist of corresponding software modules, which can be stored in random access memory (RAM), flash memory, read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), registers, hard disks, portable hard disks, CD-ROMs, or any other form of storage medium known in the art. An exemplary storage medium is coupled to the processor, enabling the processor to read information from and write information to the storage medium. Of course, the storage medium can also be a component of the processor. The processor and the storage medium can reside in an ASIC.

[0164] In the above embodiments, implementation can be achieved, in whole or in part, through software, hardware, firmware, or any combination thereof. When implemented in software, it can be implemented, in whole or in part, as a computer program product. A computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the flow or function according to the embodiments of this application is generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in or transmitted through a computer-readable storage medium. The computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (e.g., coaxial cable, fiber optic, digital subscriber line (DSL)) or wireless (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center that integrates one or more available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., DVD), or a semiconductor medium (e.g., solid-state disk (SSD)).

[0165] It is understood that the various numerical designations used in the embodiments of this application are merely for descriptive convenience and are not intended to limit the scope of the embodiments of this application.

Claims

1. A method for identifying damaged land, characterized in that, The method includes: Long-term remote sensing image data is acquired based on the spatial range and monitoring period of the area to be investigated. The long-term remote sensing image data includes blue light, green light, red light, near infrared, and short-wave infrared bands. A vegetation index curve covering the entire year is generated based on the long-term remote sensing image data; the vegetation index curve is used to characterize the annual vegetation growth phenology characteristics under the latest surface conditions of the current year. The vegetation index curve is subjected to traversal analysis to determine the maximum value index and the minimum value index within the year; the period corresponding to the maximum value index within the year represents the peak growth season of vegetation, and the period corresponding to the minimum value index within the year represents the non-growth season of vegetation. The annual amplitude vegetation index (IAAVI) is determined based on the annual maximum value index and the annual minimum value index. The annual amplitude vegetation index is a relative quantity formed by the difference between the annual maximum value index and the annual minimum value index, and is used to quantify the growth fluctuation range of the surface vegetation in the latest state in the current year. Detect the interannual IAAVI variation trajectory, identify and correct outliers in the interannual IAAVI variation trajectory; the interannual IAAVI variation trajectory includes the variation trajectory of IAAVI values ​​in multiple years; Damaged land was identified based on the revised interannual IAAVI variation trajectory. The generation of vegetation index curves covering the entire year's time series based on the long-term remote sensing image data includes: By selecting a vegetation index type, and performing band calculations on each image based on the long-term remote sensing image data, a vegetation index time series of the selected vegetation index type is obtained. The vegetation index is used to indicate the land cover type and vegetation change characteristics. By fitting the time series of vegetation indices for each pixel, the temporal change trajectory of the corresponding pixel is obtained. Breakpoint detection is performed based on the time-series change trajectory, and the annual fitting coefficient is determined based on the breakpoint detection results; A vegetation index curve covering the entire year's time series is generated based on the time-series change trajectory and annual fitting coefficient. The vegetation index curve is used to characterize the annual vegetation growth phenology characteristics under the latest surface conditions of the current year. The step of detecting breakpoints based on the time-series change trajectory and determining the annual fitting coefficients based on the breakpoint detection results includes: Identify breakpoints in the time-series change trajectory, whereby the breakpoints indicate the transition of a pixel from one state to another in the time series; Extract the occurrence time of the breakpoint, the start and end times of the segments before and after the breakpoint, and the fitting coefficients; Using years as the time unit, if there are no breaks within a year, the annual fitting coefficient is taken from the fitting coefficient of the unique segment corresponding to that year; if there are breaks within a year, the annual fitting coefficient is updated to the fitting coefficient of the last segment; the fitting coefficient is used to characterize the latest state of the surface after the abrupt event at the break point occurs.

2. The method according to claim 1, characterized in that, The step of performing band-wise calculations on the long-term remote sensing image data to obtain vegetation index values ​​includes: Based on the type, stage and surface characteristics of land damage-reclamation in the area to be investigated, vegetation index types are selected, including normalized vegetation index, enhanced vegetation index, comprehensive forest characteristic index and / or normalized degradation fraction index. The long-term remote sensing image data is processed image by image to determine the vegetation index value corresponding to the vegetation index type.

3. The method according to any one of claims 1-2, characterized in that, The vegetation index curve includes breakpoints. The step of performing a traversal analysis on the vegetation index curve to determine the maximum value index and the minimum value index within the year also includes: correcting the maximum and minimum value index within the year of the representative remote sensing index based on the updated annual fitting coefficient. The maximum and minimum value index within the year includes the maximum value index and the minimum value index within the year.

4. The method according to any one of claims 1-2, characterized in that, The method of identifying damaged land based on the annual amplitude vegetation index change trajectory includes: The annual amplitude vegetation index change trajectory is detected to identify "pseudo-breakpoints"; the "pseudo-breakpoints" are outliers that are unrelated to actual surface changes. The outliers are corrected to obtain the corrected interannual IAAVI variation trajectory; Based on the corrected interannual IAAVI variation trajectory, each pixel in the area to be investigated is analyzed to identify damaged land.

5. The method according to any one of claims 1-2, characterized in that, The method of identifying damaged land based on the annual amplitude vegetation index change trajectory includes: Using the IAAVI change trajectory of all pixels in the area to be investigated during the monitoring period as input samples, the Otsu automatic thresholding method is adopted. Based on the principle of maximizing inter-class variance, the damage judgment threshold is determined by automatically searching for the optimal segmentation point. The damage judgment threshold is used to determine the surface damage status of each pixel in different years. Based on the damage assessment threshold and the IAAVI change trajectory, it is determined whether damage has occurred and the year of damage is identified.

6. The method according to any one of claims 1-2, characterized in that, The method further includes: The annual IAAVI variation trajectory analysis was further performed on the damaged pixels in the area to be investigated to determine whether reclamation was implemented and to identify the year of reclamation implementation. If the determination is to implement reclamation, assess the state after reclamation to determine whether the reclamation is sufficient; Based on the land damage and reclamation status of each pixel, a spatial distribution map of damaged land is generated.

7. The method according to any one of claims 1-2, characterized in that, The method further includes: Based on the harmonic fitting trajectory of the remote sensing index and the interannual IAAVI variation trajectory, combined with the damage and reclamation status and the adequacy of reclamation, the land in the area to be investigated is classified into the following types: undamaged land; damaged land that has not been reclaimed; damaged land that has not been fully reclaimed; and damaged land that has been fully reclaimed.

8. A computing device, comprising: At least one memory for storing programs; At least one processor is configured to execute a program stored in the memory, wherein when the program stored in the memory is executed, the processor is configured to perform the method as described in any one of claims 1-2.

Images