A method and system for constructing a vegetation drought resistance index

Abstract

The application discloses a vegetation drought resistance index construction method and system, which comprises the following steps: acquiring multi-source time series data and unifying spatial and time resolutions, constructing a precipitation comprehensive index, a runoff comprehensive index, a soil moisture index and an underground water drought index, and generating a vertical composite drought index VCDI by fusing the R-Vine Copula model; identifying a drought event and extracting an event duration, severity and intensity; performing normalization processing on a kernel normalized difference vegetation index, a land surface temperature and the VCDI, constructing a three-dimensional feature space, and defining and calculating a vegetation drought resistance index VDRI by using a Euclidean distance; calculating vegetation drought recovery time and recovery force by using VDRI abnormal values and standardized abnormal values, and identifying a time lag influence of drought on vegetation ecological functions by using lag cross-correlation analysis. The method accurately depicts a drought process and a comprehensive influence of the drought process on vegetation, fills a research blank in the field of comprehensive evaluation of drought-vegetation response, and provides a scientific basis for ecological restoration and drought adaptation management.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of ecological remote sensing and hydro-meteorological monitoring, and particularly relates to a vegetation drought resistance index construction method and system. BACKGROUND

[0002] As a complex natural disaster, the influence of drought spreads along the hydrological cycle path, and is manifested in turn as insufficient precipitation (meteorological drought), reduced runoff (hydrological drought), and decreased soil moisture (agricultural drought). Subsequently, drought conditions will cause stress to vegetation growth and damage to ecosystem functions, leading to ecological drought. As a basic component of the terrestrial ecosystem, the health of vegetation is directly constrained by the availability of water resources. Therefore, elucidating the mechanism of the influence of drought on vegetation and accurately assessing the response and resistance of vegetation to drought have become the key to current research. Therefore, developing a precise monitoring method for vegetation drought resistance is crucial for ensuring regional ecological safety and developing effective climate change adaptation strategies in both theoretical and practical aspects. SUMMARY

[0003] The present application provides a vegetation drought resistance index construction method, which comprises the following steps:

[0004] Step 1: Obtain multi-source time series data of the study area and unify them to the same spatial and temporal resolution;

[0005] Step 2: Construct the precipitation comprehensive index xSPI, the runoff comprehensive index xSRI, the soil moisture index SWI, the groundwater drought index GGDI, and the normalized difference vegetation index kNDVI;

[0006] Step 3: Calculate the edge distribution of xSPI, xSRI, SWI, and GGDI, use the Copula model to fit the time series grid data to calculate the vertical composite drought index VCDI, evaluate the goodness of fit of the Copula fitting results of different structures, and select the optimal model to obtain the monthly scale vertical composite drought index VCDI;

[0007] Step 4: Obtain a set of drought events, and calculate the number of drought events, the duration of the events, the severity, and the intensity;

[0008] Step 5: According to the drought events identified in step 4, further calculate the vegetation drought resistance index within the time range of each drought event, normalize kNDVI, land surface temperature LST, and VCDI in the same direction, and define the vegetation drought resistance index VDRI based on the Euclidean distance of the three-dimensional feature space. The larger the VDRI, the stronger the vegetation drought resistance, and the smaller the VDRI, the weaker the vegetation drought resistance.

[0009] Further, the multi-source time series data in step 1 includes precipitation, runoff, soil moisture SM, land water storage anomaly TWSA based on GRACE satellite and GLDAS land water elements, normalized difference vegetation index kNDVI, land surface temperature LST.

[0010] Further, in step 2, the Gamma distribution is used to fit the precipitation sequence and standardized to obtain the multi-time scale standardized runoff index SRI, and the precipitation comprehensive index xSPI is obtained based on principal component analysis PCA fusion, and the calculation formula is as follows:

[0011] (1)

[0012] (2)

[0013] In the formula, wherein represents the probability density function PDF of represents the Gamma function of , and and respectively represent the shape parameter and the scale parameter, and both are greater than 0; represents the precipitation in the precipitation time series data; represents the cumulative distribution function CDF of the PDF fitted to the time series precipitation data; represents the inverse cumulative distribution function of the standard normal distribution;

[0014] Similarly, the multi-time scale standardized runoff index SRI is calculated for the runoff sequence, and the runoff comprehensive index xSRI is obtained by PCA fusion, and the calculation formula is as follows:

[0015] (3)

