Drought assessment method considering phenological characteristics of vegetation

Abstract

The invention discloses a drought assessment method considering phenological characteristics of vegetation, relates to the technical field of ecological engineering assessment, and solves the problem that the actual influence of drought on the key growth stage of the vegetation cannot be accurately reflected in the vegetation drought condition assessment condition. The method includes the following steps that S1, Savitzky-Golay filtering and enhanced vegetation index (EVI) are adopted to identify a vegetation growth period SOS and a vegetation growth withered and yellow period EOS, and then a growth season time window of each pixel is determined; and S2, extracting a vegetation condition index VCI, a temperature condition index TCI and a vegetation health index VHI in the growing season pixel by pixel in combination with enhanced vegetation index EVI data and land surface temperature LST data, calculating a mean value, and constructing a growing season health index GSHI. According to the method, vegetation phenology is brought into drought monitoring, the influence of drought on vegetation can be evaluated more comprehensively and dynamically, and the method has an important value for improving scientificity and accuracy of regional agricultural risk management and ecological protection.

Description

Technical Field

[0001] The invention relates to the technical field of ecological engineering assessment, and in particular to a drought assessment method taking into account vegetation phenological characteristics. Background Art

[0002] Traditional drought monitoring mainly relies on meteorological observation data from discrete sites to construct indicators such as the Standardized Precipitation Index (SPI) and the Palmer Drought Index (PDSI). These indicators can reflect regional drought conditions to a certain extent, but they are limited by the coverage of observation sites and the duration of continuous observation. With the development of remote sensing technology, satellite remote sensing has become a powerful tool for vegetation drought monitoring and assessment by providing spatially clear and temporally continuous vegetation conditions and temperature stress information with its advantages of large range, multi-phase, and low cost. A variety of drought indices have been constructed, such as the Vegetation Condition Index (VCI) and the Temperature Condition Index (TCI), which are widely used to characterize drought conditions at different spatial scales and have shown effectiveness in capturing the onset, severity, and scope of drought events. Among them, the Vegetation Health Index (VHI), which integrates vegetation conditions and temperature information, has received more attention.

[0003] LP is the Loess Plateau, and LP vegetation is the vegetation of the Loess Plateau. Previous researchers used meteorological and remote sensing data to analyze the spatiotemporal characteristics of LP drought. Based on meteorological factors such as precipitation and temperature, Gao et al. used the Standardized Precipitation Evapotranspiration Index (SPEI) to explore the spatiotemporal evolution characteristics of LP; Wang et al. discussed the correlation between LP drought and atmospheric circulation based on SPI and SPEI; Liu et al. explored the changing trend and drought intensity of LP drought revealed by SPEI and SPI at different time scales. In recent years, remote sensing data has been widely used in drought monitoring in this area. Li et al. used the Vegetation Condition Index (VCI) to evaluate the drought conditions of LP and confirmed the applicability of the index in this area; Wang et al. analyzed the interannual and seasonal dynamics of LP drought based on the Temperature Vegetation Drought Index (TVDI). These studies reveal the spatiotemporal differentiation of LP drought, indicating that meteorological and remote sensing drought indices are suitable for monitoring the dynamic changes of LP drought.

[0004] However, there are still some shortcomings. Existing studies usually evaluate vegetation drought conditions based on calendar years or empirically defined growing seasons (such as April-October), which may not accurately reflect the actual impact of drought on key vegetation growth stages. LP has a vast territory, complex and diverse topography and vegetation, and significant climate differences, resulting in obvious spatial differentiation of vegetation phenology. The key periods in the development of different vegetation types and regional phenology are ignored, which may lead to estimation biases in the impact of drought on vegetation growth stages. In addition, the spatiotemporal differentiation of vegetation drought depends not only on climate fluctuations, but also on factors such as topography, hydrology, soil, and human activities. However, the current understanding of the formation mechanism of the spatiotemporal pattern of vegetation drought in the LP region is still relatively weak; therefore, it does not meet existing needs. In this regard, we propose a drought assessment method that considers vegetation phenological characteristics. Summary of the invention

[0005] The purpose of the present invention is to provide a drought assessment method taking into account vegetation phenological characteristics, so as to solve the problem that the existing vegetation drought condition assessment mentioned in the above background technology may not accurately reflect the actual impact of drought on the key growth stage of vegetation.

[0006] To achieve the above object, the present invention provides the following technical solution: a drought assessment method considering vegetation phenological characteristics, comprising the following steps:

[0007] S1: Savitzky-Golay filtering and enhanced vegetation index (EVI) are used to identify the vegetation growth period (SOS) and the vegetation growth period (EOS), and then determine the growing season time window of each pixel;

[0008] S2: Combine the enhanced vegetation index EVI data and the land surface temperature LST data to extract the vegetation condition index VCI, temperature condition index TCI and vegetation health index VHI during the growing season pixel by pixel and calculate the mean to construct the growing season health index GSHI;

[0009] S3: Combine the Standardized Precipitation Evapotranspiration Index (SPEI) and the Net Primary Productivity (NPP) data to calculate the NPP deviation (NPPd) to verify the validity of the constructed growing season health index (GSHI);

[0010] S4: Calculate drought frequency, quantitatively evaluate vegetation drought susceptibility of different intensities, and use a combination of Theil-Sen slope and Mann-Kendall analysis tests to analyze the long-term variation and significance level of the growing season health index (GSHI) at the pixel scale;

