A method and system for evaluating urban-rural gradient vegetation resilience based on remote sensing data

Abstract

This invention discloses a method and system for assessing vegetation resilience across urban and rural gradients based on remote sensing data. The method includes: dividing the study area into urban and rural gradients based on a multi-period urban boundary and buffer zone expansion method; performing quality control, sensor consistency correction, and vegetation index calculation on remote sensing surface reflectance data; generating annual time series of vegetation resilience for each urban and rural gradient zone by calculating the first-order autocorrelation coefficient on a sliding window based on the deseasoned and detrended vegetation index residual sequence; identifying dynamic changes and turning points of vegetation resilience under different urban and rural gradients using trend testing and piecewise regression; and constructing a partial least squares regression model based on multi-source driving factors to determine the dominant driving factors of vegetation resilience changes under different urban and rural gradients. This invention enables continuous monitoring, dynamic identification, and trend attribution of vegetation resilience across urban and rural gradients, and can be used for urban and rural ecological risk early warning, ecological restoration assessment, and ecological governance decision support.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing, geographic information, and ecology, specifically to a method and system for assessing the resilience of urban and rural gradient vegetation based on remote sensing data. Background Technology

[0002] As urbanization continues, urban built-up areas are expanding into surrounding rural areas, creating a clear urban-rural gradient pattern among urban core areas, new urban areas, urban fringe areas, suburbs, rural fringe areas, and rural background areas. Under these different urban-rural gradients, vegetation ecosystems are affected by multiple factors, including changes in thermal environment, water conditions, land use transformation, and increased human disturbance, resulting in significant spatial differences in their growth status, stability, and resilience. Vegetation resilience, as a crucial indicator of a vegetation system's ability to resist external disturbances and maintain or restore its original functional state, is of great significance for urban ecological risk early warning, green space system optimization, ecological restoration assessment, and national land space ecological governance.

[0003] Existing urban vegetation monitoring methods largely rely on static or state-based indicators such as vegetation indices, vegetation cover, green space area, or land use type. While these indicators can reflect vegetation greenness and spatial distribution, they struggle to reveal the resilience and stability of vegetation systems under long-term disturbances. Although some methods utilize remote sensing time series analysis to study vegetation change trends, they typically focus on identifying greening or browning trends, lacking a comprehensive characterization of dynamic changes in vegetation resilience, phased transitions, and differences in urban-rural gradients. Furthermore, existing regional ecological resilience assessment methods do not fully consider the shaping effect of urban expansion on spatial gradients, making it difficult to distinguish the differences in vegetation resilience changes across urban and rural gradients.

[0004] Furthermore, vegetation resilience changes are driven by multiple factors, including climate background, climate variability, and the intensity of urbanization in neighboring areas. Current technologies typically use remote sensing vegetation time series, urban-rural gradient zoning, climate environmental factors, and land cover change information in a fragmented manner, lacking a unified spatial scale for calculating vegetation resilience, identifying trend reversals, and attributing dominant driving factors. Therefore, it is necessary to propose a method and system for assessing urban-rural gradient vegetation resilience based on remote sensing data. This would enable continuous monitoring, dynamic identification, and attribution of drivers of vegetation resilience under different urban-rural gradients, providing technical support for urban-rural ecological risk early warning, ecological restoration assessment, and ecological governance decision-making. Summary of the Invention

[0005] This invention addresses the shortcomings of existing technologies by providing a method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data, comprising the following steps: Collect urban boundary data from multiple periods and combine them with the buffer zone expansion method to divide the urban-rural gradient of the study area; Remote sensing surface reflectance data of the study area were collected, and the raw data were subjected to quality control and sensor consistency correction. The enhanced vegetation index (EVI) monthly time series was further obtained through band calculation. Based on environmental data of the study area, the annual mean and intra-annual variability were calculated and resampled to a spatial resolution consistent with the remote sensing vegetation index. At the same time, the urbanization intensity of the neighborhood was calculated based on the land cover data of the study area, and the annual mean, intra-annual variability and neighborhood urbanization intensity were used as multi-source driving factors. The generated monthly EVI time series were subjected to deseasoning, detrending, and sliding window lag first-order autocorrelation analysis to generate annual vegetation resilience time series for each urban-rural gradient zone. Based on the annual time series of vegetation resilience, the Mann-Kendall trend test and piecewise linear regression were used to identify the dynamic changes and turning points of vegetation resilience under different urban-rural gradients. By utilizing multiple driving factors and annual time series data on vegetation resilience, a partial least squares regression model was constructed to identify the dominant driving factors for changes in vegetation resilience under different urban-rural gradients.