[0016] represents the cumulative distribution function CDF of the PDF fitted to the time series precipitation data.

[0017] Further, in step 2, the soil moisture index SWI is constructed by exponential smoothing for the monthly scale soil moisture, and the calculation formula is as follows:

[0018] (4)

[0019] In the formula, represents the soil water content index of the th month, represents the soil water content of the th month, represents the current month when calculating the soil water content index SWI, and represents the soil water content of a certain month in the past,​ denotes a time lag constant, which is used to reflect the time scale, i is an integer index, which is used to uniquely identify each historical month participating in the calculation of the current month soil moisture index;

[0020] After removing the changes of soil, snow, vegetation canopy and surface water, the groundwater storage anomaly GWSA is obtained, and the groundwater drought index GGDI is calculated, and the calculation formula is as follows:

[0021] (5)

[0022] (6)

[0023] (7)

[0024] In the formula, the soil moisture storage anomaly SMSA, the snow water equivalent anomaly SWEA, the canopy water storage anomaly CWSA and the surface water storage anomaly SWSA are all from the global land data assimilation system GLDAS land surface model, denotes the groundwater storage anomaly GWSA value of the i-th month of the j-th year, denotes the average value of GWSA of the j-th month of all years, denotes the groundwater storage deviation GSD of the i-th month of the j-th year; and denote the average value and the standard deviation of the GSD of the i-th month of the j-th year, respectively, denotes the GRACE groundwater drought index of the i-th month of the j-th year; The kernel normalized difference vegetation index kNDVI is constructed by using the Gaussian kernel function, and the specific calculation formula is as follows: (8) In the formula, denotes the sensitivity of the normalized difference vegetation index to sparse or dense vegetation areas, denotes the near-infrared band data of the remote sensing image,

[0025] denotes the red band data of the remote sensing image.

[0026] (9)

[0027] In the formula, denotes the sensitivity of the normalized difference vegetation index to sparse or dense vegetation areas, denotes the near-infrared band data of the remote sensing image, denotes the red band data of the remote sensing image.

[0028] Further, the Copula model is used to fit the time series grid data in step 3 to calculate VCDI, and for random variables and their joint cumulative distribution function ​​​, there is a Copula function satisfies:

[0029] (10)

[0030] where n represents the number of variables; the goodness of fit (GOF) of the Copula model is evaluated by the root mean square error (RMSE), Akaike information criterion (AIC) and Bayesian information criterion (BIC), and the formulas are as follows:

[0031] (11)

[0032] (12)

[0033] (13)

[0034] wherein, represents the maximum likelihood value of the Copula model, represents the number of parameters in the model, represents the sample size used to fit the Copula model; and respectively represent the fitted value and the fitted expected value, and after comparing multiple Copula models, the model with the highest goodness of fit is used.

[0035] Further, in step 4, the drought threshold is set according to the value of VCDI, and the threshold is set as , and are constants, and if the vertical compound drought index value of the t-th month satisfies , it is recorded as a potential drought event; for a drought event with a duration of 1 month, if , it is considered as an independent drought event and is retained, otherwise it is not counted; for adjacent drought events with an interval of 1 month, if , they will be combined into one drought event;

[0036] The specific calculation formulas of the number of drought events, event duration, severity and intensity are as follows:

[0037] (14)

[0038] (15)

[0039] (16)

[0040] wherein, represents the total duration of the identified drought events, Total number of drought events, Duration of drought events , and Index of the start time and end time of drought events , respectively; Total severity of all drought events in a certain period, Severity of drought events , is an indicator function, taking the value 1 when , otherwise taking the value 0, Drought intensity in a certain period.

[0041] Further, the calculation formula of the vegetation drought resistance index is as follows:

[0042] (17)

[0043] (18)

[0044] (19)

[0045] In the formula, , and respectively represent the normalized monthly kNDVI, LST and VCDI, and the subscripts min and max respectively represent the minimum value and the maximum value; by establishing a three-dimensional feature space of , and , the VDRI is calculated using the Euclidean distance, and the calculation formula is as follows:

[0046] (20)

[0047] In the formula, in the three-dimensional feature space , and respectively represent the maximum value of , the minimum value of and the minimum value of .

[0048] Further, it also includes calculating the VDRI anomaly value A VDRI to reflect the change of the vegetation drought resistance, using the VDRI standardized anomaly value SA VDRI to determine the time required for the vegetation drought resistance to recover to the best state, and the vegetation drought recovery force is quantified by the difference between the SA VDRI of the most severe drought month and the best recovery month, and the formula is as follows:

[0049] (21)

[0050] (22)

[0051] (23)

[0052] wherein, represents the recovery of vegetation after experiencing drought, represents the A VDRI of the first month, represents the SA VDRI of the first month, represents the VDRI value of the first month, and respectively represent the mean and standard deviation of VDRI in the time range of selected data, and respectively represent the SA VDRI of the month when drought is most serious and the month when recovery is best.

[0053] Further, the method further comprises identifying the lag relationship between VDRI and the vegetation ecological function index by using lag cross-correlation, and calculating the lag response month of the vegetation ecological function index to VDRI, and the formula is as follows:

[0054] (24)

[0055] (25)

[0056] (26)

[0057] wherein, represents the maximum lag month, represents the maximum lag correlation between VDRI and the vegetation ecological function index after the kth month, and the greater the absolute value of the numerical value, the higher the correlation between VDRI and the vegetation ecological function index, represents the covariance between and , and are respectively the standard deviations of the two sequences, and respectively represent the lag correlation coefficient and the lag time.

[0058] The application further provides a vegetation drought resistance index construction system for implementing the vegetation drought resistance index construction method.

[0059] Alternatively, the system comprises a readable storage medium, and the readable storage medium stores a computer program.

[0060] Compared with the prior art, the application has the following advantages:

[0061] (1) The application unifies multiple drought indexes such as meteorology (xSPI), runoff (xSRI), soil moisture (SWI), and groundwater (GGDI) to the same spatial and temporal resolution, and fuses them through R-Vine Copula, so as to retain the respective edge distribution characteristics, accurately depict the nonlinear, asymmetric, and tail correlation, construct a vertical comprehensive drought index VCDI, significantly improve the comprehensiveness, accuracy, and robustness of drought monitoring, and have strong regional adaptability.

[0062] (2) In view of the problems that the single-threshold run theory is prone to event fragmentation and weak event misjudgment in drought event identification, the application introduces a multi-threshold system and adjacent event merging and short-time weak event elimination criteria, can effectively suppress the identification error caused by data disturbance fluctuation, ensure the continuity and reliability of the time sequence of drought events, and can simultaneously output the event duration, severity, and intensity, so as to provide high-quality event data for subsequent analysis and application of the vegetation drought resistance.

[0063] (3) The application normalizes and directionally unifies three indexes of kNDVI, LST, and VCDI with different physical meanings, constructs a three-dimensional feature space, defines a vegetation drought resistance index VDRI based on the Euclidean distance to the "extreme drought worst state point", takes into account the three-dimensional information of vegetation greenness, thermal stress, and water condition, has an explicit numerical range and an intuitive physical meaning, supports weight adaptive allocation, and can stably reflect the vegetation drought resistance under different vegetation types and regional conditions.

[0064] (4) The application uses VDRI anomaly values (A VDRI ) and standardized anomaly values (SA VDRI ) to identify the best recovery time of vegetation after a drought event, quantifies the vegetation drought recovery force by the difference between the most severe drought month and the best recovery month SA VDRI , and reveals the time lag influence of drought on vegetation function and the long-term evolution trend through lag cross-correlation analysis, so as to provide a scientific basis and time window for accurately formulating ecological restoration measures and drought adaptation management strategies.

[0065] (5) The present application accurately quantifies the drought resistance of vegetation, analyzes the influence of drought on the ecological function of vegetation, quantifies the recovery time and recovery force of vegetation after drought, and the research results provide scientific basis and data support for evaluating the stability of an ecological system. BRIEF DESCRIPTION OF DRAWINGS

[0066] In order to more clearly illustrate the technical solutions in the present application or prior art, the drawings needed to be used in the embodiments or prior art description will be briefly introduced as follows. Obviously, the drawings in the following description are some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor on the basis of these drawings.

[0067] Figure 1 A flow chart of a vegetation drought resistance index construction method according to an embodiment of the present application.

[0068] Figure 2 A spatial distribution map of kNDVI in a study area according to an embodiment of the present application.

[0069] Figure 3 A spatial distribution map of VCDI according to an embodiment of the present application.

[0070] Figure 4 A spatial distribution map of drought intensity in a study area according to an embodiment of the present application.

[0071] Figure 5 A spatial distribution map of VDRI in a study area according to an embodiment of the present application.