[0011] S5: Based on the growing season health index (GSHI), the long-term spatiotemporal characteristics of vegetation growing season drought are revealed, and the optimal geographic detector (OPGD) and the multi-scale geographically weighted regression (MGWR) model are used to quantitatively analyze the impact of meteorological, topographic, soil and human activities on growing season drought.

[0012] Preferably, the enhanced vegetation index EVI data and land surface temperature LST data are derived from MODIS / 006 / MOD13Q1 and MODIS / 006 / MOD11A1 data products of the Google Earth Engine platform, the temporal and spatial resolution of the enhanced vegetation index EVI data is 16 days and 500 meters, the enhanced vegetation index EVI data is preprocessed on the Google Earth Engine platform for quality control, cloud masking, image stitching, and study area clipping, and the enhanced vegetation index EVI data is provided with a quality control file QC, and only high-quality pixels with a quality control file QC value = 0 are screened and retained;

[0013] Pixels with cloud cover, emissivity error greater than 0.04, and LST data error greater than 3 K in the surface temperature were marked as invalid and reconstructed.

[0014] Preferably, the Savitzky-Golay filter is a data smoothing algorithm based on local polynomial least squares fitting, which performs smoothing on the original time series of the enhanced vegetation index EVI to obtain a smoothed EVI time series, and the algorithm is:

[0015]

[0016] Where Y is the original value of the Enhanced Vegetation Index EVI data, is the enhanced vegetation index EVI value after filtering, C i is the coefficient of the filter for the i-th enhanced vegetation index (EVI) value, j is the original sequence index of the enhanced vegetation index (EVI), N is the smoothing window size (2m+1), and m is the half-width of the window.

[0017] Preferably, the maximum value of the smoothed EVI time series is EVI max , the minimum value of the smoothed EVI time series is EVI min , EVI max With EVI min The difference is the EVI amplitude, and the SOS of the vegetation growth period and the EOS of the vegetation growth period are obtained by a dynamic threshold method;

[0018] The time point when the smoothed EVI time series exceeds 20% of the EVI amplitude is the SOS period of vegetation growth, and the time point when the smoothed EVI time series exceeds 80% of the EVI amplitude is the EOS period of vegetation growth.

[0019] The period between the SOS period (when vegetation grows green) and the EOS period (when vegetation grows yellow) is the vegetation growing season.

[0020] Preferably, the growing season health index GSHI is constructed as follows:

[0021] Step 1: Based on the information of the vegetation growth period (SOS) and the vegetation growth period (EOS), determine the annual vegetation growth season range;

[0022] Step 2: Using the 16-day composite product data of the Enhanced Vegetation Index (EVI) and the 8-day composite product data of the Land Surface Temperature (LST), the vegetation condition index (VCI), the temperature condition index (TCI) and the vegetation health index (VHI) are calculated for each 16-day period.

[0023] The temperature condition index TCI algorithm is:

[0024] The vegetation condition index VCI algorithm is:

[0025] The algorithm of vegetation health index VHI is: VHI = α × VCI + (1-α) × TCI;

[0026] Among them, i, j and k represent the i-th pixel, j-th month and k-th year respectively; max and min represent the multi-year maximum and minimum values ​​of the corresponding pixel respectively; α represents the weight factor, which is generally fixed to 0.5 in practical applications;

[0027] Step 3: Traverse the Julian day DOY interval corresponding to the vegetation growing season window pixel by pixel, extract the growing season VHI subset from the annual VHI time series, and calculate its average value to obtain the annual growing season health index GSHI of the grid cell. The algorithm is:

[0028]

[0029] GSDI y,j VHI represents the health index of pixel j in the growing season of year y; i represents the i-th VHI value within the growing season range of pixel j; DOY(SOS j ) and DOY(EOS j ) are the Julian days corresponding to the start and end dates of the growing season of pixel j; DOY(t i) is the Julian day corresponding to the i-th VHI value; n is the total number of VHI sequences in the growing season of pixel j.

[0030] Step 4: Refer to the vegetation condition index VCI, temperature condition index TCI and vegetation health index and combine them with the growing season health index GSHI to classify the drought level. .

[0031] Preferably, the SPEI is a meteorological drought index that comprehensively considers precipitation and potential evapotranspiration; NPP is a productivity index that characterizes vegetation growth conditions, which refers to the total carbon fixed by vegetation through photosynthesis per unit time and area minus the remaining amount after autotrophic respiration consumption. The formula for calculating NPP deviation (NPPd) is:

[0032]

[0033] NPP ij is the NPP value of the jth pixel in the i-th year, and are the multi-year mean NPP and multi-year standard deviation of pixel j, respectively;

[0034] NPPd characterizes the degree of deviation of NPP of each pixel relative to the multi-year variability, reflecting the intensity of the impact of extreme events such as drought on vegetation.

[0035] Preferably, the validity of the growing season health index GSHI constructed by the verification is verified using the Pearson correlation coefficient, and the calculation formula of the Pearson correlation coefficient is:

[0036]

[0037] where x i is the drought index of the i-th year, y i is the NPP or SPEI-12 value of the i-th year, n is the number of years, and the Pearson correlation coefficient ranges from [-1,1]. The larger the absolute value, the stronger the correlation.

[0038] Preferably, the calculation formula for drought frequency is:

[0039]

[0040] Where f represents the drought frequency (%), n represents the number of moderate and above droughts during the entire study period, and N represents the entire duration of the study period;

[0041] The Theil-Sen slope calculation formula of the growing season health index GSHI time series is:

[0042]

