A method and system for identifying the risk of critical ecosystem transformation based on multi-source remote sensing vegetation indicators

Abstract

A method and system for identifying the critical transformation risk of an ecosystem based on multi-source remote sensing vegetation indicators are disclosed. The method includes: S10: acquiring and preprocessing vegetation time-series data of the study area; S20: performing signal separation based on multiple time-series decomposition methods; S30: selecting the best separation result from the multiple time-series decomposition methods; S40: calculating statistical indicators for early warning of critical transformation; S50: performing linear trend analysis on the time series of early warning indicators for each pixel, and calculating the slope b of the indicator value changing over time using the least squares method; S60: directly summing the trend values ​​b of each indicator obtained in step S50 to obtain the pixel-scale consistency score CS; S70: constructing a multi-indicator covariance matrix; calculating the standard deviation MS of the Mahalanobis distance time series; S80: constructing a critical transformation risk index R, and judging the critical transformation risk of the ecosystem in the study area based on R.

Description

Technical Field

[0001] This invention relates to the field of environmental monitoring and remote sensing technology, specifically to a method and system for identifying critical transformation risks in ecosystems based on multi-source remote sensing vegetation indicators. Background Technology

[0002] Ecosystems may undergo abrupt degradation (critical transformation) under the influence of climate change and human activities, such as grassland desertification and forest decline. In order to provide early warning before such transformations occur, the academic community has proposed the Critical Slowing Down (CSD) theory, which suggests that when a system approaches a critical point, its resilience decreases, which can be identified through changes in statistical characteristics.

[0003] In existing technologies, the identification of critical turning points mainly includes the following methods: Single time-series statistical index methods, which calculate statistics such as autocorrelation coefficient (AR1), variance, skewness, and kurtosis of remote sensing time series to determine whether a critical deceleration signal exists. Single remote sensing variable analysis methods, where existing studies typically select a specific vegetation variable (such as NDVI or LAI) as a representative of the ecosystem state, assuming that this variable can fully reflect changes in system dynamics. Independent interpretation of indicators methods, where different statistical indicators are often used independently or simply superimposed, lacking a systematic analysis of the differences in the dynamic hierarchy reflected by different indicators. The assumption of unified spatial scale, where most studies conduct analyses at a fixed spatial resolution, rarely considering the differences in the spatial structure expression of different statistical indicators.

[0004] The aforementioned methods have shown some shortcomings in practice. For example, a single statistical indicator is insufficient to comprehensively characterize system dynamics. Different statistical indicators reflect different dynamic information, but existing methods often select only one (such as AR1 or variance), leading to a one-sided characterization of system instability. They neglect the behavioral differences of the same remote sensing variable under different statistical indicators. Existing studies often focus on "whether a critical signal appears," ignoring the spatial and statistical behavioral differences of the same vegetation variable under different statistical indicators, resulting in insufficient understanding of critical mechanisms. There is a lack of a comparability analysis framework between different remote sensing variables. Structural variables (such as LAI), functional variables (such as SIF), and water variables (such as VOD) have significant differences in physical meaning, but existing methods lack a unified framework for systematically comparing their critical statistical behaviors. Spatial clustering is insufficiently characterized. Most studies only conduct pixel-level statistical analysis, failing to quantitatively assess the spatial clustering degree and structural differences of high-risk areas. The interpretation of results relies on experience and lacks systematic quantitative support. The interpretation of the spatial distribution of critical signals largely depends on empirical judgment, lacking a unified quantitative analysis process. Summary of the Invention

[0005] Therefore, to address the shortcomings of existing technologies, this invention proposes a method and system for identifying critical transition risks in ecosystems based on multi-source remote sensing vegetation indicators. This method systematically distinguishes and quantifies the differences in ecosystem dynamics reflected by different critical statistical indicators within the same vegetation remote sensing variable; constructs a collaborative analysis method for multiple statistical indicators to avoid misjudgments caused by a single indicator; quantitatively describes the risk distribution structure and clustering characteristics corresponding to different statistical indicators at a spatial scale; and establishes a unified critical transition analysis framework applicable to multiple vegetation remote sensing variables.

[0006] More specifically, according to one aspect of the present invention, a method for identifying critical transformation risk in ecosystems based on multi-source remote sensing vegetation indicators is provided, comprising:

[0007] S10: Obtain vegetation time series data of the study area and perform preprocessing. The data includes leaf area index (LAI), normalized vegetation index (NDVI), solar induced chlorophyll fluorescence (SIF), and vegetation optical depth (VOD). The preprocessing includes missing value imputation, outlier handling, and time scale unification.

[0008] S20: Signal separation is performed based on multiple time series decomposition methods to obtain the residual time series (Ri) after removing long-term trends and seasonal cycles. t The various time series decomposition methods include the seasonal-trend decomposition method (STL), the structural time series analysis method (STA), and the state-space model method (SSM).

[0009] S30: Based on the residual time series from step S20, select one of the various time series decomposition methods as the subsequent analysis method;

[0010] S40: Calculation of early warning statistical indicators for critical transformation, including calculating early warning statistical indicators for the ecosystem using a sliding window method on the residual time series of the method selected in step S30, to obtain the time-varying indicator series Z. T The indicators include first-order autocorrelation coefficient (AR(1)), variance (Var), skewness, and kurtosis.

[0011] S50: Perform linear trend analysis on the time series of early warning indicators for each pixel, and use the least squares method to calculate the slope b of the indicator value changing over time, which is used to characterize the direction and intensity of changes in ecosystem stability, including:

[0012] Based on the time series of early warning indicators for each pixel {Z1, Z2, ..., Z} T}, that is, Z T Construct a linear regression model Where Zt is the index value, t is the time sequence number, a is the intercept, b is the slope, and εt is the residual term;

[0013] The slope b is calculated using the least squares method, where b>0 indicates that the index increases with time and the risk of decreased stability increases; b<0 indicates that the index decreases with time and stability recovers; b=0 indicates that there is no significant change in the long term.

[0014] S60: The trend values ​​b of each indicator obtained in step S50 are directly summed to obtain the pixel-scale consistency score CS.

[0015] S70: The initial ecosystem state of the study area is used as the baseline state, that is, the baseline state is the first data point within the same time window, based on the indicator sequence Z. T Construct a multi-index covariance matrix; calculate the Mahalanobis distance of each pixel in the target period relative to the baseline state in the multi-index feature space, and further calculate the standard deviation MS of the Mahalanobis distance time series;

[0016] S80: Construct a critical transformation risk index R, and determine the critical transformation risk of the ecosystem in the study area based on R.

[0017]

[0018] Where R(x,y) is the risk index of pixel (x,y), CS(x,y) is the multi-indicator consistency score of the pixel, and MS(x,y) is the standard deviation of Mahalanobis distance of the pixel; w1 and w2 are the fusion weights, min(CS) and max(CS) are the minimum and maximum values ​​of the multi-indicator consistency score in the study area, respectively, and min(MS) and max(MS) are the minimum and maximum values ​​of the standard deviation of Mahalanobis distance in the study area, respectively.

[0019] According to an embodiment of the present invention, in step S10, missing values ​​are filled using linear interpolation and moving average methods; outlier handling includes identifying outliers using the Z-score method or box plot method, and performing interpolation or neighborhood averaging substitution.

[0020] According to an embodiment of the present invention, in step S30, the reconstruction performance of different time series decomposition methods is quantitatively compared by calculating the root mean square error (RMSE) of the residual time series; RMSE is used as an evaluation index to quantitatively compare the stability of the residual signals of different methods, and then a method is selected based on the comparison results.

[0021] According to an embodiment of the present invention, in step S40, the sliding window length is set to 12 months, the window step size is 12 months, and the residual time series R of the method selected in step S3 is... tA sliding window with a length of 12 months and a step size of 12 months is used; the samples within a certain window are denoted as: {r1, r2, ..., r t ,…,r n}, n=12

[0022] Calculate each early warning indicator for each window separately, using the following formula:

[0023] The formula for calculating AR(1) is as follows: ,in This represents the mean value within the window.

[0024] The formula for calculating variance (Var) is as follows:

[0025] Skewness is calculated as follows: , where s is the standard deviation of each value within the window.

[0026] The formula for calculating kurtosis is as follows: .

[0027] According to an embodiment of the present invention, kurtosis employs excessive kurtosis. excess =Kurtosis-3.

[0028] According to an embodiment of the present invention, S50 further includes using a statistical significance test on the trend, and performing weight reduction processing on the trend results that fail the significance test, so as to reduce the impact of random noise on the subsequent risk assessment results.

[0029] According to an embodiment of the present invention, the weight reduction process includes reducing the slope b by weight, as shown in the following formula:

[0030] Where w is the weight, which varies with significance p.

[0031] .

[0032] According to an embodiment of the present invention, S40 further includes plotting the spatial distribution of risk based on the critical transformation risk index.

[0033] According to another aspect of the present invention, an ecosystem critical transformation risk identification system based on multi-source remote sensing vegetation indicators is also provided, comprising:

[0034] The vegetation time series data acquisition and preprocessing module is used to acquire and preprocess vegetation time series data of the study area. The data includes leaf area index (LAI), normalized vegetation index (NDVI), solar induced chlorophyll fluorescence (SIF), and vegetation optical depth (VOD). The preprocessing includes missing value imputation, outlier handling, and time scale unification.

[0035] The signal separation module is used to separate signals based on various time series decomposition methods to obtain residual time series (Ri) after removing long-term trends and seasonal cycles. t The various time series decomposition methods include the seasonal-trend decomposition method (STL), the structural time series analysis method (STA), and the state-space model method (SSM).

[0036] The method comparison and optimization module, based on the various time series decomposition methods and their residual time series described in the signal separation module, selects one method and its residual time series for subsequent analysis;

[0037] The critical transition early warning statistical index calculation module is used to calculate the early warning statistical index for critical transitions. This includes using a sliding window method to calculate the ecosystem's early warning statistical index based on the residual time series of the methods selected by the method comparison and optimization module, resulting in the time-varying warning index sequence Z. T The indicators include first-order autocorrelation coefficient (AR(1)), variance (Var), skewness, and kurtosis.

[0038] The linear trend analysis module is used to perform linear trend analysis on the time series of early warning indicators for each pixel. It uses the least squares method to calculate the slope b of the indicator value over time, which is used to characterize the direction and intensity of changes in ecosystem stability. b>0 indicates that the indicator increases over time, and the risk of stability decline increases; b<0 indicates that the indicator decreases over time, and stability recovers; b=0 indicates no significant change in the long term.

[0039] The Consistency Score (CS) calculation module is used to directly sum the trend values ​​b of each indicator obtained from the linear trend analysis module to obtain the pixel-scale consistency score (CS).