[0072] Figure 6 A spatial distribution map of vegetation drought recovery force according to an embodiment of the present application.

[0073] Figure 7 A spatial distribution map of the lag time of the influence of drought on the physiological function of vegetation in a study area according to an embodiment of the present application. DETAILED DESCRIPTION

[0074] In order to make the objects, technical solutions and advantages of the present application clearer, the technical solutions in the present application will be described clearly and completely in combination with the drawings and embodiments in the present application. Obviously, the described embodiments are some embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the present application.

[0075] Embodiment 1

[0076] As Figure 1As shown, the embodiment of the present application provides a vegetation drought resistance index construction method, which comprises the following steps:

[0077] Step 1, the present application selects a sub-basin as a research area, obtains multi-source time series data of the research area and unifies to the same spatial and temporal resolution, and the data adopted includes precipitation (Precipitation, P), runoff (Runoff, R), soil moisture (Soil moisture, SM), terrestrial water storage anomalies (Terrestrial water storage anomalies, TWSA) based on GRACE satellite, GLDAS land water element, kernel normalized difference vegetation index (kernel normalized difference vegetation index, kNDVI), land surface temperature (Land surface temperature, LST) and the like, and re-projection and resampling are performed.

[0078] Step 2, the Gamma distribution is adopted to fit the precipitation sequence and standardized to obtain the standardized precipitation index (Standardized precipitation index, SPI) of multiple time scales, and the precipitation comprehensive index (xSPI) is obtained based on principal component analysis (principal components analysis, PCA) fusion; the standardized runoff index (Standardized runoff index, SRI) of multiple time scales is calculated by the same method for the runoff sequence, and the runoff comprehensive index (xSRI) is obtained by PCA fusion; the monthly soil moisture is constructed by exponential smoothing to obtain the soil water index (Soil water index, SWI); the changes of soil, snow, vegetation canopy and surface water are removed by GRACE-TWSA to obtain the groundwater storage anomalies (Groundwater storage anomalies GWSA), and the GRACE groundwater drought index (GRACE groundwater drought index, GGDI) is calculated; the kernel normalized difference vegetation index (kNDVI) is constructed by using the Gaussian kernel function.

[0079] The Gamma distribution is adopted to fit the precipitation sequence and standardized to obtain the SPI of time scales of 1 month, 3 months, 6 months and 12 months, the xSPI is obtained based on principal component analysis (PCA) fusion of the SPI of 4 time scales, and the calculation formula of the SPI is as follows:

[0080] (1)

[0081] (2)

[0082] wherein denotes the probability density function (PDF) of denotes the gamma function of and denote the shape and scale parameters, respectively, and are both greater than 0. denotes the precipitation amount in the precipitation time series data. denotes the cumulative distribution function (CDF) of the PDF fitted to the time series precipitation data. denotes the inverse cumulative distribution function of the standard normal distribution. Furthermore, four SPI datasets with different time scales (SPI-1, SPI-3, SPI-6, and SPI-12) are integrated by principal component analysis (PCA) to effectively represent short-term and long-term precipitation anomalies and water shortage conditions, and the output results are recorded as xSPI.

[0083] The Gamma distribution is used to fit the precipitation series and standardized to obtain SRI with time scales of 1 month, 3 months, 6 months, and 12 months. xSRI is obtained by fusing the SRI of the four time scales based on principal component analysis (PCA). The calculation formula of SRI is as follows:

[0084] (3)

[0085] denotes the cumulative distribution function (CDF) of the PDF fitted to the time series precipitation data. denotes the inverse cumulative distribution function of the standard normal distribution. Furthermore, four SRI datasets with different time scales (SPI-1, SPI-3, SPI-6, and SPI-12) are integrated by principal component analysis (PCA) to effectively represent short-term and long-term runoff anomaly conditions, and the output results are recorded as xSRI.

[0086] The monthly soil moisture is used to construct the soil moisture index (SWI) by exponential smoothing, and the calculation formula is as follows:

[0087] (4)

[0088] wherein denotes the soil water content index of the th month. denotes the soil water content of the th month. denotes the current month when calculating the soil water content index (SWI), and denotes the soil water content of a certain past month. The time lag constant is used to reflect the time scale, and i is an integer index used to uniquely identify each historical month involved in the calculation of the soil moisture index for the current month.