[0043] GSHI i and GSHI j is the GSHI value of the i-th year and the j-th year in the time series data; Median() is the median function; β is the Sen slope of the time series. When β>0, the growing season health index GSHI shows an increasing trend; when β<0, the growing season health index GSHI shows a decreasing trend;

[0044] The significance of the growing season health index GSHI trend was tested by the Mann-Kendall method when the time series length n ≥ 10, and the calculation formula of the statistic S is:

[0045]

[0046] where x i and x j are the i-th and j-th items of the time series respectively, and sgn(θ) is the sign function, defined as:

[0047]

[0048] When n ≥ 10, the test statistic Z after standardization of S can be calculated according to the following formula:

[0049]

[0050]

[0051] The value range of the statistic Z is (-∞, +∞). At a given significance level α, when |Z|>u 1-α / 2 , it indicates that there is a significant change in the study sequence at the α level.

[0052] Preferably, the working method of the optimal geographic detector OPGD comprises the following steps:

[0053] Step 1: Select the optimal parameters. The discretization of continuous data by geographic detectors is manually set, which may be affected by inaccurate discretization and subjective factors. On this basis, the optimal parameter geographic detector OPGD model combines different discretization classification methods and classification numbers of continuous variables, and the discretization parameter combination corresponding to the highest explanatory power of the variable is selected as the optimal one.

[0054] Step 2: Factor and interaction detection. Based on the selection of optimal parameters, single factor detection and two-factor interaction detection in the geographic detector are called. The algorithm is:

[0055]

[0056] Where q represents the explanatory power of the factor, and the range of q is (0,1). The closer it is to 1, the greater the explanatory power. h represents the explanatory variable or the stratification of the explanatory variable (Strata). Nh and N are the number of units in stratum h and the whole area. σ 2 and σh represent the variances of Y values ​​in layer h and the whole area respectively. SSW and SST are the sum of the variances within the layer and the total variance of the whole area respectively.

[0057] Preferably, the multi-scale geographically weighted regression MGWR model analyzes the driving mechanism of vegetation drought, and its formula is:

[0058] where y i is the dependent variable; x i j represents the explanatory variable; represents the bandwidth used to calibrate the j-th conditional relationship, and the bandwidth is b w The jth local regression coefficient of i ,v i ) is the spatial geographic location of the ith sample point; ε is the intercept and error term of the model.

[0059] Compared with the prior art, the present invention has the following beneficial effects:

[0060] 1. The spatial differentiation of drought in the growing season of the present invention is mainly affected by factors such as meteorology, topography, soil and human activities. Precipitation is the key meteorological factor that dominates the regional drought pattern. Topography and soil further affect the spatial heterogeneity of drought by regulating water redistribution and vegetation growth. The intensity of human activities has shown an upward trend in recent years and has become an important factor affecting drought. The multi-scale geographically weighted regression MGWR model reveals the spatial differentiation law and scale characteristics of the intensity of key influencing factors.

[0061] 2. The present invention incorporates vegetation phenology into drought monitoring, which can more comprehensively and dynamically assess the impact of drought on vegetation. It is of great value for improving the scientificity and accuracy of regional agricultural risk management and ecological protection. Considering the profound impact of climate change on vegetation phenology and ecological security, the development of vegetation drought monitoring methods that integrate phenological rhythms also has important application prospects on a global scale. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 It is a schematic diagram of the overall process of the present invention;

[0063] Figure 2 Schematic diagram of the process for constructing the growing season health index GSHI;

[0064] Figure 3 Schematic diagram of the interannual variation curves of various drought indices.

[0065] Figure 4Schematic diagram of the spatial distribution of the long-term correlation between each index and NPPd and SPEI;

[0066] Figure 5 Schematic diagram of the temporal and spatial distribution of GSHI changes in the LP region

[0067] Figure 6 Schematic diagram of the interactive detection results of drought influencing factors in four time sections. DETAILED DESCRIPTION

[0068] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0069] See also Figures 1 to 6 , an embodiment provided by the present invention: a drought assessment method considering vegetation phenological characteristics, comprising the following steps:

[0070] S1: Savitzky-Golay filtering and enhanced vegetation index (EVI) are used to identify the vegetation growth period (SOS) and the vegetation growth period (EOS), and then determine the growing season time window of each pixel;

[0071] S2: Combine the enhanced vegetation index EVI data and the land surface temperature LST data to extract the vegetation condition index VCI, temperature condition index TCI and vegetation health index VHI during the growing season pixel by pixel and calculate the mean to construct the growing season health index GSHI;

[0072] S3: The NPP deviation NPPd is calculated by combining the standardized precipitation evapotranspiration index SPEI index and the net primary productivity NPP data to verify the validity of the constructed growing season health index GSHI. The validity of the constructed growing season health index GSHI is verified using the Pearson correlation coefficient. The calculation formula of the Pearson correlation coefficient is:

[0073]

[0074] where x i is the drought index of the i-th year, y i is the NPP or SPEI-12 value of the i-th year, n is the number of years, and the value range of the Pearson correlation coefficient is [-1,1]. The larger the absolute value, the stronger the correlation;

[0075] S4: Calculate drought frequency, quantitatively evaluate vegetation drought susceptibility of different intensities, and use a combination of Theil-Sen slope and Mann-Kendall analysis tests to analyze the long-term variation and significance level of the growing season health index (GSHI) at the pixel scale;

[0076] The calculation formula for drought frequency is:

[0077]

[0078] Where f represents the drought frequency (%), n represents the number of moderate and above droughts during the entire study period, and N represents the entire duration of the study period;