[0006] Furthermore, urban boundary data for the study area in three representative historical years were obtained to characterize the phased changes in urban expansion, denoted as follows: U t1 , U t2 and U t3 , respectively representing the early city boundary, the city boundary at the stage of expansion, and the recent city boundary; among them, U t1 The internal area is divided into the city core area, U t1 outside, U t2 The area within this area will be designated as a new urban district. U t2 outside, U t3 The area within this zone is designated as the urban fringe; U t3 As the boundary of the urban area, a buffer zone is constructed to delineate the rural areas, ensuring that the first buffer zone meets the following requirements:

[0007] In the formula, Area(·) represents a function for calculating the area of ​​a spatial object; Ω represents the spatial extent of the study area; Ω \ U t3 Indicates the non-urban spatial extent of the study area; x Dist represents any spatial location outside the urban area; x , Ut3 )express x to U t3 The shortest distance; d Indicates the buffer distance; A u This indicates the recent urban area within the study region; based on the determined... d As a benchmark buffer distance for expanding rural areas outward from urban areas, the distance will be... U t3 0- d , d -2 d and 2 d -3 d The non-urban areas are divided into suburban areas, rural fringe areas, and rural background areas.

[0008] Furthermore, remote sensing surface reflectance data of the study area were collected. After quality control and sensor consistency correction, surface reflectance in the blue, red, and near-infrared bands was extracted, and the Enhanced Vegetation Index (EVI) was calculated. The EVI was calculated according to the following formula:

[0009] In the formula, r NIR , r Red and r Blue These represent the surface reflectance in the near-infrared, red, and blue light bands, respectively, after sensor consistency correction. G , C 1. C 2 and L The preset coefficients are used; the median of the EVI observations within the same month of the same year is synthesized using pixels as the unit to generate a pixel-by-pixel monthly EVI time series.

[0010] Furthermore, the environmental data includes temperature, precipitation, vapor pressure deficit, and soil moisture; each environmental variable is statistically analyzed annually, and its annual mean and intra-annual variability are calculated to represent the climate background factor and climate variability factor, respectively, and further resampled to a spatial resolution consistent with the Enhanced Vegetation Index (EVI) data; for any environmental variable... X , its first y Annual average X mean,y and intra-annual variability X cv,y They are represented as follows:

[0011]

[0012] In the formula, X y Representing environment variables X In the y The annual time series consists of multiple observations within the year, with Mean(·) representing the mean calculation function and Std(·) representing the standard deviation calculation function. Based on land cover data, the changes in impervious surfaces within the study area at different times are identified, and the proportion of pixels that transition from impervious to impervious surfaces within a preset neighborhood window is calculated as the neighborhood urbanization intensity; the neighborhood urbanization intensity... NUI y Represented as:

[0013] In the formula, N y,imp This indicates the number of pixels within a preset neighborhood window that transition from non-impermeable to impermeable relative to a reference year. N y,total This represents the total number of valid pixels within the preset neighborhood window.

[0014] Furthermore, the generated monthly EVI time series undergoes deseasoning and detrending processing to obtain the EVI residual time series. A sliding time window is set with a preset year as the time step, and the lagged first-order autocorrelation coefficient AC1 of the EVI residual time series is calculated within each sliding time window. The lagged first-order autocorrelation coefficient is used as a vegetation resilience characterization index to generate annual vegetation resilience time series for each urban-rural gradient zone. Specifically, for any given monthly EVI time series, the time... t The month to which it belongs is recorded as m(t) The multi-year average EVI value for that month is expressed as:

[0015] In the formula, EVI s Indicates time s EVI monthly value, m(s) Indicates time s The month to which it belongs, EVI m(t) Indicates time t The multi-year average EVI value for the same month; based on the multi-year average EVI value, the original monthly EVI time series is deseasoned to obtain the deseasoned EVI:

[0016] In the formula, EVI t Indicates timet EVI monthly value, EVI' t Indicates time t The deseasoned EVI is obtained; further, the deseasoned EVI is linearly detrended to obtain the trend term:

[0017] In the formula, Trending t Indicates time t The corresponding long-term trend item, α and β The coefficients are the linear trend fitting coefficients; based on the deseasoned EVI and the trend term, the EVI residual values ​​are calculated:

[0018] In the formula, R t Indicates time t The EVI residual values ​​are derived from the values ​​at each time point. R t The EVI residual time series are constructed in chronological order. For any given year y Corresponding sliding time window W y First-order autocorrelation coefficient AC1 y Represented as:

[0019] Where Corr(·) represents the correlation coefficient calculation function, R t Indicates time t The EVI residual value, R t 1 indicates its first-order lag residual value. W y Indicates the year y The corresponding sliding time window; AC1 y The calculation results are summarized according to urban-rural gradient zones to obtain the annual time series of vegetation resilience for each urban-rural gradient zone.