[0089] By removing soil, snow, canopy, and surface water variations using TWSA retrieved from the GRACE satellite, groundwater storage anomalies (GWSA), soil moisture storage anomalies (SMSA), snow water equivalent anomalies (SWEA), canopy water storage anomalies (CWSA), and surface water storage anomalies (SWSA) were obtained from the Global Land Data Assimilation System (GLDAS) land surface model. The calculation formulas are as follows:

[0090] (5)

[0091] (6)

[0092] (7)

[0093] In the formula, Indicates the first Year Monthly GWSA value. This represents the average GWSA for month j across all years. Indicates the first Year The deviation in groundwater storage in a given month. and They represent the first Year The mean and standard deviation of monthly GSD. Indicates the first Year The GRACE groundwater drought index for the month.

[0094] The kernel-normalized difference vegetation index (kNDVI) is constructed using a Gaussian kernel function, and the specific calculation formula is as follows:

[0095] (8)

[0096] In the formula, This indicates the index's sensitivity to areas with sparse or dense vegetation. This represents near-infrared band data from remote sensing images. This represents the red band data of the remote sensing image. When When kNDVI is calculated, the formula can be simplified to the following form:

[0097] (9)

[0098] Taking kNDVI as an example, the spatial distribution map of kNDVI in the study area is as follows:Figure 2 The maximum value of kNDVI in July 2022 is 0.75, and the minimum value is 0. It can be seen from the figure that the low value is mainly concentrated in the north of the region, and the greenness of the vegetation is low. The kNDVI in the west of the region is higher than that in the east, which shows that the greenness of the local vegetation is higher.

[0099] Step 3, construction of vertical composite drought index (VCDI): calculate the marginal distribution of xSPI, xSRI, SWI and GGDI, use Copula model to fit and calculate VCDI for 4 sets of time series grid data, evaluate the goodness of fit of different structure of Copula fitting results according to RMSE, AIC and BIC, select the optimal model to obtain monthly scale VCDI.

[0100] The Copula model is used to fit and calculate VCDI for 4 sets of time series grid data. For random variables and their joint cumulative distribution function , there is a Copula function satisfying:

[0101] (10)

[0102] Where n is the number of variables. The goodness of fit (GOF) of the Copula model is evaluated by root mean square error (RMSE), Akaike information criterion (AIC) and Bayesian information criterion (BIC), and the formula is as follows:

[0103] (11)

[0104] (12)

[0105] (13)

[0106] In the formula, is the maximum likelihood value of the Copula model, is the number of parameters in the model, is the sample size used to fit the Copula model; and are the fitted value and the fitted expected value, respectively. After comparing multiple Copula models (Gaussian, Cayton, Frank, Gumbel and R-Vine), the R-Vine Copula model with the highest goodness of fit is used.

[0107] The spatial distribution map of VCDI in the study area is as follows Figure 3The figure can reflect the area affected by drought in July 2022. The maximum value of VCDI in this area is 2.5, and the minimum value is -0.3. It can be seen that the degree of drought in this area increases from south to north, especially in the northwest part of the region, where most of the VCDI values are negative, indicating that the area is more severely affected by drought.

[0108] Step 4: Identify drought events based on the multi-threshold run theory (MTRT): First, set the drought identification threshold and event identification rules, merge adjacent small events, and eliminate short and weak events to obtain the drought event set. Specifically, set the drought threshold according to the value of VCDI, and set the threshold 、 and as 0, -1 and -2, respectively. If the VCDI of the t-th month satisfies , it is recorded as a potential drought event. For a drought event with a duration of 1 month, if , it is considered as an independent drought event and retained, otherwise it is not counted. For adjacent drought events with an interval of 1 month, if , they will be merged into one drought event. Then calculate the number of drought events, event duration, severity and intensity, and the specific calculation formulas are as follows:

[0109] (14)

[0110] (15)

[0111] (16)

[0112] In the formula, represents the total duration of the identified drought events, represents the total number of drought events. represents the duration of the drought event , and represent the start time and end time indices of the drought event , respectively. represents the sum of the severity of all drought events in this period. represents the severity of the drought event . is an indicator function that takes the value 1 when , otherwise it takes the value 0. represents the drought intensity in this period; according to the identified regional drought characteristics, it can be applied to further explore the causes of changes in vegetation drought resistance.