[0079] The Theil-Sen slope calculation formula of the growing season health index GSHI time series is:

[0080]

[0081] GSHI i and GSHI j is the GSHI value of the i-th year and the j-th year in the time series data; Medain() is the median function; β is the Sen slope of the time series. When β>0, the growing season health index GSHI shows an increasing trend; when β<0, the growing season health index GSHI shows a decreasing trend;

[0082] The significance of the growing season health index GSHI trend was tested by the Mann-Kendall method when the time series length n ≥ 10, and the calculation formula of the statistic S is:

[0083]

[0084] where x i and x j are the i-th and j-th items of the time series respectively, and sgn(θ) is the sign function, defined as:

[0085]

[0086] When n ≥ 10, the test statistic Z after standardization of S can be calculated according to the following formula:

[0087]

[0088]

[0089] The value range of the statistic z is (-∞, +∞). Under a given significance level α, when |Z|>u 1-α / 2 , it indicates that there is a significant change in the study sequence at the α level.

[0090] S5: Based on the growing season health index GSHI, the long-term spatiotemporal characteristics of vegetation growing season drought are revealed, and the optimal geographic detector OPGD and the multi-scale geographical weighted regression MGWR model are used to quantitatively analyze the impact of meteorological, topographic, soil and human activities on growing season drought. The above-mentioned enhanced vegetation index EVI data and land surface temperature LST data are derived from the MODIS / 006 / MOD13Q1 and MODIS / 006 / MOD11A1 data products of the Google Earth Engine platform. The spatiotemporal resolution of the enhanced vegetation index EVI data is 16 days and 500 meters. The enhanced vegetation index EVI data is pre-processed on the Google Earth Engine platform for quality control, cloud mask, image stitching, and study area cropping. The enhanced vegetation index EVI data comes with a quality control file QC, and only high-quality pixels with a quality control file QC value = 0 are selected and retained;

[0091] The cloud cover of the LST data, the emissivity error greater than 0.04, and the pixels with the LST data error greater than 3K are marked as invalid and reconstructed. The Savitzky-Golay filter is a data smoothing algorithm based on local polynomial least squares fitting. The original time series of the enhanced vegetation index EVI is smoothed to obtain the smoothed EVI time series. The algorithm is:

[0092]

[0093] Where Y is the original value of the Enhanced Vegetation Index EVI data, is the enhanced vegetation index EVI value after filtering, C i is the coefficient of the filter for the i-th enhanced vegetation index (EVI) value, j is the original sequence index of the enhanced vegetation index (EVI), N is the smoothing window size (2m+1), and m is the half-width of the window.

[0094] The maximum value of the smoothed EVI time series is EVI max , the minimum value of the smoothed EVI time series is EVI min , EVI max With EVI min The difference is the EVI amplitude. The SOS of the greening period of vegetation growth and the EOS of the withering and yellowing period of vegetation growth are obtained through the dynamic threshold method. The time point when the smoothed EVI time series exceeds 20% of the EVI amplitude is the SOS of the greening period of vegetation growth, the time point when the smoothed EVI time series exceeds 80% of the EVI amplitude is the EOS of the withering and yellowing period of vegetation growth, and the time period between the SOS of the greening period of vegetation growth and the EOS of the withering and yellowing period of vegetation growth is the vegetation growing season.

[0095] The growing season health index GSHI is constructed as follows:

[0096] Step 1: Based on the information of the vegetation growth period (SOS) and the vegetation growth period (EOS), determine the annual vegetation growth season range;

[0097] Step 2: Using the 16-day composite product data of the Enhanced Vegetation Index (EVI) and the 8-day composite product data of the Land Surface Temperature (LST), the vegetation condition index (VCI), the temperature condition index (TCI) and the vegetation health index (VHI) are calculated for each 16-day period.

[0098] The temperature condition index TCI algorithm is:

[0099] The vegetation condition index VCI algorithm is:

[0100] The vegetation health index VHI algorithm is: VHI ijk =α×VCI ijk +(1-α)×TCI ijk

[0101] Among them, i, j and k represent the i-th pixel, j-th month and k-th year respectively; max and min represent the multi-year maximum and minimum values ​​of the corresponding pixel respectively; α represents the weight factor, which is generally fixed to 0.5 in practical applications;

[0102] Step 3: Traverse the Julian day DOY interval corresponding to the vegetation growing season window pixel by pixel, extract the growing season VHI subset from the annual VHI time series, and calculate its average value to obtain the annual growing season health index GSHI of the grid cell. The algorithm is:

[0103]

[0104] GSDI y,j VHI represents the health index of pixel j in the growing season of year y; i represents the i-th VHI value within the growing season range of pixel j; DOY(SOS j ) and DOY(EOS j ) are the Julian days corresponding to the start and end dates of the growing season of pixel j; DOY(t i ) is the Julian day corresponding to the i-th VHI value; n is the total number of VHI sequences in the growing season of pixel j.

[0105] Step 4: Refer to the vegetation condition index VCI, temperature condition index TCI and vegetation health index and combine with the growing season health index GSHI to classify the drought level;

[0106] The growing season health index GSHI can be divided into five levels of drought level, and the level icons are shown in the following table (Table 1):

[0107] Table 1

[0108]

[0109] Net primary productivity (NPP) refers to the carbon fixed by vegetation through photosynthesis minus the carbon consumed by autotrophic respiration. SPEI is a meteorological drought index that comprehensively considers precipitation and potential evapotranspiration. SPEI includes SPEI-12. The formula for calculating NPP deviation (NPPd) is:

[0110]

[0111] NPP ij is the NPP value of the jth pixel in the i-th year, and are the multi-year mean NPP and multi-year standard deviation of pixel j, respectively;

[0112] NPPd characterizes the degree of deviation of NPP of each pixel relative to the multi-year variability, reflecting the intensity of the impact of extreme events such as drought on vegetation.

[0113] The working method of the optimal geographic detector OPGD includes the following steps:

[0114] Step 1: Select the optimal parameters. The discretization of continuous data by geographic detectors is manually set, which may be affected by inaccurate discretization and subjective factors. On this basis, the optimal parameter geographic detector OPGD model combines different discretization classification methods and classification numbers of continuous variables, and the discretization parameter combination corresponding to the highest explanatory power of the variable is selected as the optimal one.

[0115] Twelve indicators that may affect drought differences were selected from five dimensions of meteorology, topography, soil, human activities and transportation location. The OPGD model was used to reveal the explanatory power of different factors on the spatial differentiation of vegetation drought. Time section data from 2005, 2010, 2015 and 2020 were selected at intervals of 5 years to quantitatively evaluate the influence of the 12 factors on GSHI differentiation. The following table (Table 2) was obtained. Table 2 is a table of driving factors affecting the growing season health index GSHI, thereby screening out the dominant driving factors.

[0116] Table 2

[0117]

[0118] Step 2: Factor and interaction detection. Based on the selection of optimal parameters, single factor detection and two-factor interaction detection in the geographic detector are called. The algorithm is:

[0119]

[0120]

[0121] Where q represents the explanatory power of the factor, and the range of q is (0,1). The closer it is to 1, the greater the explanatory power. h represents the explanatory variable or the stratification of the explanatory variable (Strata). Nh and N are the number of units in stratum h and the whole area. σ 2 and σh represent the variance of Y value of layer h and the whole area respectively, SSW and SST are the sum of variance within layer and the total variance of the whole area respectively. The criteria for distinguishing different types of interaction are shown in the following table (Table 3):

[0122] Table 3

[0123]

[0124] The multi-scale geographically weighted regression (MGWR) model analyzes the driving mechanism of vegetation drought, and its formula is:

[0125] where y i is the dependent variable; x i j represents the explanatory variable; represents the bandwidth used to calibrate the j-th conditional relationship, and the bandwidth is b w The jth local regression coefficient of i ,v i ) is the spatial geographic location of the ith sample point; ε is the intercept and error term of the model.

[0126] Based on MODIS-EVI, the spatiotemporal patterns of the LP vegetation greening period (SOS) and yellowing period (EOS) from 2003 to 2022 were analyzed. The spatial differentiation of LP vegetation phenology was obvious, showing a zonal law that changes with the climate gradient.

[0127] The spatial pattern of SOS generally shows a trend of gradual postponement from southeast to northwest. The SOS in the southeast is concentrated in late March to early April, while the SOS in the northwest is postponed to late April to early May. Due to the differences in regional water and heat conditions, the spring temperature in the southeast rises earlier and the precipitation is abundant, which promotes the early greening of vegetation. However, the spring temperature in the northwest is lower and the precipitation is scarce, which inhibits the growth process of vegetation. The interannual variation of SOS is characterized by large fluctuations from 2003 to 2008, relative stability from 2009 to 2017, and an early trend after 2018. Among them, the interannual fluctuation range in the north is greater than that in the south.

[0128] EOS has different changing trends in different geographical regions. EOS in most regions is concentrated between October and November. The vegetation withers earlier in the northeastern region and later in the southern region. This may be due to the early drop in temperature and reduced precipitation in the north in autumn, which causes vegetation to enter dormancy earlier. In the south, the autumn climate conditions are relatively good. The interannual variability of EOS is more stable than SOS. Considering that the autumn climate fluctuations are relatively small, the impact on vegetation growth is limited. However, in some years such as 2008 and 2012, the EOS is generally early, which may be related to the extreme climate events in the autumn of that year.

[0129] In order to verify the effectiveness of the constructed GSHI for LP vegetation drought monitoring, NPPd and SPEI-12 were selected as reference indicators. Figure 3 The interannual variation curves of various drought indices (VCI, TCI, VHI and GSHI) and NPPd and SPEI-12 are shown. In general, the drought indices can well capture the evolution of vegetation drought in LP in the past 20 years, and the variation trends are generally consistent. For example, troughs appeared in 2005 and 2009, indicating that LP experienced more severe droughts in these two periods, which is consistent with the actual situation.

[0130] There are some differences in the changing trends of TCI, especially after 2012, when TCI showed a clear downward trend, while other indices were relatively stable, indicating that in the context of climate warming, only considering the temperature factor may not accurately reflect the health of vegetation. In contrast, the changes of GSHI and NPPd are more synchronized and show higher consistency, which shows that by comprehensively considering the vegetation growth status and temperature stress information, and combining the vegetation phenology differences at the pixel scale, GSHI can more comprehensively and accurately reveal the actual impact of drought on vegetation. To quantitatively evaluate the performance of GSHI in characterizing vegetation drought, the long-term correlation coefficients of each index with NPPd and SPEI-12 were further calculated. The long-term correlation coefficients are shown in the following table (Table 4):

[0131] Table 4

[0132]

[0133] The results showed that GSHI had the highest correlation with the two reference indicators, reaching 0.803 and 0.635, respectively, which was better than VCI, TCI and VHI;