[0020] Furthermore, the Mann-Kendall trend test was used to calculate the trend statistic of the annual time series of vegetation resilience to determine the direction and significance of changes in vegetation resilience. S Represented as:

[0021] Where sign(·) represents the sign function,n Indicates the length of the annual time series of vegetation resilience. AC1 i and AC1 j They represent the first i Year and the j The annual vegetation resilience index value; when S A value greater than 0 indicates an upward trend in AC1, corresponding to a decrease in vegetation resilience; when... S When it is less than 0, it indicates that AC1 is decreasing, which corresponds to an increase in vegetation resilience; Piecewise linear regression was performed on the annual average vegetation resilience time series of each urban-rural gradient zone to determine the year that minimizes the piecewise regression residual as the turning point of the dynamic change in vegetation resilience; the piecewise linear regression is expressed as:

[0022] in, AC1 y Indicates the year y The vegetation resilience index value; β 0、 β 1 and β 2 is the regression coefficient; Y b Indicates the year of transition to be identified. Indicates when y Greater than Y b Time to take Otherwise, take 0; e y Represents the regression residuals; determined by minimizing the sum of squared residuals. Y b .

[0023] Furthermore, the multi-source driving factors are matched annually with the annual time series of vegetation resilience, and the multi-source driving factors are standardized. Using the annual time series of vegetation resilience as the response variable, and climate background factors (i.e., annual mean), climate variability factors (i.e., intra-annual variability), and neighborhood urbanization intensity as explanatory variables, a partial least squares regression model is constructed. Based on the partial least squares regression model, the variable projection importance values ​​of each explanatory variable are calculated, and the explanatory variable with the highest variable projection importance value is determined as the dominant driving factor of vegetation resilience change. The direction of the dominant driving factor's effect on vegetation resilience change is determined according to the corresponding regression coefficient sign. Furthermore, the types and directions of the dominant driving factors are statistically analyzed according to urban-rural gradient zones to obtain the attribution results of vegetation resilience change under different urban-rural gradients. The partial least squares regression model is expressed as follows:

[0024] in,Y This represents the annual time series of vegetation resilience. X This represents the matrix of explanatory variables consisting of climate background factors, climate variability factors, and the intensity of urbanization in the surrounding area. B Represents the regression coefficient matrix. E This represents the residual matrix.

[0025] Furthermore, after constructing the partial least squares regression model, the variable projection importance values ​​of each explanatory variable are calculated. For the , j Each explanatory variable has a projected importance value. VIP j Represented as:

[0026] In the formula, p Indicates the number of explanatory variables. A This indicates the number of latent variables extracted by the partial least squares regression model. SSY a Indicates the first a One latent variable to the response variable Y Explanation of the sum of squares, w ja Indicates the first j The explanatory variable in the th... a The weights on each latent variable will VIP The largest explanatory variable was identified as the dominant driving factor of vegetation resilience change, and its direction of action was determined based on the sign of the regression coefficient corresponding to the explanatory variable. When the regression coefficient of the explanatory variable is positive, it indicates that it changes in the same direction as the lagged first-order autocorrelation coefficient AC1, which corresponds to a decrease in vegetation resilience. When the regression coefficient of the explanatory variable is negative, it indicates that it changes in the opposite direction to AC1, which corresponds to an increase in vegetation resilience.

[0027] The present invention also provides a system for assessing the resilience of urban and rural gradient vegetation based on remote sensing data, comprising: a processor and a memory, wherein the memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute a method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data as described in the above technical solution.

[0028] The present invention also provides a computer-readable storage medium, comprising: a readable storage medium storing a computer program, wherein when the computer program is executed, it implements a method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data as described in the above technical solution.

[0029] Compared with the prior art, the present invention has the following advantages: This invention uses urban-rural gradient zoning as an analytical framework, combining urban expansion processes with vegetation resilience assessment. It can distinguish spatial differences in vegetation resilience under different urbanization intensities. By performing deseasoning, detrending, and sliding window lag first-order autocorrelation analysis on the time series of remote sensing vegetation indices, continuous and dynamic monitoring of vegetation resilience can be achieved, overcoming the limitations of traditional greenness assessments in reflecting ecosystem recovery capacity. Furthermore, by combining trend testing, piecewise regression, and partial least squares regression models, it can identify trends, turning points, and dominant driving factors in vegetation resilience changes, improving the interpretability and application value of the assessment results. This provides technical support for urban and rural ecological risk early warning, ecological restoration, and ecological governance decision-making. Attached Figure Description

[0030] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0031] Figure 1 This is a flowchart of the urban-rural gradient vegetation resilience assessment method according to an embodiment of the present invention.

[0032] Figure 2 This is a spatial distribution map of urban and rural gradients in the study area of ​​this invention embodiment.