[0113] The spatial distribution of drought intensity in the study area is shown inFigure 4 The maximum value of drought intensity in the region is 1.99 and the minimum value is 1.32. It can be seen from the figure that the regions with higher drought intensity are mainly distributed in the middle of the sub-basin. The regions with lower drought intensity are mainly distributed in the north and part of the southeast, and the spatial distribution is different from VCDI.

[0114] Step 5, Vegetation Drought Resistance Index (VDRI) construction: According to the drought events identified in step 4, the vegetation drought resistance index in each drought event time range is further calculated. Specifically, first, kNDVI, LST, and VCDI are normalized in the same direction, and the distance from the sample to the "weakest condition of vegetation drought resistance" is defined as VDRI based on the Euclidean distance of the three-dimensional feature space. Normalize kNDVI, LST, and VCDI in the Dongting Lake Basin, the formula is as follows:

[0115] (17)

[0116] (18)

[0117] (19)

[0118] In the formula, , and represent the normalized monthly kNDVI, LST, and VCDI, respectively. The subscripts max and min represent the maximum value and the minimum value, respectively. According to its physical meaning, it can be inferred that lower kNDVI and VCDI and higher LST indicate that vegetation is under greater drought stress. On the contrary, higher kNDVI and VCDI and lower LST indicate that vegetation has higher resistance under drought conditions. By establishing a three-dimensional feature space of , and , the Euclidean distance is used to calculate VDRI, and the calculation formula is as follows:

[0119] (20)

[0120] In the formula, the maximum value of , and in the three-dimensional feature space represents the minimum value of , the minimum value of , and the minimum value of , respectively. The value range of VDRI is 0 to , the larger the VDRI, the stronger the vegetation drought resistance, and the smaller the VDRI, the weaker the vegetation drought resistance.

[0121] The spatial distribution map of VDRI in the study area is as follows: Figure 5 As shown, the VDRI accurately represents the spatial distribution of vegetation drought resistance in this sub-basin in July 2022. The VDRI is lower in the northeastern part of the region, indicating lower drought resistance compared to the vegetation in the southern part.