[0134] like Figure 4The spatial distribution map shown shows that the correlation between GSHI and NPPd and SPEI-12 shows strong spatial consistency throughout the LP, and the proportion of regional areas that pass the 0.05 significance test is more than 90%. The higher regional correlation of GSHI is attributed to the comprehensive consideration of the spatiotemporal differences in vegetation phenology, which effectively captures the impact of drought in the critical growth period of vegetation, and can more comprehensively and accurately reveal the compound stress effect of drought on the vegetation growing season than a single index. The correlation between GSHI and NPPd and SPEI-12 shows strong spatial consistency throughout the LP, and the proportion of regional areas that pass the 0.05 significance test is more than 90%. The higher regional correlation of GSHI is attributed to the comprehensive consideration of the spatiotemporal differences in vegetation phenology, which effectively captures the impact of drought in the critical growth period of vegetation, and can more comprehensively and accurately reveal the compound stress effect of drought on the vegetation growing season than a single index;

[0135] like Figure 5 As shown in the figure, from the spatial and temporal distribution of GSHI, drought in the LP vegetation growing season shows obvious regional differences. In general, the degree of drought gradually intensifies from the southeast to the northwest, forming a pattern of "southeast humid-central semi-humid-northwest semi-arid". The southeast region has abundant precipitation, good heat conditions, high GSHI and low drought risk; the central region has a relatively dry climate and a slight decrease in GSHI; the northwest region belongs to the semi-arid climate zone, with generally low GSHI, and is a high-incidence area for vegetation drought. In addition, around 2005 and 2009 were high-incidence periods of drought. After 2010, the GSHI showed an overall fluctuating upward trend, indicating that the drought risk in the LP vegetation growing season has generally decreased. However, in individual years such as 2013 and 2015, severe and extreme droughts rebounded due to the occurrence of local droughts.

[0136] Four time sections, 2005, 2010, 2015 and 2020, were selected, and the OPGD model was used to analyze the explanatory power of different factors on the spatial differentiation of GSHI, as shown in the following table (Table 5):

[0137] Table 5

[0138]

[0139] From the perspective of influencing factors, meteorological factors have the strongest explanatory power for the spatial pattern of vegetation drought, among which the q value of mean annual precipitation (PRC) is the highest, averaging 0.34 in the four periods, indicating that moisture conditions are the key meteorological factors that dominate the regional drought pattern, followed by mean annual temperature (TME) and solar radiation (SRD), with q values ​​of around 0.25. Among topographic factors, elevation (ELV) and slope (SLP) have more significant effects, with q values ​​of 0.24 and 0.05, respectively. Soil sand content (SSM) and land use type (LUT) have comparable explanatory power, with q values ​​of around 0.20. In contrast, socioeconomic factors such as population density (POP), night lights (NLC) and transportation location factors such as distance to the city (DTC) and distance to the water system (DTR) have relatively weak explanatory power, with q values ​​mostly below 0.10;

[0140] To further reveal the impact of the interaction between different driving factors on the spatial differentiation of GSHI, the interaction detection function was used to analyze the explanatory power of the combination of meteorological, topographic, soil, and socio-economic factors, such as Figure 6 As shown, the interactive detection results of the main influencing factors of the four time sections of 2005, 2010, 2015 and 2020 are displayed;

[0141] Overall, the interaction q values ​​of most factor combinations are higher than the single factor q values, and there is an obvious enhancement effect between different factors. Among them, the interaction between PRC and other factors is the most significant. The combination of PRC with TME, ELV, and LUT all show nonlinear enhancement effects, and the interaction q values ​​are above 0.5. This shows that the matching of moisture and heat conditions, the coordinated distribution of precipitation and topography, and the coupling of water-heat combination with underlying surface type are key factors leading to regional drought. In addition, the explanatory power of SRD combined with temperature, precipitation, and wind speed (WND) is also high, and the q values ​​are all above 0.4, showing a double factor enhancement effect. This may be because solar radiation affects the surface energy balance, thereby regulating regional water evaporation and vegetation transpiration, and acts together with other meteorological elements on the drought process.

[0142] From the perspective of temporal dynamics, the interactions of some factor combinations vary from year to year. For example, the interactions between TME and LUT, and between PRC and ELV increased significantly between 2005 and 2015, but weakened in 2020. This is related to human interventions such as land use changes and the implementation of vegetation restoration projects in different periods. In contrast, the interactions between meteorological elements such as PRC and TME, and WND are relatively stable, indicating that the impact of climate change on regional drought patterns is long-term and continuous. It is worth noting that the interactions between socioeconomic factors such as population density (POP) and night lights (NLC) and natural factor combinations are generally weak. This may be because socioeconomic activities are mainly limited to local areas such as towns, while natural elements such as meteorology and topography dominate the drought pattern at the regional scale.

[0143] From the perspective of temporal dynamics, the explanatory power of most factors varies from year to year, but the overall trend of change is consistent. The influence of meteorological and topographic factors is relatively stable, while the role of socioeconomic factors has increased. For example, the q value of PRC remained above 0.28 in all four periods; the influence of ELV decreased slightly; while the q values ​​of POP and DTC showed an upward trend, indicating that the impact of human activities on the drought process in local areas is intensifying;

[0144] In order to further reveal the scale dependence of key influencing factors on the spatial differentiation of GSHI, a spatial regression analysis of the four main controlling factors of annual precipitation (PRC), annual solar radiation (SRD), elevation (ELV) and soil sediment content (SSM) was carried out using an 8km×8km fishing net and the MGWR model to explore the spatiotemporal variation characteristics of the influence intensity of different factors. The results are shown in the following table (Table 6, Table 6 is a schematic table for comparing the performance evaluation indicators of the GWR and MGWR models):