[0033] Figure 3 This is a graph showing the interannual variation of vegetation resilience along the urban-rural gradient in the study area of ​​this invention embodiment.

[0034] Figure 4 This is a map showing the area proportion of the dominant driving factors for vegetation resilience trends in the study area of ​​this invention embodiment. Detailed Implementation

[0035] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings and embodiments. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0036] Example 1 like Figure 1 As shown, this embodiment of the invention provides a method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data, comprising the following steps: Step 1: Collect urban boundary data from multiple periods and combine it with the buffer expansion method to divide the urban-rural gradient of the study area.

[0037] This embodiment takes the Yangtze River Delta urban agglomeration as the study area, and extracts urban boundary data for the study area in 1990, 2000, and 2018 from the Global UrbanBoundaries (GUB) dataset, denoted as follows: U 1990 , U 2000 and U 2018 Divide the city's internal gradient zones according to the time sequence of urban expansion: U 1990 The coverage area is divided into the urban core area, U 2000 Compared to U 1990 The newly added urban area is designated as a new urban district. U 2018 Compared to U 2000 The newly added urban area is designated as the urban fringe.

[0038] Furthermore, with U 2018 As the boundary of the urban area, a buffer zone is constructed outwards to delineate the rural areas, ensuring that the first buffer zone meets the following requirements:

[0039] In the formula, Area(·) represents a function for calculating the area of ​​a spatial object; Ω represents the spatial extent of the study area; Ω \ U 2018 Indicates the non-urban spatial extent of the study area; x Dist represents any spatial location outside the urban area; x , U 2018 )express x to U 2018 The shortest distance; d Indicates the buffer distance; A u This indicates the recent urban area within the study region. Based on the determined... d As a benchmark buffer distance for expanding rural areas outward from urban areas, the distance will be... U 2018 0- d , d -2 d and 2 d -3 d The non-urban areas are divided into suburban areas, rural fringe areas, and rural background areas. Figure 2This is a spatial distribution map of the urban-rural gradient in the study area.

[0040] Step 2: Collect remote sensing surface reflectance data of the study area, perform quality control and sensor consistency correction on the raw data, and further obtain the monthly time series of Enhanced Vegetation Index (EVI) through band calculation.

[0041] In this embodiment, the remote sensing surface reflectance data uses 30 m spatial resolution surface reflectance images acquired by Landsat 5, Landsat 7, and Landsat 8 sensors from 2004 to 2024, and is processed in batches on the Google Earth Engine cloud platform. First, the remote sensing images covering the study area are quality controlled, and images with cloud cover greater than 80% or solar zenith angle greater than 60° are removed. For the remote sensing images retained after image-level screening, the quality of pixels is further judged based on the quality control bands inherent in the images. Low-quality pixels marked as clouds, cloud shadows, snow, fill values, radiation saturation, or other abnormal observation states are identified and removed, and only clear-sky effective pixels that are not affected by the above factors and whose reflectance in each target band is within the effective range are retained.

[0042] Because the Operational Land Imager (OLI) sensor on the Landsat 8 differs in spectral response from the Thematic Mapper™ sensor on the Landsat 5 and the Enhanced Thematic Mapper Plus (ETM+) sensor on the Landsat 7, this embodiment further performs consistency correction on the surface reflectance data acquired by different sensors. Specifically, using the linear conversion relationship between the surface reflectance of the Landsat 7 ETM+ and the Landsat 8 OLI disclosed in existing literature, and taking the Landsat 7 ETM+ as a reference, the surface reflectance of the Landsat 8 OLI in the blue, red, and near-infrared bands is converted into reflectance comparable to the corresponding bands of the ETM+. The conversion relationship is as follows:

[0043]

[0044]

[0045] in, r Blue,OLI , r Red,OLI and r NIR,OLI These represent the surface reflectance in the blue, red, and near-infrared bands acquired by the Landsat 8 OLI sensor, respectively. rBlue , r Red and r NIR These represent the surface reflectance of the corresponding band after sensor consistency correction.

[0046] After consistency correction, the surface reflectance in the blue, red, and near-infrared bands is extracted, and the Enhanced Vegetation Index (EVI) is calculated using the following formula:

[0047] In the formula, G , C 1. C 2 and L This is a preset coefficient. In this embodiment, G Take 2.5, C Take 6 from 1. C 2. Take 7.5. L Take 1. After calculating the EVI based on the effective clear sky observations, the median of the EVI observations within the same month of the same year is synthesized in pixels to generate a pixel-by-pixel monthly time series of EVI from 2004 to 2024.