[0122] Step 6, calculate VDRI anomalies (A VDRI To reflect changes in vegetation drought resistance, standardized outliers (SA) are used. VDRI To determine the time required for vegetation drought resistance to recover to its optimal state, the average time between the month of most severe drought and the month of optimal recovery is calculated using the SA (Self-Rating) curve. VDRI The difference is used to quantify the drought resilience of vegetation, and the formula is as follows:

[0123] (twenty one)

[0124] (twenty two)

[0125] (twenty three)

[0126] In the formula, This indicates the resilience of vegetation after drought. Indicates the first A month VDRI , Indicates the first SA of one month VDRI , Indicates the first VDRI value over a month. and These represent the mean and standard deviation of the VDRI over the selected data time range, respectively. and SA represents the months with the most severe drought and the months with the best recovery, respectively. VDRI value.

[0127] Spatial distribution map of vegetation drought resilience in the study area is shown below. Figure 6 As shown in the figure, the maximum vegetation drought resilience in this region from 2000 to 2022 was 1.21, and the minimum was 0.51. The figure shows that areas with lower values ​​are mainly distributed in the northeastern part of the sub-basin, while areas with higher vegetation drought resilience are mainly distributed in the southern part of the sub-basin.

[0128] Step 7: Use lag cross-correlation to identify the lag relationship between VDRI and vegetation ecological function indicators, and calculate the lag response months of vegetation ecological function indicators to VDRI, as shown in the following formula:

[0129] (24)

[0130] (25)

[0131] (26)

[0132] wherein, represents the maximum lag month, and in the present application, the value is set to 12. represents the maximum lag correlation between the VDRI and the vegetation ecological function parameter after the kth month, and the greater the absolute value of the value, the higher the correlation between the VDRI and the vegetation ecological function parameter. represents the covariance between and and are the standard deviations of the two sequences respectively. and respectively represent the lag correlation coefficient and the lag time.

[0133] The spatial distribution map of the lag time of the influence of drought on the physiological function of vegetation in the study area is shown in FIG. 1. Figure 7 The figure can reflect the lag month of the influence of drought on vegetation in the region. Among them, the lag month in the west of the region is smaller, indicating that the ecological function of the vegetation in the region is immediately affected under the condition of drought. The spatial distribution of the lag time of the influence of the vegetation in the east of the region is not uniform, which may be due to the agricultural irrigation and other activities, which may cause the difference in the response of the vegetation to drought.

[0134] In specific implementation, the method proposed in the technical scheme of the present application can be automatically run by a computer software technology, and the system device of the method, such as a computer readable storage medium storing the corresponding computer program of the technical scheme of the present application and a computer device including the running of the corresponding computer program, should also be within the protection scope of the present application.

[0135] Embodiment 2

[0136] Based on the same inventive concept, the present application also provides a vegetation drought resistance index construction system, which comprises a processor and a memory, the memory is used for storing program instructions, and the processor is used for calling the program instructions in the memory to execute the vegetation drought resistance index construction method as described above.

[0137] Embodiment 3

[0138] Based on the same inventive concept, the present application also provides a vegetation drought resistance index construction system, which comprises a readable storage medium, and the readable storage medium stores a computer program, and the computer program is executed to realize the vegetation drought resistance index construction method as described above.

[0139] The specific embodiments described herein are merely illustrative of the spirit of the application. Various modifications or changes in the specific embodiments described herein can occur to those skilled in the art to which the application pertains without departing from the spirit of the application, and it is understood that such modifications or changes are to be considered as within the scope of the application as defined by the appended claims.

Claims

1. A method for constructing a vegetation drought resistance index, characterized by, Comprising the following steps: Step 1, obtaining multi-source time series data of the study area and unifying to the same spatial and temporal resolution; Step 2, constructing precipitation composite index xSPI, runoff composite index xSRI, soil moisture index SWI, groundwater drought index GGDI and normalized difference vegetation index kNDVI; Step 3, calculating the edge distribution of xSPI, xSRI, SWI and GGDI, fitting the time series grid data by Copula model to calculate the vertical composite drought index, evaluating the goodness of fit of the Copula fitting results of different structures, and selecting the optimal model to obtain the monthly scale vertical composite drought index VCDI; In step 3, the Copula model is used to fit the time series raster data to calculate VCDI. For random variables and their joint cumulative distribution function , there is a Copula function satisfying: (10) Wherein n refers to the number of variables; the goodness of fit GOF of the Copula model is evaluated by root mean square error RMSE, Akaike information criterion AIC and Bayesian information criterion BIC, and the formula is as follows: (11) (12) (13) wherein, represents the maximum likelihood value of the Copula model, represents the number of parameters in the model, represents the sample size used to fit the Copula model; and respectively represent the fitted value and the fitted expected value, and the model with the highest goodness of fit is adopted after comparing a plurality of Copula models. Step 4, obtaining a set of drought events, and calculating the number of drought events, event duration, severity and intensity; Step 5, according to the drought events identified in step 4, further calculating the vegetation drought resistance index in the time range of each drought event, and performing directional consistent normalization on kNDVI, land surface temperature LST and VCDI, defining vegetation drought resistance index VDRI based on three-dimensional feature space Euclidean distance, the larger the VDRI, the stronger the vegetation drought resistance, and the smaller the VDRI, the weaker the vegetation drought resistance; The calculation formula of vegetation drought resistance index is as follows: (17) (18) (19) wherein, , and denote the normalized monthly kNDVI, LST and VCDI, respectively, and the subscripts min and max denote the minimum and maximum, respectively; by establishing a three-dimensional feature space of , and , the VDRI is calculated using the Euclidean distance, and the calculation formula is as follows: (20) wherein the maximum value of , and represent the maximum value of the minimum value of and the minimum value of , respectively.

2. The method of constructing a vegetation drought resistance index according to claim 1, wherein: The multi-source time series data in step 1 includes precipitation, runoff, soil moisture SM, land water storage anomaly TWSA based on GRACE satellite and GLDAS land water element, normalized difference vegetation index kNDVI and land surface temperature LST.

3. The method for constructing a vegetation drought resistance index as described in claim 1, characterized in that: In step 2, the precipitation sequence is fitted by Gamma distribution and standardized to obtain multi-time scale standardized runoff index SRI, and the precipitation composite index xSPI is obtained by fusion based on principal component analysis PCA, and the calculation formula is as follows: (1) (2) wherein denotes the probability density function PDF of denotes the gamma function of and denote a shape parameter and a scale parameter, respectively, and are both greater than 0; denotes the precipitation amount in the precipitation time series data; denotes the cumulative distribution function CDF of the PDF fitted to the time series precipitation data; denotes the inverse cumulative distribution function of the standard normal distribution; Similarly, the multi-time scale standardized runoff index SRI is calculated for the runoff sequence, and the runoff composite index xSRI is obtained by PCA fusion, and the calculation formula is as follows: (3) Cumulative distribution function CDF representing the PDF fitted to the time series precipitation data.

4. The method for constructing a vegetation drought resistance index as described in claim 1, characterized in that: In step 2, the monthly scale soil moisture is constructed by exponential smoothing to obtain soil moisture index SWI, and the calculation formula is as follows: (4) wherein represents the soil moisture index of the i-th month, represents the soil moisture index of the i-th month, represents the soil moisture of the i-th month, represents the soil moisture of the i-th month, represents the current month when the soil moisture index SWI is calculated, and represents the soil moisture of a certain month in the past, represents a time lag constant for reflecting a time scale, i is an integer index for uniquely identifying each historical month participating in the calculation of the soil moisture index of the current month. After eliminating the changes of soil, snow, vegetation canopy and surface water, the groundwater storage anomaly GWSA is obtained, and the groundwater drought index GGDI is calculated, and the calculation formula is as follows: (5) (6) (7) In the formula, the soil moisture storage anomaly (SMSA), snow cover water equivalent anomaly (SWEA), canopy water storage anomaly (CWSA), and surface water storage anomaly (SWSA) are all derived from the Global Land Data Assimilation System (GLDAS) land surface model. Indicates the first Year The monthly groundwater storage anomaly GWSA value, This represents the average GWSA for month j across all years. Indicates the first Year Monthly groundwater storage deviation GSD; and They represent the first Year The mean and standard deviation of monthly GSD, Indicates the first Year The monthly GRACE groundwater drought index; The Gaussian kernel function is used to construct the kernel normalized difference vegetation index kNDVI, and the specific calculation formula is as follows: (8) In the formula, represents the sensitivity of the normalized difference vegetation index to sparse or dense vegetation areas, represents near-infrared band data of the remote sensing image, represents red light band data of the remote sensing image.

5. The method of constructing a vegetation drought resistance index according to claim 1, wherein: In step 4, the drought threshold is set according to the numerical value of VCDI, set threshold , and are constants, when the vertical composite drought index value of the t-th month satisfies , it is recorded as a potential drought event; for a drought event with a duration of 1 month, if it is considered as an independent drought event and retained, otherwise it is not counted; for adjacent drought events with an interval of 1 month, if they will be combined into one drought event; The specific calculation formula of the number of drought events, event duration, severity and intensity is as follows: (14) (15) (16) In the formula, This indicates the total duration of the identified drought events. This represents the total number of drought events. Indicates drought event Duration, and These represent drought events. The index of the start and end times; This represents the sum of the severity of all drought events within a certain period. Indicates drought event The severity It is an indicator function, when The value is 1 if the condition is met, and 0 otherwise. This indicates the intensity of drought during a certain period.

6. The method of constructing a vegetation drought resistance index according to claim 1, wherein: Also included is the calculation of the VDRI anomaly A VDRI to reflect the change in vegetation drought resistance, using the VDRI standardized anomaly SA VDRI to determine the time required for vegetation drought resistance to return to the optimal state by the difference between the most severe drought month and the optimal recovery month SA VDRI of vegetation drought recovery, the formula is as follows: (21) (22) (23) wherein represents the resilience of vegetation after experiencing drought, represents the first month A VDRI , represents the first month SA VDRI , represents the first month VDRI value, and respectively represent the average and standard deviation of VDRI in the time range of selected data, and respectively represent the SA VDRI of the month when drought is most severe and the month when recovery is best.

7. The method of constructing a vegetation drought resistance index according to claim 1, wherein: It also includes identifying the lag relationship of VDRI and vegetation ecological function index by using lag cross-correlation, calculating the lag response month of vegetation ecological function index to VDRI, and the formula is as follows: (24) (25) (26) In the formula, represents the maximum lag month, represents that there is a maximum lag correlation between VDRI and the ecological function parameter of vegetation after the kth month, and the greater the absolute value of the numerical value, the higher the correlation between VDRI and the ecological function parameter of vegetation, represents the covariance between , and are the standard deviations of the two sequences respectively, and respectively represent the lag correlation coefficient and the lag time.

8. A vegetation drought resistance index construction system characterized by comprising: a vegetation drought resistance index construction program; and a computer. It includes 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 vegetation drought resistance index construction method in any one of claims 1-7.

Images