[0040] The module for calculating the standard deviation (MS) of the Mahalanobis distance time series is used to construct a multi-indicator covariance matrix based on the initial ecosystem state of the study area as the baseline state and the time-varying index sequence (ZT). It calculates the Mahalanobis distance of each pixel in the target period relative to the baseline state in the multi-indicator feature space and further calculates the standard deviation (MS) of the Mahalanobis distance time series.

[0041] The critical transformation risk index R module is used to construct the critical transformation risk index R and to determine the critical transformation risk of the ecosystem in the study area based on R.

[0042]

[0043] Where R(x,y) is the risk index of pixel (x,y), CS(x,y) is the multi-indicator consistency score of the pixel, and MS(x,y) is the standard deviation of Mahalanobis distance of the pixel; w1 and w2 are the fusion weights, min(CS) and max(CS) are the minimum and maximum values ​​of the multi-indicator consistency score in the study area, respectively, and min(MS) and max(MS) are the minimum and maximum values ​​of the standard deviation of Mahalanobis distance in the study area, respectively.

[0044] Compared with existing technologies, the present invention has the following advantages: it improves the reliability of identifying critical turning points in ecosystems; it avoids misjudgments caused by single statistical indicators; it reveals the dynamic differences of the same remote sensing variable under different statistical indicators; it enhances the objectivity and quantification of the spatial distribution interpretation of critical risks; and it is applicable to multi-regional and multi-scale ecological monitoring applications. Attached Figure Description

[0045] Figure 1 This is a flowchart illustrating the method for identifying critical transformation risks in ecosystems based on multi-source remote sensing vegetation indicators, according to the implementation plan.

[0046] Figure 2 A schematic diagram of an ecosystem critical transformation risk identification system based on multi-source remote sensing vegetation indicators, according to the implementation plan;

[0047] Figure 3 This is a schematic diagram of the electronic device structure according to an embodiment of the present invention;

[0048] Figure 4 The image shows the RMSE detection results of the residuals obtained by various time series decomposition methods in the ecosystem critical transformation risk identification method based on multi-source remote sensing vegetation indicators according to the implementation plan.

[0049] Figure 5 A consistency score map of NDVI in the ecosystem critical transformation risk identification method based on multi-source remote sensing vegetation indicators according to the implementation plan;

[0050] Figure 6 Mahalanobis distance data for NDVI in the ecosystem critical transformation risk identification method based on multi-source remote sensing vegetation indicators according to the implementation plan;

[0051] Figure 7 The image shows the critical transformation risk identification results of LAI based on the multi-source remote sensing vegetation index-based ecological system critical transformation risk identification method according to the implementation plan.

[0052] Figure 8 A graph showing the verification results of drought frequency based on the implementation plan;

[0053] Figure 9 This is the result of long-term moisture verification based on the implementation plan. Detailed Implementation

[0054] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0055] It should be understood that the seasonal-trend decomposition method (STL), structural time series analysis method (STA), state-space model method (SSM), AR1, variance, skewness and kurtosis, Kalman filtering, Mahalanobis distance stability index and other methods and tools involved in this invention are known. Therefore, this invention focuses on how to combine these model tools to realize the ecosystem critical transformation risk identification method based on multi-source remote sensing vegetation indicators.

[0056] Figure 1 This is a flowchart illustrating the ecosystem critical transformation risk identification method based on multi-source remote sensing vegetation indicators, according to the implementation plan. Please refer to [link / reference]. Figure 1 The implementation plan's method for identifying critical transformation risks in ecosystems based on multi-source remote sensing vegetation indicators may include the following steps.

[0057] S1: Acquisition and Preprocessing of Vegetation Time Series Data

[0058] Monthly time-series data of various vegetation remote sensing indicators within the study area were acquired, including but not limited to leaf area index (LAI), normalized difference vegetation index (NDVI), solar-induced chlorophyll fluorescence (SIF), and vegetation optical depth (VOD). Data sources are shown in Table 1 below.

[0059] index Data Product Examples Spatial resolution Time range Example of acquisition method Leaf Area Index (LAI) MODIS LAI(MOD15A2H / MCD15A3H) 500 m–1 km 2000 to present NASA Earthdata or LAADS DAAC Normalized Difference Vegetation Index (NDVI) MODIS NDVI (MOD13A2 / MYD13A2) or AVHRR NDVI 1 km 1980s to present (depending on the product) NASA LP DAAC, NOAACLASS Sun-induced chlorophyll fluorescence (SIF) GOSIF, TROPOMI SIF, OCO-2 SIF derivative products 0.05° or higher Since the 2000s Harvard GOSIF Database ESA / TROPOMI Data Center Vegetation optical depth (VOD) SMAP VOD, AMSR-E / 2 VOD (LPRM) 0.25° 2002 to present NASA NSIDC, VODCA Global Convergence Products

[0060] Then, the time series data is processed by imputing missing values, handling outliers, and standardizing the time scale to ensure the integrity and comparability of the time series.

[0061] More specifically, when observation gaps or missing data exist, the following methods can be used: time series interpolation: linear interpolation, spline interpolation; moving average method (3-month or 5-month window smoothing), etc. In this example, linear interpolation and moving average are used to fill relatively small gaps, and interpolation based on seasonal patterns is used for large gaps. To avoid anomalous jumps caused by sensor errors or aerosol contamination, outliers (such as those exceeding ±3σ) can be identified using the Z-score method or box plot method, and then interpolation or neighborhood averaging can be used to replace these outliers. Since different sensor products have different time steps and projection formats, standardization is required: including time scale standardization (all samples are resampled to a monthly scale using maximum or mean synthesis) and spatial matching (if necessary), reprojecting to a unified coordinate system and resolution.