[0145] Table 6

[0146]

[0147] The results show that the model determination coefficient R of the MGWR model in the four periods is 2 All of them are above 0.6, and the AICc value is better than that of GWR, indicating that the model has good explanatory power, as shown in the following table (Table 7, Table 7 is the model fitting effect of MGWR and the statistical characteristics of the main influencing factors):

[0148] Table 7

[0149]

[0150]

[0151] The ranking of the mean values ​​of the regression coefficients of the factors in each period is generally consistent, from large to small, they are PRC, ELV, SRD and SSM;

[0152] As shown in the following table (Table 8, Table 8 is a schematic table of factor scale bandwidth estimation of the MGWR model):

[0153] Table 8

[0154]

[0155] From the factor scale bandwidth estimated by MGWR, the scale bandwidth of PRC and SRD is much larger than that of ELV and SSM, indicating that the influence of meteorological factors has a relatively large spatial scale, while the effects of topography and soil factors are more limited to the local scale. This means that the spatial pattern of LP vegetation drought is jointly affected by multi-scale factors, the regional climate background dominates the overall pattern, and local environmental factors regulate spatial heterogeneity. In addition, the scale bandwidth of SRD and SSM shows a decreasing trend year by year, indicating that the impact scale of these two factors has gradually become localized in recent years. This is related to the increase in local heterogeneity caused by the intensification of human activities such as vegetation cover changes and land use adjustments.

[0156] It will be apparent to those skilled in the art that the invention is not limited to the details of the exemplary embodiments described above and that the invention can be implemented in other specific forms without departing from the spirit or essential features of the invention. Therefore, the embodiments should be considered exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description, and it is intended that all variations falling within the meaning and scope of the equivalent elements of the claims be included in the invention. Any reference numeral in a claim should not be considered as limiting the claim to which it relates.

Claims

1. A drought assessment method considering vegetation phenological characteristics, characterized in that: The steps include: S1: Savitzky-Golay filtering and enhanced vegetation index (EVI) are used to identify the vegetation growth period (SOS) and the vegetation growth period (EOS), and then determine the growing season time window of each pixel; S2: Combine the enhanced vegetation index EVI data and the land surface temperature LST data to extract the vegetation condition index VCI, temperature condition index TCI and vegetation health index VHI during the growing season pixel by pixel and calculate the mean to construct the growing season health index GSHI; S3: Combine the Standardized Precipitation Evapotranspiration Index (SPEI) and the Net Primary Productivity (NPP) data to calculate the NPP deviation (NPPd) to verify the validity of the constructed growing season health index (GSHI); S4: Calculate drought frequency, quantitatively evaluate vegetation drought susceptibility of different intensities, and use a combination of Theil-Sen slope and Mann-Kendall analysis tests to analyze the long-term variation and significance level of the growing season health index (GSHI) at the pixel scale; S5: Based on the growing season health index (GSHI), the long-term spatiotemporal characteristics of vegetation growing season drought are revealed, and the optimal geographic detector (OPGD) and the multi-scale geographically weighted regression (MGWR) model are used to quantitatively analyze the impact of meteorological, topographic, soil and human activities on growing season drought.

2. A drought assessment method considering vegetation phenological characteristics according to claim 1, characterized in that: The enhanced vegetation index EVI data and land surface temperature LST data are derived from MODIS / 006 / MOD13Q1 and MODIS / 006 / MOD11A1 data products on the Google Earth Engine platform. The temporal and spatial resolutions of the enhanced vegetation index EVI data are 16 days and 500 meters. The enhanced vegetation index EVI data is preprocessed on the Google Earth Engine platform for quality control, cloud masking, image stitching, and study area clipping. The enhanced vegetation index EVI data is provided with a quality control file QC, and only high-quality pixels with a quality control file QC value of 0 are screened and retained; Pixels with cloud cover, emissivity error greater than 0.04, and LST data error greater than 3 K in the surface temperature were marked as invalid and reconstructed.

3. A drought assessment method considering vegetation phenological characteristics according to claim 1, characterized in that: The Savitzky-Golay filter is a data smoothing algorithm based on local polynomial least squares fitting. The original time series of the enhanced vegetation index EVI is smoothed to obtain the smoothed EVI time series. The algorithm is: Where Y is the original value of the Enhanced Vegetation Index EVI data, is the enhanced vegetation index EVI value after filtering, C i is the coefficient of the filter for the i-th enhanced vegetation index (EVI) value, j is the original sequence index of the enhanced vegetation index (EVI), N is the smoothing window size (2m+1), and m is the half-width of the window.

4. The drought assessment method considering vegetation phenological characteristics according to claim 1, characterized in that: The maximum value of the smoothed EVI time series is EVI max , the minimum value of the smoothed EVI time series is EVI min , EVI max With EVI min The difference is the EVI amplitude, and the SOS of the vegetation growth period and the EOS of the vegetation growth period are obtained by a dynamic threshold method; The time point when the smoothed EVI time series exceeds 20% of the EVI amplitude is the SOS period of vegetation growth, and the time point when the smoothed EVI time series exceeds 80% of the EVI amplitude is the EOS period of vegetation growth. The period between the SOS period (when vegetation grows green) and the EOS period (when vegetation grows yellow) is the vegetation growing season.