[0048] Step 3: Collect environmental and land cover data for the study area, process temperature, precipitation, vapor pressure deficit and soil moisture into annual mean and intra-annual variability, and resample to a spatial resolution consistent with the remote sensing vegetation index. At the same time, calculate the urbanization intensity of the neighborhood based on the land cover data.

[0049] This embodiment uses environmental data from Terraclimate and SMCI 1.0 to statistically analyze temperature, precipitation, vapor pressure deficit, and soil moisture on an annual basis. The annual mean and intra-annual variability for 2004–2024 are calculated, representing climate background factors and climate variability factors, respectively. The data is further resampled to a spatial resolution consistent with the Enhanced Vegetation Index (EVI) data. For any environmental variable… X , its first y Annual average X mean,y and intra-annual variability X cv,y They are represented as follows:

[0050]

[0051] In the formula, X y Representing environment variables X In the yThe annual time series consists of multiple observations within the year, with Mean(·) representing the mean calculation function and Std(·) representing the standard deviation calculation function. This embodiment uses land cover data from the China Land Cover Dataset (CLCD), with 2004 as the base year, and identifies, pixel by pixel, the land cover changes from impermeable to impermeable surfaces in subsequent years. When calculating the neighborhood urbanization intensity, a neighborhood window is constructed centered on each target pixel. In this embodiment, the neighborhood window is a 3×3 pixel moving window. The proportion of pixels within the moving window that change from impermeable to impermeable surfaces is calculated as the neighborhood urbanization intensity; this neighborhood urbanization intensity... NUI y Represented as:

[0052] In the formula, N y,imp This represents the number of pixels within the neighborhood window that transition from impermeable to impermeable relative to the reference year. N y,total This represents the total number of valid pixels within the neighborhood window.

[0053] Step 4: Perform deseasoning, detrending, and sliding window lag first-order autocorrelation analysis on the monthly EVI time series generated in Step 2 to generate annual vegetation resilience time series for each urban and rural gradient zone.

[0054] In this embodiment, the monthly EVI time series generated in step 2 is deseasonalized to eliminate the influence of monthly seasonal cycles; subsequently, the deseasonalized EVI time series is detrended to eliminate the influence of long-term trends, resulting in the EVI residual time series used for vegetation resilience calculation. Specifically, for any monthly EVI time series, the time... t The month to which it belongs is recorded as m(t) The multi-year average EVI value for that month is expressed as:

[0055] In the formula, EVI s Indicates time s EVI monthly value, m(s) Indicates time s The month to which it belongs, EVI m(t) Indicates time t The multi-year average EVI value for the same month. Based on the multi-year average EVI value, the original monthly EVI time series is deseasoned to obtain the deseasoned EVI:

[0056] In the formula, EVI t Indicates time t EVI monthly value, EVI' t Indicates time t The deseasonalized EVI. Further linear detrending processing is performed on the deseasonalized EVI to obtain the trend term:

[0057] In the formula, Trending t Indicates time t The corresponding long-term trend item, α and β The coefficients represent the linear trend fitting coefficients. Based on the deseasoned EVI and the trend term, the EVI residuals are calculated:

[0058] In the formula, R t Indicates time t The EVI residual values. (Based on data from each time point.) R t The EVI residual time series are constructed in chronological order.

[0059] In this embodiment, a 5-year sliding time window and a 1-year sliding step are used to calculate the lag first-order autocorrelation coefficient AC1 of the EVI residual time series within each sliding time window. This lag first-order autocorrelation coefficient is then used as a characterization index of vegetation resilience. For any given year... y Corresponding sliding time window W y First-order autocorrelation coefficient AC1 y Represented as:

[0060] Where Corr(·) represents the correlation coefficient calculation function, R t Indicates time t The EVI residual value, R t 1 indicates its first-order lag residual value. W y Indicates the year y The corresponding sliding time window. AC1 ,y The calculation results are summarized according to urban-rural gradient zones to obtain the annual time series of vegetation resilience for each urban-rural gradient zone.

[0061] Step 5: Based on the vegetation resilience time series generated in Step 4, the Mann-Kendall trend test and piecewise linear regression are used to identify the dynamic changes and turning points of vegetation resilience under different urban-rural gradients.

[0062] The Mann-Kendall trend test was used to calculate the trend statistic of the annual time series of vegetation resilience to determine the direction and significance of changes in vegetation resilience. S Represented as:

[0063] Where sign(·) represents the sign function, n Indicates the length of the annual time series of vegetation resilience. AC1 i and AC1 j They represent the first i Year and the j The annual vegetation resilience index value; when S A value greater than 0 indicates an upward trend in AC1, corresponding to a decrease in vegetation resilience; when... S When it is less than 0, it indicates that AC1 is decreasing, which corresponds to an increase in vegetation resilience.