[0062] Finally, time series data for each remote sensing indicator (e.g., a pixel) were obtained, including time series data for leaf area index (LAI), normalized difference vegetation index (NDVI), solar induced chlorophyll fluorescence (SIF), and vegetation optical depth (VOD).

[0063] S2: Signal separation based on multiple time series decomposition methods

[0064] To separate long-term trends, seasonal periodic components, and short-term disturbance signals from vegetation time series, the vegetation time series of each pixel is processed using several time series decomposition methods, namely, the seasonal-trend decomposition method (STL), the structural time series analysis method (STA), and the state-space model method (SSM). These methods are known in themselves, and therefore this invention will only provide a brief description.

[0065] Seasonal-Trend Decomposition Method (STL): Based on Locally Weighted Regression (LOESS), this method performs nonparametric decomposition of time series, decomposing the original series into a trend term, a seasonal term, and a residual term. The seasonal period is set to 12 months, and a robust estimation model is used to reduce the impact of extreme values ​​on the decomposition results.

[0066] Structural Time Series Analysis (STA) method: Construct a structural time series model containing local level terms and annual periodic seasonal terms, and estimate the model parameters through state-space representation and Kalman filtering and maximum likelihood estimation to obtain the corresponding residual terms.

[0067] State-space model (SSM) method: By constructing a recursive relationship between system state and observation, it uses extended Kalman filter to perform dynamic state estimation of time series and extracts the deviation between observation and predicted state as residual term.

[0068] The above three decomposition methods all provide three types of information separated from each time series, as shown in Table 2 below:

[0069] Element significance Trend Item Represents long-term vegetation changes (such as climate-driven or ecological restoration). Seasonal items The inherent periodic changes in vegetation (temperature, light, and water regulation) residuals Short-term disturbance signals reflect rapid responses to drought, pests and diseases, fires, and human intervention.

[0070] Ultimately, the implementation example uses the residual time series after removing long-term trends and seasonal cycles as an indicator of the rapid dynamic response of the ecosystem to external disturbances.

[0071] S3: Optimal Selection of Time Series Decomposition Methods

[0072] The reconstruction performance of different time series decomposition methods was quantitatively compared by calculating the root mean square error (RMSE) of the residual time series. The time series decomposition method with the lowest overall reconstruction error and highest residual stability among multiple vegetation indices was selected as the preferred method for subsequent analysis.

[0073] The implementation plan aims to quantitatively compare the stability of residual signals obtained by three time series decomposition methods: STL, STA, and SSM. The evaluation metric used is RMSE (root mean square error).

[0074] RMSE (Root Mean Square Error) is a statistical indicator that measures the magnitude of error and is commonly used to evaluate prediction errors or the intensity of abnormal fluctuations.

[0075] For a certain pixel's residual sequence r1, r2...r n Its RMSE is defined as: RMSE = sqrt((r1) 2 +r2 2 +…+r n 2 ) / n)

[0076] The larger the residual or the more drastic the fluctuation, the higher the RMSE, indicating that the short-term disturbance separated by this method contains large irregular fluctuations; the lower the RMSE, the more stable the residual signal and the less noise.

[0077] In this study scenario, the residuals are the remaining parts after removing the trend and seasonal terms from the time series decomposition model. Therefore, RMSE can reflect: the model's explanatory power for seasonality and trend, short-term changes in vegetation status or noise intensity, and the difference between the results of different models.

[0078] Each remote sensing index (such as Lai, ndvi, etc.) is processed using three time series decomposition methods (STL, STA, SSM) to obtain residual data corresponding to the three methods. Then, the optimal data for each remote sensing index is obtained based on the aforementioned RMSE detection. For example, in this embodiment, the residual sequence output by the STL decomposition method is preferably used. Figure 4 This is a graph showing the RMSE detection results of residuals obtained by various time series decomposition methods in a study area according to an embodiment of the present invention.

[0079] S4: Calculation of statistical indicators for early warning of critical transformation

[0080] Based on the critical deceleration theory, the residual time series obtained in step S3 is used to calculate the early warning statistical indicators of the ecosystem using a sliding window method. The early warning statistical indicators include: first-order autocorrelation coefficient (AR(1)); variance (Var); skewness; and kurtosis. The sliding window length can be set to 12 months, and the window step size can be 12 months, to obtain the evolution sequence Z of each statistical indicator over time. T .

[0081] The residual sequence r obtained in step S3 t (For example, on a monthly scale), a sliding window with a length of 12 months and a step size of 12 months is used.

[0082] Let the samples within a certain window be denoted as: {r1, r2, ..., r n}, n=12

[0083] Calculate the above four early warning statistical indicators for each window to obtain the early warning indicator sequence that changes over time.

[0084] The calculation formulas for each early warning statistical indicator are as follows:

[0085] The first-order autocorrelation coefficient (AR(1)) reflects the system's recovery speed to external disturbances; a larger value indicates a slower recovery. Its calculation formula is as follows:

[0086] ,in This represents the mean value within the window.

[0087] Variance (Var) measures the magnitude of system state fluctuations, and its calculation formula is as follows:

[0088] ;

[0089] Skewness reflects the degree to which a distribution is skewed relative to its mean, and it is calculated as follows:

[0090] , where s is the standard deviation of each value within the window;

[0091] Kurtosis measures the flatness or narrowness of a probability distribution, and its formula is as follows:

[0092] ;

[0093] Kurtosis can also be used. excess The calculation formula is as follows:

[0094] Kurtosis excess =Kurtosis-3.

[0095] The specific calculation process for the above formulas can be implemented using Python code.

[0096] S5: Early Warning Indicator Trend Analysis

[0097] Linear trend analysis was performed on the time series of early warning indicators for each pixel, and the slope of the indicator value changing over time was calculated using the least squares method to characterize the direction and intensity of changes in ecosystem stability.

[0098] More specifically, based on the time series of each early warning indicator {Z1, Z2, ..., Z} for each pixel T}, that is, Z T Construct a linear regression model Where Zt is the index value, t is the time sequence number (window number), a is the intercept, b is the slope, and εt is the residual term;

[0099] Then the slope b can be calculated using the least squares method:

[0100] ,in ;

[0101] The obtained slope b represents the changing trend of early warning indicators:

[0102] Where b>0 indicates that the risk of decreased stability increases with time; for example, for typical critical deceleration indicators such as AR1 and variance, an increase means a slower system recovery and increased volatility. b<0 indicates that the indicator decreases with time, and stability recovers. b=0 indicates no significant change in the long term. The larger the |b| value, the stronger the trend and the more significant the change. Therefore, the slope b can be used as a scaled quantitative indicator to measure the system's approach to the critical inflection point.

[0103] To avoid the influence of random fluctuations on the judgment, statistical significance tests can be used on the trend. For example, linear regression t-test or Mann-Kendall trend test (non-parametric, more commonly used in ecological remote sensing) can be used. Let the test p value be p. If p < 0.05: the trend is significant. If p ≥ 0.05: the trend may be caused by random noise.

[0104] Trends that fail the significance test can be weighted to reduce the impact of random noise on subsequent risk assessment results. Specifically, to reduce the interference of insignificant trends in subsequent risk assessment, this invention proposes weighting the slope of the trend. The slope b of the weighted trend is... * The calculation formula is as follows:

[0105] Where w is the weight, which varies with significance p: .

[0106] When p = 0.05, w is approximately 1, indicating no weakening; when p = 0.10, w is 0.5; and when p ≥ 0.5, the larger the value, the closer it is to 0.

[0107] S6: Multi-indicator Consistency Score Calculation

[0108] To integrate spatial information from different early warning indicators, a multi-indicator consistency score (ConsistencyScore) is constructed. After weight trend weakening in time s5, the trend values ​​of the four types of indicators are directly summed to obtain the pixel-scale consistency score. The higher the consistency score, the greater the likelihood that the trends of multiple indicators in that region simultaneously point to a decline in system stability.

[0109] More specifically, after the above processing, for each pixel, the trend slope b of the four types of early warning indicators has been obtained. k AR1 trend, variance trend, skewness trend, kurtosis trend. Then sum the trends of each indicator:

[0110] ,

[0111] Where b k (x,y) represents the trend slope of the k-th indicator at pixel (x,y).

[0112] The consistency score is used to assess whether the trend signals of different indicators support each other: when all trends are positive and the amplitude is large, it indicates that the system state is changing significantly in the direction of potential instability; when the trend direction is scattered or the amplitude is small, it is considered that the system degradation trend is not obvious; when the overall trend is negative, it can be regarded as the restoration of ecological stability.

[0113] By integrating the four indicators, this indicator can strengthen common trends and reduce the noise impact of individual indicators.

[0114] Figure 5 This is a consistency score map of NDVI in an ecosystem critical transformation risk identification method based on multi-source remote sensing vegetation indicators according to an embodiment of the present invention.

[0115] S7: Multivariate State Deviation Analysis Based on Mahalanobis Distance