5. The drought assessment method considering vegetation phenological characteristics according to claim 1, characterized in that: The growing season health index GSHI is constructed as follows: Step 1: Based on the information of the vegetation growth period (SOS) and the vegetation growth period (EOS), determine the annual vegetation growth season range; Step 2: Using the 16-day composite product data of the Enhanced Vegetation Index (EVI) and the 8-day composite product data of the Land Surface Temperature (LST), the vegetation condition index (VCI), the temperature condition index (TCI) and the vegetation health index (VHI) are calculated for each 16-day period. The temperature condition index TCI algorithm is: The vegetation condition index VCI algorithm is: The vegetation health index VHI algorithm is: VHI ijk =α×VCI ijk +(1-α)×TCI ijk ; Among them, i, j and k represent the i-th pixel, j-th month and k-th year respectively; max and min represent the multi-year maximum and minimum values ​​of the corresponding pixel respectively; α represents the weight factor, which is generally fixed to 0.5 in practical applications; Step 3: Traverse the Julian day DOY interval corresponding to the vegetation growing season window pixel by pixel, extract the growing season VHI subset from the annual VHI time series, and calculate its average value to obtain the annual growing season health index GSHI of the grid cell. The algorithm is: Among them, GSDI y,j VHI represents the health index of pixel j in the growing season of year y; i represents the i-th VHI value within the growing season range of pixel j; DOY(SOS j ) and DOY(EOS j ) are the Julian days corresponding to the start and end dates of the growing season of pixel j; DOY(t i ) is the Julian day corresponding to the i-th VHI value; n is the total number of VHI sequences in the growing season of pixel j. Step 4: Refer to the vegetation condition index VCI, temperature condition index TCI and vegetation health index and combine them with the growing season health index GSHI to classify the drought level.

6. The drought assessment method considering vegetation phenological characteristics according to claim 1, characterized in that: The SPEI is a meteorological drought index that comprehensively considers precipitation and potential evapotranspiration; NPP is a productivity index that characterizes vegetation growth conditions, which refers to the total carbon fixed by vegetation through photosynthesis per unit time and area minus the remaining amount after autotrophic respiration consumption. The formula for calculating NPP deviation (NPPd) is: NPP ij is the NPP value of the jth pixel in the i-th year, and are the multi-year mean NPP and multi-year standard deviation of pixel j, respectively; NPPd characterizes the degree of deviation of NPP of each pixel relative to the multi-year variability, reflecting the intensity of the impact of extreme events such as drought on vegetation.

7. The drought assessment method considering vegetation phenological characteristics according to claim 1, characterized in that: The validity of the growing season health index GSHI constructed by the verification is verified using the Pearson correlation coefficient, and the calculation formula of the Pearson correlation coefficient is: where x i is the drought index of the i-th year, y i is the NPP or SPEI-12 value of the i-th year, n is the number of years, and the Pearson correlation coefficient ranges from [-1,1]. The larger the absolute value, the stronger the correlation.

8. The drought assessment method considering vegetation phenological characteristics according to claim 1, characterized in that: The calculation formula for drought frequency is: Where f represents the drought frequency (%), n represents the number of moderate and above droughts during the entire study period, and N represents the entire duration of the study period; The Theil-Sen slope calculation formula of the growing season health index GSHI time series is: GSHI i and GSHI j is the GSHI value of the i-th year and the j-th year in the time series data; Median() is the median function; β is the Sen slope of the time series. When β>0, the growing season health index GSHI shows an increasing trend; when β<0, the growing season health index GSHI shows a decreasing trend; The significance of the growing season health index GSHI trend was tested by the Mann-Kendall method when the time series length n ≥ 10, and the calculation formula of the statistic S is: where x i and x j are the i-th and j-th items of the time series respectively, and sgn(θ) is the sign function, defined as: When n ≥ 10, the test statistic Z after standardization of S can be calculated according to the following formula: The value range of the statistic Z is (-∞, +∞). Under a given significance level α, when |Z|>u 1-α / 2 , it indicates that there is a significant change in the study sequence at the α level.

9. The drought assessment method considering vegetation phenological characteristics according to claim 1, characterized in that: The working method of the optimal geographic detector OPGD comprises the following steps: Step 1: Select the optimal parameters. The discretization of continuous data by geographic detectors is manually set, which may be affected by inaccurate discretization and subjective factors. On this basis, the optimal parameter geographic detector OPGD model combines different discretization classification methods and classification numbers of continuous variables, and the discretization parameter combination corresponding to the highest explanatory power of the variable is selected as the optimal one. Step 2: Factor and interaction detection. Based on the selection of optimal parameters, single factor detection and two-factor interaction detection in the geographic detector are called. The algorithm is: Where q represents the explanatory power of the factor, and the range of q is (0,1). The closer it is to 1, the greater the explanatory power. h represents the explanatory variable or the stratification of the explanatory variable (Strata). Nh and N are the number of units in stratum h and the whole area. σ 2 and σh represent the variances of Y values ​​in layer h and the whole area respectively. SSW and SST are the sum of the variances within the layer and the total variance of the whole area respectively.

10. The drought assessment method considering vegetation phenological characteristics according to claim 1, characterized in that: The multi-scale geographically weighted regression (MGWR) model analyzes the driving mechanism of vegetation drought, and its formula is for: where y i is the dependent variable; x ij represents the explanatory variable; represents the bandwidth used to calibrate the j-th conditional relationship, and the bandwidth is b w The jth local regression coefficient of i ,v i ) is the spatial geographic location of the i-th sample point; ε is the intercept and error term of the model.

Images