[0064] Further, piecewise linear regression was performed on the annual average vegetation resilience time series of each urban-rural gradient zone to determine the year with the smallest piecewise regression residual as the turning point of the dynamic change in vegetation resilience; the piecewise linear regression is expressed as:

[0065] in, AC1 y Indicates the year y The vegetation resilience index value; β 0、 β 1 and β 2 is the regression coefficient; Y b Indicates the year of transition to be identified. Indicates when y Greater than Y b Time to take Otherwise, take 0; e y Represents the regression residuals; determined by minimizing the sum of squared residuals. Y b . Figure 3The figure shows the interannual variation curves of vegetation resilience along the urban-rural gradient in the study area. AC1 is inversely related to vegetation resilience. The vegetation resilience of the suburban gradient consistently changed from an increasing trend to a decreasing trend, with the turning point around 2012.

[0066] Step 6: Using the multi-source driving factors collected in Step 3 (including climate background factors, climate variability factors and neighboring urbanization intensity) and the vegetation resilience time series generated in Step 4, construct a partial least squares regression model to determine the dominant driving factors of vegetation resilience changes under different urban-rural gradients.

[0067] Multi-source driving factor data were matched annually with vegetation resilience time series, and the multi-source driving factors were standardized. A partial least squares regression model was constructed using the annual time series of vegetation resilience as the response variable and climate background factors, climate variability factors, and neighboring urbanization intensity as explanatory variables. Based on the partial least squares regression model, the variable projection importance values ​​of each explanatory variable were calculated. The explanatory variable with the highest variable projection importance value was identified as the dominant driving factor for vegetation resilience change, and the direction of the dominant driving factor's effect on vegetation resilience change was determined according to the corresponding regression coefficient sign. Furthermore, the types and directions of the dominant driving factors were statistically analyzed according to urban-rural gradient zones to obtain the attribution results of vegetation resilience change under different urban-rural gradients. The partial least squares regression model is expressed as follows:

[0068] in, Y This represents the annual time series of vegetation resilience. X This represents the matrix of explanatory variables consisting of climate background factors, climate variability factors, and the intensity of urbanization in the surrounding area. B Represents the regression coefficient matrix. E Let represent the residual matrix. After constructing the partial least squares regression model, the variable projection importance values ​​of each explanatory variable are further calculated. For the , j Each explanatory variable has a projected importance value. VIP j Represented as:

[0069] In the formula, p Indicates the number of explanatory variables. A This indicates the number of latent variables extracted by the partial least squares regression model. SSY a Indicates the first a One latent variable to the response variable Y Explanation of the sum of squares, w ja Indicates the first j The explanatory variable in the th... aWeights on each latent variable. VIP j The larger the value, the higher the value. j The greater the contribution of each explanatory variable to the model's explanatory power, the better. VIP The largest explanatory variable was identified as the dominant driving factor of vegetation resilience changes, and its direction of action was determined based on the sign of its regression coefficient. When the regression coefficient of the explanatory variable is positive, it indicates that it changes in the same direction as the lagged first-order autocorrelation coefficient AC1, corresponding to a decrease in vegetation resilience; when the regression coefficient of the explanatory variable is negative, it indicates that it changes in the opposite direction to AC1, corresponding to an increase in vegetation resilience. Figure 4 The map shows the area proportion of the dominant driving factors for vegetation resilience trends in the study area. The annual average temperature and the annual average soil moisture are the two most important dominant driving factors for vegetation resilience trends, accounting for 41.3% and 15.2% of the vegetation resilience changes, respectively.

[0070] Example 2 Based on the same inventive concept, the present invention also provides an urban-rural gradient vegetation resilience assessment system based on remote sensing data, including a processor and a memory. The memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute the urban-rural gradient vegetation resilience assessment method based on remote sensing data as described above.

[0071] Example 3 Based on the same inventive concept, the present invention also provides a computer-readable storage medium, including a readable storage medium on which a computer program is stored, wherein when the computer program is executed, it implements the above-described method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data.

[0072] In specific implementation, the method proposed in the technical solution of this invention can be automatically executed by those skilled in the art using computer software technology. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of this invention and computer equipment including the computer program running the corresponding computer program, should also be within the protection scope of this invention.