[0116] Using the initial ecosystem state of the study area as the baseline state (i.e., the baseline state being the first period's data within the same time window), a multi-indicator covariance matrix is ​​constructed based on the early warning indicator sequence ZT. The Mahalanobis distance of each pixel in the target period relative to the baseline state in the multi-indicator feature space is calculated, and the standard deviation (MS) of the Mahalanobis distance time series is further calculated to quantify the overall deviation of the ecosystem state in the multi-dimensional indicator space. It should be understood that Mahalanobis distance and the covariance matrix are known in this field and will only be briefly explained below.

[0117] Construct a multivariate matrix Xt that changes over time based on the values ​​of each indicator:

[0118]

[0119] Calculate the covariance matrix based on the observed data:

[0120] At the same time, a small perturbation term is added to ensure numerical stability:

[0121] Then, using the first point in time as the reference state:

[0122] Calculate the Mahalanobis distance for each time point t:

[0123] Then, the standard deviation of the Mahalanobis distance time series is calculated:

[0124] A higher standard deviation indicates a larger shift in the ecological state across multiple indicators at different time points, reflecting reduced system stability; a lower standard deviation indicates weaker overall fluctuations in the multivariate signal, suggesting a relatively stable system state. This indicator can utilize multiple indicators simultaneously and consider the correlation structure between them, making it more robust than a single indicator.

[0125] Figure 6 This is a Mahalanobis distance map of NDVI in an ecosystem critical transformation risk identification method based on multi-source remote sensing vegetation indicators according to an embodiment of the present invention.

[0126] S8: Construction and Spatial Representation of the Critical Transition Risk Index

[0127] By integrating the consistency score of multiple indicators with the standard deviation of Mahalanobis distance, a pixel-scale ecological critical transformation risk index is constructed. Areas with higher risk values ​​indicate that the ecosystem is more likely to be approaching critical transformation or experiencing a significant decline in stability.

[0128] More specifically, on a per-pixel basis, two types of indicators are used as the fusion basis: multi-indicator consistency score (CS) and Mahalanobis distance standard deviation (MS). CS characterizes whether the trends of different statistical indicators spatially point to a decline in ecosystem stability; a higher CS indicates stronger trend consistency, which implies greater risk. MS reflects the overall fluctuation range of multiple indicators over time; a higher MS indicates a greater deviation of the system state from the reference state, which implies higher risk.

[0129] The critical transformation risk index R is constructed as follows:

[0130]

[0131] Where R(x,y) is the risk index of pixel (x,y), CS(x,y) is the multi-indicator consistency score of the pixel, and MS(x,y) is the standard deviation of Mahalanobis distance of the pixel; w1 and w2 are fusion weights, w1+w2=1, which can be set according to the importance of the indicators determined by experience, for example, w1 =w2=0.5; min(CS) and max(CS) are the minimum and maximum values ​​of the multi-indicator consistency score in the study area, respectively; min(MS) and max(MS) are the minimum and maximum values ​​of the standard deviation of Mahalanobis distance in the study area, respectively.

[0132] The resulting risk index R(x,y) takes values ​​in the range [0,1]. Pixels with high risk values ​​(R→1) indicate that the ecosystem may be approaching a critical transition point or that its stability is significantly declining; medium- and low-risk values ​​(R≈0) indicate that the system is relatively stable or in a state of recovery. This index can be used for spatial distribution mapping, risk hotspot identification, regional ecological management and intervention decisions, and trend analysis of subsequent time series monitoring, etc.

[0133] Figure 2 This is a schematic diagram of an ecosystem critical transformation risk identification system based on multi-source remote sensing vegetation indicators, according to the implementation plan. Figure 2As shown, the system includes: a vegetation time-series data acquisition and preprocessing module 210, used to acquire and preprocess vegetation time-series data of the study area, the data including leaf area index (LAI), normalized difference vegetation index (NDVI), solar-induced chlorophyll fluorescence (SIF), and vegetation optical depth (VOD); the preprocessing includes missing value imputation, outlier handling, and uniform time scale processing; and a signal separation module 220, used to perform signal separation based on multiple time-series decomposition methods to obtain the residual time series (Rt) after removing long-term trends and seasonal cycles, the multiple time-series decomposition methods including seasonal... The system includes three time series decomposition methods: Situational Time Series (STL), Structural Time Series Analysis (STA), and State-Space Model (SSM). A method comparison and selection module 230, based on the residual time series from the signal separation module, selects one method and its residual time series from the various time series decomposition methods for subsequent analysis. A critical transition early warning statistical index calculation module 240 calculates critical transition early warning statistical indicators, including using a sliding window method to calculate ecosystem early warning statistical indicators from the residual time series of the methods selected by the method comparison and selection module, obtaining a time-varying index sequence ZT. This includes the first-order autocorrelation coefficient (AR(1)), variance (Var), skewness, and kurtosis; the linear trend analysis module 250 is used to perform linear trend analysis on the time series of early warning indicators for each pixel, and uses the least squares method to calculate the slope b of the indicator value over time, which is used to characterize the direction and intensity of changes in ecosystem stability. b>0 indicates that the indicator increases over time, and the risk of stability decline increases; b<0 indicates that the indicator decreases over time, and stability recovers; b=0 indicates that there is no significant change in the long term; the consistency score CS calculation module 260 is used to integrate the linear trend analysis module The trend values ​​b of each indicator are directly summed to obtain the pixel-scale consistency score CS; the standard deviation MS of the Mahalanobis distance time series is calculated by module 270, which uses the initial ecosystem state of the study area as the baseline state, that is, the baseline state is the first (year) data in the same time window, to construct a multi-indicator covariance matrix; the Mahalanobis distance of each pixel in the target period relative to the baseline state in the multi-indicator feature space is calculated, and the standard deviation MS of the Mahalanobis distance time series is further calculated; the critical transformation risk index R module 280 is used to construct the critical transformation risk index R, and to judge the critical transformation risk of the ecosystem in the study area based on R.

[0134] Figure 3 This is a schematic diagram of the electronic device structure according to an embodiment of the present invention, such as... Figure 3 As shown, the electronic device includes a processor 310, a memory 320, an input device 330, and an output device 340; the number of processors 310 in the electronic device can be one or more. Figure 3 Taking a processor 310 as an example; the processor 310, memory 320, input device 330, and output device 340 in the electronic device can be connected via a bus or other means. Figure 3 Take a bus connection as an example.

[0135] The memory 320, as a computer-readable storage medium, can be used to store software programs, computer-executable programs, and modules, such as the program instructions / modules corresponding to the ecosystem critical transformation risk identification method based on multi-source remote sensing vegetation indicators in this embodiment of the invention (e.g., vegetation time series data acquisition and preprocessing module 210; signal separation module 220; method comparison and optimization module 230; critical transformation early warning statistical indicator calculation module 240; linear trend analysis module 250; consistency score CS calculation module 260; Mahalanobis distance time series standard deviation MS calculation module 270; and critical transformation risk index R module 280). The processor 310 executes various functional applications and data processing of the electronic device by running the software programs, instructions, and modules stored in the memory 320, thereby realizing the aforementioned ecosystem critical transformation risk identification method based on multi-source remote sensing vegetation indicators.

[0136] The following section uses an ecosystem in Northeast China as an example to further verify and illustrate the invention through a verification scheme that studies the correlation between critical risk responses of different vegetation indicators and external environmental stresses. Specifically, the effectiveness of the method is verified by constructing a correspondence between climate moisture gradients and critical risk responses.

[0137] This validation uses the pixel as the basic unit of analysis. Through binning statistics, response curve construction, and cross-index comparison, it systematically analyzes the correlation between the critical risk response of different vegetation indicators and climate stress, specifically including:

[0138] For the two climate variables, multi-year average precipitation (Pmean) and multi-year frequency of meteorological drought (Drought-freq), equidistant binning was performed according to the numerical intervals to ensure that each bin contains a sufficient number of pixels to guarantee statistical reliability. At the same time, the distribution characteristics of the variable values ​​were taken into account, so that the bin intervals can not only cover the climate gradient of the entire study area, but also accurately capture the key climate threshold intervals.

[0139] For each climate variable bin, the Critical Deceleration Index (CTI) of the corresponding pixel is extracted. The median statistical method is used to calculate the typical value of the CTI in each bin, eliminating the interference of extreme values ​​on the statistical results. This allows us to obtain representative characteristics of the critical risk response of vegetation systems under different climate background intervals and establish the correspondence between climate gradient and critical risk response.

[0140] Based on the binning statistical results, two-dimensional relationship curves between CTI and Pmean, and between CTI and Drought-freq, were plotted to visually present the changing patterns of critical deceleration indices under different climate stress gradients. Simultaneously, for four remote sensing indices—leaf area index (LAI), normalized difference vegetation index (NDVI), solar-induced chlorophyll fluorescence (SIF), and vegetation optical depth (VOD)—the aforementioned relationship curves were constructed. Through cross-indicator comparisons of curve shape, trends, and threshold ranges, the sensitivity differences of different vegetation state variables to external water stress were systematically assessed, revealing the risk accumulation patterns of various vegetation indices.

[0141] This validation scheme, by constructing a quantitative correlation model between climate stress and critical deceleration characteristics, deepens the understanding of vegetation critical risk response from spatial pattern representation to causal mechanism analysis. The validation results are available in [link to validation report]. Figure 8 and Figure 9 The results show that the risk index systematically increases with increasing drought frequency and decreases with increasing multi-year average precipitation, indicating a causal consistency between the risk identification results and external environmental stress. This consistency suggests that the critical risk signal extracted in this invention has an ecophysical basis, rather than being a statistically random fluctuation. Furthermore, comparative analysis of different vegetation indicators reveals the differences in the sensitivity of various state variables to water stress, further demonstrating the mechanistic explanatory power and application reliability of the method presented in this invention.

[0142] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for identifying critical transformation risk in ecosystems based on multi-source remote sensing vegetation indicators, characterized in that, include: S10: Obtain vegetation time series data of the study area and perform preprocessing. The data includes leaf area index (LAI), normalized vegetation index (NDVI), solar induced chlorophyll fluorescence (SIF), and vegetation optical depth (VOD). The preprocessing includes missing value imputation, outlier handling, and time scale unification. S20: Signal separation is performed based on multiple time series decomposition methods to obtain the residual time series R after removing long-term trends and seasonal cycles. t The various time series decomposition methods include the seasonal-trend decomposition method STL, the structural time series analysis method STA, and the state-space model method SSM. S30: Based on the various time series decomposition methods and their residual time series described in step S20, select one method and its residual time series for subsequent analysis; S40: Calculation of early warning statistical indicators for critical transformation, including calculating early warning statistical indicators for the ecosystem using a sliding window method on the residual time series of the method selected in step S30, to obtain the time-varying warning indicator sequence Z. T The indicators include first-order autocorrelation coefficient AR(1), variance Var, skewness Skewness, and kurtosis. S50: Perform linear trend analysis on the time series of early warning indicators for each pixel, and use the least squares method to calculate the slope b of the indicator value changing over time, which is used to characterize the direction and intensity of changes in ecosystem stability, including: Based on the time series of early warning indicators for each pixel {z1, z2, ..., z t }, that is, Z T Construct a linear regression model , where z t Here, t represents the index value, a represents the time sequence, b represents the intercept, and ε represents the slope. t For residual terms; The slope b is calculated using the least squares method, where b>0 indicates that the index increases with time and the risk of decreased stability increases; b<0 indicates that the index decreases with time and stability recovers; b=0 indicates that there is no significant change in the long term. The trend was tested for statistical significance, and the trend results that failed the significance test were weighted down to reduce the impact of random noise on the subsequent risk assessment results. Weight reduction processing includes reducing the weight of slope b, as shown in the following formula: Where w is the weight, which varies with significance p: ； S60: Directly sum the trend values ​​b* of each indicator obtained in step S50 to obtain the pixel-scale consistency score CS. S70: Using the initial ecosystem state of the study area as the baseline state, based on the early warning indicator sequence Z in step S40. T Construct a multi-index covariance matrix; calculate the Mahalanobis distance of each pixel in the target period relative to the baseline state in the multi-index feature space, and further calculate the standard deviation MS of the Mahalanobis distance time series; S80: Construct a critical transformation risk index R, and judge the critical transformation risk of the ecosystem in the study area based on R; ； Where R(x,y) is the risk index of pixel (x,y), CS(x,y) is the multi-indicator consistency score of the pixel, and MS(x,y) is the standard deviation of Mahalanobis distance of the pixel; w1 and w2 are the fusion weights, min(CS) and max(CS) are the minimum and maximum values ​​of the multi-indicator consistency score in the study area, respectively, and min(MS) and max(MS) are the minimum and maximum values ​​of the standard deviation of Mahalanobis distance in the study area, respectively.

2. The method for identifying critical transformation risk of ecosystems based on multi-source remote sensing vegetation indicators according to claim 1, characterized in that: In step S10, missing values ​​are filled using linear interpolation and moving average methods; outlier handling includes identifying outliers using Z-score or box plot methods and replacing them with interpolation or neighborhood average.

3. The method for identifying critical transformation risk of ecosystems based on multi-source remote sensing vegetation indicators according to claim 1, characterized in that: In step S30, the reconstruction performance of different time series decomposition methods is quantitatively compared by calculating the root mean square error (RMSE) of the residual time series; RMSE is used as an evaluation index to quantitatively compare the stability of the residual signals of different methods, and then the method is selected based on the comparison results.

4. The method for identifying critical transformation risk of ecosystems based on multi-source remote sensing vegetation indicators according to claim 1, characterized in that: In step S40, the sliding window length is set to 12 months, and the window step size is 12 months. The residual time series R from the method selected in step S30 is then processed. t Using a sliding window with a length of 12 months and a step size of 12 months, the samples within a certain window are denoted as: {r1, r2, ..., r t ,…,r n }, n=12; Calculate each early warning indicator for each window separately, using the following formula: The formula for calculating AR(1) is as follows: ,in This represents the mean value within the window. The formula for calculating variance (Var) is as follows: ; Skewness is calculated as follows: , where s is the standard deviation of each value within the window; The formula for calculating kurtosis is as follows: .

5. The method for identifying critical transformation risk of an ecosystem based on multi-source remote sensing vegetation indicators according to claim 4, characterized in that: Kurtosis uses excessive kurtosis. excess =Kurtosis-3.

6. The method for identifying critical transformation risk of ecosystems based on multi-source remote sensing vegetation indicators according to claim 1, characterized in that: S40 also includes mapping the spatial distribution of risk based on the critical transformation risk index.

7. A system for identifying critical transformation risks in ecosystems based on multi-source remote sensing vegetation indicators, characterized in that, include: The vegetation time series data acquisition and preprocessing module is used to acquire and preprocess vegetation time series data of the study area. The data includes leaf area index (LAI), normalized difference vegetation index (NDVI), solar induced chlorophyll fluorescence (SIF), and vegetation optical depth (VOD). The preprocessing includes missing value imputation, outlier handling, and time scale unification. The signal separation module is used to separate signals based on various time series decomposition methods to obtain the residual time series R after removing long-term trends and seasonal cycles. t The various time series decomposition methods include the seasonal-trend decomposition method STL, the structural time series analysis method STA, and the state-space model method SSM. The method comparison and optimization module, based on the various time series decomposition methods and their residual time series described in the signal separation module, selects one method and its residual time series for subsequent analysis; The critical transition early warning statistical index calculation module is used to calculate the critical transition early warning statistical index, including the residual time series of the methods selected by the method comparison and optimization module. Sliding window method Calculate the statistical indicators for early warning of ecosystems to obtain the indicator sequence Z that changes over time. T The indicators include first-order autocorrelation coefficient AR(1), variance Var, skewness Skewness, and kurtosis. The linear trend analysis module is used to perform linear trend analysis on the time series of early warning indicators for each pixel. It uses the least squares method to calculate the slope b of the indicator value over time, which is used to characterize the direction and intensity of changes in ecosystem stability. b>0 indicates that the indicator increases over time, and the risk of stability decline increases; b<0 indicates that the indicator decreases over time, and stability recovers; b=0 indicates no significant change in the long term. The trend was tested for statistical significance, and the trend results that failed the significance test were weighted down to reduce the impact of random noise on the subsequent risk assessment results. Weight reduction processing includes reducing the weight of slope b, as shown in the following formula: Where w is the weight, which varies with significance p: ； The Consistency Score (CS) calculation module is used to directly sum the trend values ​​b* of each indicator obtained from the linear trend analysis module to obtain the consistency score CS at the pixel scale. The standard deviation (MS) module for Mahalanobis distance time series calculation is used to calculate the standard deviation of Mahalanobis distance time series based on the initial ecosystem state of the study area as a baseline state and the index series Z. T Construct a multi-index covariance matrix; Calculate the Mahalanobis distance of each pixel in the target period relative to the baseline state in the multi-index feature space, and further calculate the standard deviation MS of the Mahalanobis distance time series; The critical transformation risk index R module is used to construct the critical transformation risk index R and to determine the critical transformation risk of the ecosystem in the study area based on R. ； Where R(x,y) is the risk index of pixel (x,y), CS(x,y) is the multi-indicator consistency score of the pixel, and MS(x,y) is the standard deviation of Mahalanobis distance of the pixel; w1 and w2 are the fusion weights, min(CS) and max(CS) are the minimum and maximum values ​​of the consistency score of multiple indicators in the study area, respectively, and min(MS) and max(MS) are the minimum and maximum values ​​of the standard deviation of Mahalanobis distance in the study area, respectively.

Images