[0073] The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains may make various modifications or additions to the described specific embodiments or use similar methods to replace them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A method for assessing vegetation resilience based on remote sensing data in urban-rural gradient, characterized in that, include: Collect urban boundary data from multiple periods and combine them with the buffer zone expansion method to divide the urban-rural gradient of the study area; Remote sensing surface reflectance data of the study area were collected, and the raw data were subjected to quality control and sensor consistency correction. The enhanced vegetation index (EVI) monthly time series was further obtained through band calculation. Based on environmental data of the study area, the annual mean and intra-annual variability were calculated and resampled to a spatial resolution consistent with the remote sensing vegetation index. At the same time, the urbanization intensity of the neighborhood was calculated based on the land cover data of the study area, and the annual mean, intra-annual variability and neighborhood urbanization intensity were used as multi-source driving factors. The generated monthly EVI time series were subjected to deseasoning, detrending, and sliding window lag first-order autocorrelation analysis to generate annual vegetation resilience time series for each urban-rural gradient zone. Based on the annual time series of vegetation resilience, the Mann-Kendall trend test and piecewise linear regression were used to identify the dynamic changes and turning points of vegetation resilience under different urban-rural gradients. By utilizing multiple driving factors and annual time series data on vegetation resilience, a partial least squares regression model was constructed to identify the dominant driving factors for changes in vegetation resilience under different urban-rural gradients.

2. The method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data as described in claim 1, characterized in that: We obtained urban boundary data for the study area in three representative historical years to characterize the phased changes in urban expansion, denoted as follows: U t1 , U t2 and U t3 , respectively representing the early city boundary, the city boundary at the stage of expansion, and the recent city boundary; among them, U t1 The internal area is divided into the city core area, U t1 outside, U t2 The area within this area will be designated as a new urban district. U t2 outside, U t3 The area within this zone is designated as the urban fringe; U t3 As the boundary of the urban area, a buffer zone is constructed to delineate the rural areas, ensuring that the first buffer zone meets the following requirements: In the formula, Area(·) represents a function for calculating the area of ​​a spatial object; Ω represents the spatial extent of the study area; Ω \ U t3 Indicates the non-urban spatial extent of the study area; x Dist represents any spatial location outside the urban area; x , U t3 )express x to U t3 The shortest distance; d Indicates the buffer distance; A u This indicates the recent urban area within the study region; based on the determined... d As a benchmark buffer distance for expanding rural areas outward from urban areas, the distance will be... U t3 0- d , d -2 d and 2 d -3 d The non-urban areas are divided into suburban areas, rural fringe areas, and rural background areas.

3. The method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data as described in claim 1, characterized in that: Remote sensing surface reflectance data of the study area were collected. After quality control and sensor consistency correction, surface reflectance in the blue light band, red light band and near-infrared band was extracted, and the Enhanced Vegetation Index (EVI) was calculated. The EVI is calculated according to the following formula: In the formula, ρ NIR , ρ Red and ρ Blue These represent the surface reflectance in the near-infrared, red, and blue light bands, respectively, after sensor consistency correction. G , C 1. C 2 and L The preset coefficients are used; the median of the EVI observations within the same month of the same year is synthesized using pixels as the unit to generate a pixel-by-pixel monthly EVI time series.

4. The method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data as described in claim 1, characterized in that: The environmental data includes temperature, precipitation, vapor pressure deficit, and soil moisture; each environmental variable is statistically analyzed annually, and its annual mean and intra-annual variability are calculated to represent the climate background factor and climate variability factor, respectively, and further resampled to a spatial resolution consistent with the Enhanced Vegetation Index (EVI) data; for any environmental variable... X , its first y Annual average X mean,y and intra-annual variability X cv,y They are represented as follows: In the formula, X y Representing environment variables X In the y The annual time series consists of multiple observations within the year, with Mean(·) representing the mean calculation function and Std(·) representing the standard deviation calculation function. Based on land cover data, the changes in impervious surfaces within the study area at different times are identified, and the proportion of pixels that transition from impervious to impervious surfaces within a preset neighborhood window is calculated as the neighborhood urbanization intensity; the neighborhood urbanization intensity... NUI y Represented as: In the formula, N y,imp This indicates the number of pixels within a preset neighborhood window that transition from non-impermeable to impermeable relative to a reference year. N y,total This represents the total number of valid pixels within the preset neighborhood window.

5. The method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data as described in claim 1, characterized in that: The generated monthly EVI time series are deseasoned and detrended to obtain the EVI residual time series. A sliding time window is set with a preset year as the time step, and the lagged first-order autocorrelation coefficient AC1 of the EVI residual time series is calculated within each sliding time window. The lagged first-order autocorrelation coefficient is used as a vegetation resilience characterization index to generate annual vegetation resilience time series for each urban-rural gradient zone. For any given monthly EVI time series, the time... t The month to which it belongs is recorded as m(t) The multi-year average EVI value for that month is expressed as: In the formula, EVI s Indicates time s EVI monthly value, m(s) Indicates time s The month to which it belongs, EVI m(t) Indicates time t The multi-year average EVI value for the same month; based on the multi-year average EVI value, the original monthly EVI time series is deseasoned to obtain the deseasoned EVI: In the formula, EVI t Indicates time t EVI monthly value, EVI' t Indicates time t The deseasoned EVI is obtained; further, the deseasoned EVI is linearly detrended to obtain the trend term: In the formula, Trend t Indicates time t The corresponding long-term trend item, α and β The coefficients are the linear trend fitting coefficients; based on the deseasoned EVI and the trend term, the EVI residual values ​​are calculated: In the formula, R t Indicates time t The EVI residual values ​​are derived from the values ​​at each time point. R t The EVI residual time series are constructed in chronological order. For any given year y Corresponding sliding time window W y First-order autocorrelation coefficient AC1 y Represented as: Where Corr(·) represents the correlation coefficient calculation function, R t Indicates time t The EVI residual value, R t 1 indicates its first-order lag residual value. W y Indicates the year y The corresponding sliding time window; AC1 y The calculation results are summarized according to urban-rural gradient zones to obtain the annual time series of vegetation resilience for each urban-rural gradient zone.

6. The method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data as described in claim 1, characterized in that: The Mann-Kendall trend test was used to calculate the trend statistic of the annual time series of vegetation resilience to determine the direction and significance of changes in vegetation resilience. S Represented as: Where sign(·) represents the sign function, n Indicates the length of the annual time series of vegetation resilience. AC1 i and AC1 j They represent the first i Year and the j The annual vegetation resilience index value; when S A value greater than 0 indicates an upward trend in AC1, corresponding to a decrease in vegetation resilience; when... S When it is less than 0, it indicates that AC1 is decreasing, which corresponds to an increase in vegetation resilience; Piecewise linear regression was performed on the annual average vegetation resilience time series of each urban-rural gradient zone to determine the year that minimizes the piecewise regression residual as the turning point of the dynamic change in vegetation resilience; the piecewise linear regression is expressed as: in, AC1 y Indicates the year y The vegetation resilience index value; β 0、 β 1 and β 2 is the regression coefficient; Y b Indicates the year of transition to be identified. Indicates when y Greater than Y b Time to take Otherwise, take 0; ε y Represents the regression residuals; determined by minimizing the sum of squared residuals. Y b .

7. The method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data as described in claim 1, characterized in that: The multi-source driving factors were matched annually with the annual time series of vegetation resilience, and the multi-source driving factors were standardized. The annual time series of vegetation resilience was used as the response variable, and the climate background factor (annual mean), the climate variability factor (intra-annual variability), and the intensity of urbanization in the neighborhood were used as explanatory variables to construct a partial least squares regression model. Based on the partial least squares regression model, the projected importance values ​​of each explanatory variable are calculated. The explanatory variable with the highest projected importance value is identified as the dominant driving factor for vegetation resilience change. The direction of the dominant driving factor's effect on vegetation resilience change is determined according to the sign of the corresponding regression coefficient. Furthermore, the types and directions of the dominant driving factors are statistically analyzed according to urban-rural gradient zones to obtain the attribution results of vegetation resilience change under different urban-rural gradients. The partial least squares regression model is expressed as follows: in, Y This represents the annual time series of vegetation resilience. X This represents the matrix of explanatory variables consisting of climate background factors, climate variability factors, and the intensity of urbanization in the surrounding area. B Represents the regression coefficient matrix. E This represents the residual matrix.

8. The method for assessing the resilience of urban and rural gradient vegetation based on remote sensing data as described in claim 7, characterized in that: After constructing the partial least squares regression model, the variable projection importance values ​​of each explanatory variable are further calculated. For the , j Each explanatory variable has a projected importance value. VIP j Represented as: In the formula, p Indicates the number of explanatory variables. A This indicates the number of latent variables extracted by the partial least squares regression model. SSY a Indicates the first a One latent variable to the response variable Y Explanation of the sum of squares, w ja Indicates the first j The explanatory variable in the th... a The weights on each latent variable will VIP The largest explanatory variable was identified as the dominant driving factor of vegetation resilience change, and its direction of action was determined based on the sign of the regression coefficient corresponding to this explanatory variable. When the regression coefficient of the explanatory variable is positive, it indicates that it changes in the same direction as the lagged first-order autocorrelation coefficient AC1, which drives a decrease in vegetation resilience; when the regression coefficient of the explanatory variable is negative, it indicates that it changes in the opposite direction to AC1, which drives an increase in vegetation resilience.

9. A system for assessing the resilience of urban and rural gradient vegetation based on remote sensing data, characterized in that, include: The processor and memory, wherein the memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute the urban and rural gradient vegetation resilience assessment method based on remote sensing data as described in any one of claims 1-8.

10. A computer-readable storage medium, characterized in that, include: A readable storage medium storing a computer program, which, when executed, implements the method for assessing urban and rural gradient vegetation resilience based on remote sensing data as described in any one of claims 1-8.

Images