A quantitative analysis method for the impact of human activities on vegetation cover changes

Abstract

The invention discloses a quantitative analysis method for the change of vegetation coverage caused by human activities, which includes the following steps: S1, using MODIS / NDVI remote sensing data to characterize vegetation coverage and VIIRS / DNB remote sensing data to characterize human activity intensity, respectively constructing two time series, MTS and VTS, through time series preprocessing technology; S2, constructing time series segmentation technology to realize iterative segmentation of time series; S3, constructing time series merging and feature extraction technology, iteratively realizing the merging of time series through sorting angle method, and extracting time series features; S4, calculating MTS and VTS and analyzing spatial pattern by analyzing the time series features of MTS and VTS, and realizing quantitative analysis of the correlation between MTS and VTS. By integrating time series segmentation, time series merging, spatial analysis and statistics and other technologies, quantitative analysis of the impact of human activities on vegetation coverage is realized.

Description

Technical Field

[0001] The present invention relates to the field of geospatial analysis and application technology, and in particular to a quantitative analysis method for the impact of human activities on vegetation cover changes. Background Art

[0002] In recent years, the intensification of human activities has greatly increased the demand for construction land, thereby reducing vegetation coverage, such as cultivated land, forest land, grassland, etc. However, some studies have shown that vegetation coverage is gradually increasing, especially in some urban areas. Mainly due to the improvement of ecological environmental protection awareness and technology, the vegetation coverage in urban areas has increased year by year. Therefore, there is still a lack of corresponding technical methods for quantitative analysis of changes in vegetation coverage caused by human activities.

[0003] Changes in vegetation cover and human activities are a long-term dynamic process. Remote sensing data can provide relevant monitoring data for large areas and long time series. Nighttime Light (NTL) remote sensing images reflect the intensity of human nighttime activities and can represent the change information of human activities through unified numerical values.

[0004] The time series analysis method has achieved good results in vegetation coverage in large areas, such as abandoned farmland and disturbed forest land. At the same time, this method has also achieved good results in monitoring and analyzing changes in economic and social indicators. There are already a variety of mature time series analysis methods for the interannual change analysis of NDVI and NTL data, such as OLS model, STMAP model, TS model, etc. However, for the changes in human activities and vegetation coverage, there are large differences in the time when they start to change within a certain research period, and there are also different directions of change. At the same time, the changes in the time series are not necessarily gradual and continuous, and there are mutations. Therefore, time series segmentation is a necessary task to identify the temporal change pattern. Secondly, the impact of human activities on geographical and environmental processes often has a lag effect, especially the impact on vegetation cover changes. Therefore, time series segmentation is also a prerequisite for analyzing the time lag effect and finding the core change period.

[0005] Existing time series segmentation algorithms include LandTrendr (Landsat-based detection of Trends in Disturbance and Recovery algorithm), MODTrendr (MODIS-based detection of Trends in Disturbance and Recovery) algorithm, BFAST (Breaks For Additive and Seasonal Trend) algorithm, DBEST (Detecting Breakpoints and Estimating Segments in Trend) algorithm, etc. These algorithms fit the time series curve by setting different rules and thresholds, and find the inflection points of the curve change to achieve the purpose of preliminary segmentation of the original time series. Then, through iterative fitting and inflection point search, the time series curve will be segmented more accurately. However, these time series segmentation and fitting methods require a large number of pre-set parameters and thresholds. For large-scale time series analysis, there is still a lack of simplified, efficient and accurate time series analysis methods to achieve quantitative analysis of human activities on vegetation cover changes.

[0006] Therefore, the current quantitative analysis of the impact of human activities on vegetation cover still has certain technical defects:

[0007] (1) Efficient and accurate time series analysis method for remote sensing data.

[0008] (2) Analysis of time series of vegetation cover and human activities over large areas and long time periods.

[0009] (3) Quantitative analysis indicators of the impact of human activities on vegetation cover changes. Summary of the invention

[0010] In view of the problems in the related technology, the present invention proposes a quantitative analysis method for the impact of human activities on vegetation cover changes, so as to overcome the above-mentioned technical problems existing in the existing related technology.

[0011] To this end, the specific technical solution adopted by the present invention is as follows: the method comprises the following steps:

[0012] S1. Use MODIS / NDVI remote sensing data to represent vegetation cover and VIIRS / DNB remote sensing data to represent human activity intensity, and construct two time series, MTS and VTS, respectively, through time series preprocessing technology;

[0013] S2. Build time series segmentation technology to achieve iterative segmentation of time series;

[0014] S3. Build time series merging and feature extraction technology, iterate through sorting angle method to merge time series and extract time series features;

[0015] S4. By analyzing the time series characteristics of MTS and VTS, the calculation and spatial pattern analysis of MTS and VTS are carried out to achieve quantitative analysis of the correlation between MTS and VTS.

[0016] Furthermore, the two time series, MTS and VTS, are constructed respectively by using the time series preprocessing technology, including the following steps:

[0017] S11, construct VTS based on cf_cvg quality control file and perform smoothing;

[0018] S12, smoothing the MTS curve based on the Reliability layer;

[0019] S13, further smoothing the MTS and VTS through Savitzky-Golay filter;

[0020] S14, determine the threshold of the 6-neighborhood standard deviation by randomly selecting some samples and remove high-value outliers in MTS and VTS;

[0021] S15. Obtain the standard deviation threshold by calculating the standard deviation of random sample points and manually visually interpreting them.

[0022] Furthermore, the calculation expression for determining the threshold of the 6-neighborhood standard deviation by randomly selecting some samples includes:

[0023] n∈(min(0,n-3),max(n+3,73 or 140))

[0024] |x m -μ|<=5σ,if|x i -μ|>5σ

[0025] Where n represents the set of 6 neighborhoods in VTS or MTS;

[0026] x n Represents the value of each element of the 6-neighborhood set;

[0027] μ represents the average value of the 6-neighborhood set;

[0028] σ represents the standard deviation of the 6-neighborhood set;

[0029] m represents the set of valid elements in the 6-neighborhood;

[0030] card represents the effective length of set n or set m;

[0031] x m Represents the value of each element of the valid value set m, that is, |x m -μ|<=5σ;

[0032] x i represents an outlier in the time series, i.e. |x i -μ|>5σ;

[0033] Indicates an outlier x i Correction value of

[0034] Furthermore, the time series segmentation technology is constructed to achieve iterative segmentation of the time series, including the following steps:

[0035] S21, fitting the MTS and VTS curves using the least squares first-order linear regression algorithm, and calculating the vertical distance between the time series data points and the fitting curve;

[0036] S22, setting an inflection point according to the vertical distance between the fitting curve and the time series data point;

[0037] S23, dividing the current time series into two sub-time series by the inflection point, wherein the inflection point becomes the end point of the first sub-time series and the start point of the second sub-time series;

[0038] S24. Use the principle of minimizing the BIC index to determine the optimal number of time series segmentation times;

[0039] S25. Repeat the above segmentation steps to iteratively segment the MTS and VTS.

[0040] Furthermore, setting the inflection point according to the vertical distance between the fitting curve and the time series data point comprises the following steps:

[0041] S221, setting the time series data point with the largest vertical distance from the fitting curve as the inflection point;

[0042] S222: When the length of the segmented time series corresponding to the inflection point is less than 3 months, that is, the MTS is greater than 6 data points and the VTS is greater than 3 data points, the time series data point with the second largest vertical distance to the fitting curve is taken as the inflection point.

[0043] Furthermore, the time series merging and feature extraction technology is constructed, which realizes the merging of time series through the sorting angle method iteration and extracts the time series features, including the following steps:

[0044] S31, using the sorting angle method, select two sub-time series with the smallest angle between adjacent fitting lines and the same change direction, and merge them preferentially;

[0045] S32, for the newly merged time series, use the least squares first-order linear regression algorithm to fit the curve and calculate its overall trend, until any two adjacent sub-time series do not meet the requirements for merging;

[0046] S33, when the merging process is completed, new vertices and segmented sub-time series are formed, and the new vertices are fitted to obtain overall trend characteristics;

[0047] S34, for a non-head-tail vertex, the end point of the previous sub-time series and the start point of the next sub-time series may exist simultaneously on a date, and the specific value of the non-head-tail vertex is determined by calculating the average of the end point and the start point;

[0048] S35. After the segmentation and merging are completed, the descriptive features of the fitting curve are defined and obtained.

[0049] Furthermore, the descriptive features include core change period, maximum trend feature, change duration feature and change start time;

[0050] The core change period represents a sub-time series that is consistent with the overall trend change direction and has the largest absolute change value;

[0051] The maximum trend feature represents the absolute value of the change of the fitting curve within the core change period;

[0052] The change duration feature represents the change time of the core change period;

[0053] The change start time indicates the date when the core change period starts.

[0054] Furthermore, the calculation and spatial pattern analysis of MTS and VTS are performed by analyzing the time series characteristics of MTS and VTS to achieve quantitative analysis of the correlation between MTS and VTS, including the following steps:

[0055] S41. By analyzing the overall change trend characteristics of VTS, the areas with large changes in human activities are marked as V-change areas, and these areas are used as the target areas for analysis;

[0056] S42, using a Gaussian fitting method to determine the threshold of the VTS change trend;

[0057] S43, performing a time series analysis of the vegetation coverage MTS in the V-change area;

[0058] S44, calculate the time series analysis results and characteristics of MTS and VTS;

[0059] S45, randomly selecting certain sample points manually, and marking the time series related features of the sample points;

[0060] S46, analyze the spatiotemporal pattern of the obtained time series characteristics and results through the global Moran index;

[0061] S47. Calculate the Moran index of a single variable based on time series characteristics such as overall trend, maximum trend, time lag and duration;

[0062] S48. The obtained maximum change trend characteristics of MTS and VTS are classified, and the VII_MOD_level index is constructed to integrate the two types of characteristics, so as to intuitively display the relationship between MTS and VTS.

[0063] Furthermore, the calculation expressions for calculating the time series analysis results and characteristics of MTS and VTS include:

[0064] Lag_effect = Start_date MODIS -Start_date VIIRS

[0065] Duration_diff = Duration MODIS -Duration VIIRS

[0066] Where, Lag_effect represents the time lag feature;

[0067] Duration_diff represents the duration difference feature;

[0068] Start_date MODIS Indicates the starting time characteristics of the time series change of MTS;

[0069] Start_date VIIRS Indicates the starting time characteristics of the time series change of VTS;

[0070] Duration MODIS Represents the duration characteristics of the time series changes of MTS;

[0071] Duration VIIRS Represents the time series variation duration characteristics of VTS.

[0072] Furthermore, the obtained maximum change trend characteristics of MTS and VTS are classified, and the calculation expression of the comprehensive calculation of the two types of characteristics by constructing the VII_MOD_level indicator includes:

[0073]

[0074] VII_MOD_level=MOD_index×10+VIIRS_index

[0075] In the formula, MOD_index represents the classification index of the maximum change trend of MTS;

[0076] VIIRS_index represents the grading index of the maximum change trend of VTS;

[0077] VII_MOD_level represents the comprehensive grading index;

[0078] high indicates high classification;

[0079] medium means medium grade;

[0080] Low means low grade.

[0081] The beneficial effects of the present invention are: by integrating time series segmentation, time series merging, spatial analysis and statistics, quantitative analysis of the impact of human activities on vegetation cover is achieved. First, Southeast Asia is used as the study area, MODIS / NDVI remote sensing data are used to characterize vegetation cover, and VIIRS / DNB remote sensing data are used to characterize human activity intensity, and MTS and VTS long time series are constructed respectively. Secondly, time series segmentation, time series merging and other technologies are used to extract time series features such as overall trend, maximum trend, change start time, and duration. Finally, the time lag and duration difference of VTS and MTS are calculated, and their spatial patterns are analyzed to achieve quantitative analysis of human activities on vegetation cover changes, and the time series features of manual interpretation are analyzed to verify the effectiveness of the method. . BRIEF DESCRIPTION OF THE DRAWINGS

[0082] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0083] Figure 1 is a flow chart of a method for quantitatively analyzing changes in vegetation coverage due to human activities according to an embodiment of the present invention;

[0084] Figure 2 is a workflow diagram of a method for quantitatively analyzing changes in vegetation coverage due to human activities according to an embodiment of the present invention;

[0085] Figure 3 It is a schematic diagram of time series segmentation and merging in a quantitative analysis method of human activities on vegetation cover changes according to an embodiment of the present invention;

[0086] Figure 4 It is a schematic diagram of the VTS overall trend grid ratio and Gaussian fitting curve in a quantitative analysis method of human activities on vegetation cover changes according to an embodiment of the present invention;

[0087] Figure 5 It is a schematic diagram of spatial distribution of VTS overall trend in a quantitative analysis method of human activities on vegetation cover changes according to an embodiment of the present invention;

[0088] Figure 6 It is a schematic diagram of the extraction results of Chinese residential grids and regional display in a quantitative analysis method of the effect of human activities on vegetation cover changes according to an embodiment of the present invention;

[0089] Figure 7 It is a schematic diagram of the time lag analysis results of VTS and MTS in a quantitative analysis method of human activities on vegetation cover changes according to an embodiment of the present invention;

[0090] Figure 8 It is a schematic diagram of the difference analysis results of VTS and MTS duration in a quantitative analysis method of human activities on vegetation cover change according to an embodiment of the present invention;

[0091] Fig. 9 It is a schematic diagram of kernel density analysis results of time lag and duration difference in a quantitative analysis method of the effect of human activities on vegetation cover change according to an embodiment of the present invention. DETAILED DESCRIPTION

[0092] To further illustrate each embodiment, the present invention provides drawings, which are part of the disclosure of the present invention and are mainly used to illustrate the embodiments and can be used in conjunction with the relevant descriptions in the specification to explain the operating principles of the embodiments. With reference to these contents, ordinary technicians in the field should be able to understand other possible implementations and advantages of the present invention. The components in the figures are not drawn to scale, and similar component symbols are generally used to represent similar components.

[0093] According to an embodiment of the present invention, a method for quantitatively analyzing the impact of human activities on vegetation cover changes is provided.

[0094] In the present invention, Southeast Asia (SEA) is selected as the research object. Southeast Asia includes 11 countries, namely Brunei, Indonesia, Cambodia, Laos, Myanmar, Malaysia, Philippines, Singapore, Thailand, East Timor and Vietnam, etc., and its basic situation is shown in Table 1. The total area of ​​Southeast Asia is 4.4825 million square kilometers, and the population in 2017 was 649 million, which is one of the most densely populated areas in the world. Southeast Asia is also one of the regions with the fastest population growth and urbanization, with an average annual population growth rate of 2.52% and an average annual GDP growth rate of 4.59% from 2008 to 2017. At the same time, Southeast Asia belongs to the equatorial rainy climate and tropical monsoon climate, with rich vegetation and small seasonal changes, and is an ideal area for studying the changes in vegetation cover caused by human activities.

[0095] Table 1: Basic information of Southeast Asian countries in the study area

[0096]

[0097] (1) VIIRS / DNB night light remote sensing data

[0098] VIIRS / DNB is a sensor carried on the Suomi National Polar-Orbiting Partnership (Suomi-NPP) satellite. It is mainly used to detect the intensity of the earth's night lights and provides monthly composite products of VIIRS / DNB. The spatial resolution of VIIRS / DNB data is about 500m, and the product coordinate system is the WGS84 geographic coordinate system, covering the global area of ​​75°N-65°S. The unit of the data is nanoWatts / cm2 / sr, and the data accuracy reaches eight decimal places. The monthly product is synthesized using the mean of the daily effective grid values, which can effectively remove noise such as stray light, lightning, and moonlight on the earth's surface. Although the data uses VIIRS Cloud Mask (VCM) to remove some clouds, there are still invalid data with zero or negative values. A quality control file with the suffix cf_cvg is also provided to indicate the number of days without cloud cover in the monthly composite product. Finally, 73 scene monthly composite products (April 2012-April 2018) were selected to construct the VIIRS time series, including tile No. 3 (covering the area 75N / 60E to 00N / 120E) and tile No. 6 (covering the area 00N / 60E to 65S / 120E).

[0099] (2) MODIS / NDVI vegetation remote sensing data

[0100] The spatial resolution of the 16-day synthetic NDVI (normalized difference vegetation index) vegetation cover product MOD13A1 is about 463m, and the data comes from the LAADS data distribution center of NASA in the United States. This product is synthesized from MODIS daily images, so the data has the characteristics of low cloud and rain coverage and high data accuracy. The value range of NDVI is -1 to 1, and the data accuracy reaches 4 decimal places. The higher the NDVI value (>0.2), the more and denser the vegetation coverage, among which 0.2 to 0.4 represents sparse vegetation distribution and 0.4 to 0.8 represents green vegetation coverage area. To ensure the accuracy of time series analysis, only NDVI values ​​of 0.2 to 0.8 are retained. MODIS / NDVI data also provides a reliability quality control layer, with values ​​of -1, 2, and 3 indicating no data, ice and snow coverage, cloud and rain coverage, etc. Covering the entire Southeast Asia requires 19 MODIS images, namely h26v06-v07, h27v06-v09, h28v06-v09, h29v06-v09, h30v07-v09, h31v09, and h32v09. All MODIS products are first spliced ​​based on the Mosaic module of the official MRT toolkit, and then the original sinusoidal projection is converted to WGS84 geographic coordinates through the Resample module of ArcGIS, and finally the spatial resolution of MODIS / NDVI data is resampled to be consistent with VIIRS / DNB data. The present invention has obtained a total of 140 MOD13A1 data (April 7, 2012-April 24, 2018).

[0101] In summary, this invention takes Southeast Asia as the study area and 500m×500m grid as the evaluation unit to carry out relevant analysis.

[0102] The present invention is further described with reference to the accompanying drawings and specific embodiments. Figure 1-7 As shown, according to one embodiment of the present invention, a method for quantitatively analyzing the effect of human activities on vegetation cover changes is provided, the method comprising the following steps:

[0103] S1. Use MODIS / NDVI remote sensing data to represent vegetation cover and VIIRS / DNB remote sensing data to represent human activity intensity, and construct two time series, MTS and VTS, respectively, through time series preprocessing technology;

[0104] The method of constructing two time series, MTS and VTS, respectively by using the time series preprocessing technology includes the following steps:

[0105] S11, construct VTS based on cf_cvg quality control file and perform smoothing;

[0106] The construction of the VIIRS / DNB time series (VTS) is based on the cf_cvg quality control file, which specifies the number of clear observation days. If the number of clear observation days corresponding to a grid is 0, the mean of the valid observation values ​​of the 6 neighborhoods before and after the point is used as the value of the grid; if the grid is at the beginning or end of the time series, the number of neighborhoods can be less than 6; if the nearby 6 neighborhoods do not contain valid values, the observation value of the grid is set to 0. The above operations can ensure the length and smoothness of the VTS.

[0107] S12, smoothing the MTS curve based on the Reliability layer;

[0108] MODIS / NDVI time series (MTS) can also be smoothed based on the Reliability layer. If the VTS and MTS of a grid are all 0, it means that the grid has no night light data or vegetation cover data, so the grid can be deleted directly without subsequent analysis.

[0109] S13, further smoothing the MTS and VTS through the Savitzky-Golay filter can reduce outliers and retain the original data to the maximum extent. Although the low-value outliers are removed, there are still some high-value outliers (spikes), especially in VTS, there are many high-value outliers;

[0110] S14, determine the threshold of the 6-neighborhood standard deviation by randomly selecting some samples and remove high-value outliers in MTS and VTS;

[0111] The calculation expression for determining the threshold of the 6-neighborhood standard deviation by randomly selecting some samples includes:

[0112] n∈(min(0,n-3),max(n+3,73 or 140))

[0113] |x m -μ|<=5σ,if|x i -μ|>5σ

[0114] Where n represents the set of 6 neighborhoods in VTS or MTS;

[0115] x n Represents the value of each element of the 6-neighborhood set;

[0116] μ represents the average value of the 6-neighborhood set;

[0117] σ represents the standard deviation of the 6-neighborhood set;

[0118] m represents the set of valid elements in the 6-neighborhood;

[0119] card represents the effective length of set n or set m;

[0120] x m represents the numerical value of each element in the set of valid values m, that is, |x m -μ| <= 5σ.

[0121] x i represents the outlier in the time series, that is, |x i -μ| > 5σ;

[0122] represents the corrected value of the outlier x i ;

[0123] S15. Obtain the standard deviation threshold through the calculation of the standard deviation of random sample points and manual visual interpretation, that is, if the absolute value of the difference between the raster value and the average value of the 6-neighborhood is greater than 5σ, then the raster value is re-assigned to the average value of all valid values within the 6-neighborhood.

[0124] S2. Construct a time series segmentation technique to achieve iterative segmentation of the time series, including the following steps:

[0125] S21. Use the least squares first-order linear regression algorithm to fit the MTS and VTS curves (where for VTS, the date is converted to 1 - 73, and for MTS, the date is converted to 1 - 140), and calculate the vertical distance between the time series data points and the fitted curve ( Figure 3 A);

[0126] S22. Set the inflection point according to the vertical distance between the fitted curve and the time series data points, including the following steps:

[0127] S221. Set the time series data point with the largest vertical distance from the fitted curve as the inflection point;

[0128] S222. When the length of the segmented time series corresponding to the inflection point is less than 3 months, that is, MTS is greater than 6 data points and VTS is greater than 3 data points, use the time series data point with the second largest vertical distance from the fitted curve as the inflection point ( Figure 3 B).

[0129] S23. Segment the current time series into two sub-time series through the inflection point, and the inflection point becomes the end point of the previous sub-time series and the start point of the subsequent sub-time series;

[0130] S24. The number of segmentation iterations determines the number of inflection points detected. A larger number of inflection points will reduce the residual, but it will also lead to overfitting problems and greater computational complexity. Therefore, the principle of minimizing the BIC (Bayesian Information Criterion) index is used to determine the optimal number of segmentation times for the time series.

[0131] S25, repeat the above segmentation steps ( Figure 3 C) Iterative segmentation is performed on MTS and VTS. For VTS, iterative segmentation is performed at most 3 times, and it can be divided into 8 sub-time series at most; for MTS, iterative segmentation is performed at most 4 times, and it can be divided into 16 sub-time series at most.

[0132] In addition, if Figure 3 The F statistic is used to test the statistical significance of each sub-time series. The threshold of the statistical significance level is 0.05, and local segmentation will continue until the statistical significance level is reached.

[0133] The above segmentation process may have the problem of over-segmentation for some time series, so a merging operation is needed to simplify the segmentation results of the time series.

[0134] S3. Construct time series merging and feature extraction technology, iterate through the sorting angle method to merge time series and extract time series features, including the following steps:

[0135] S31, using the sorting angle method, select two sub-time series with the smallest angle between adjacent fitting lines and the same change direction, and merge them first (such as Figure 3 D. When merging sub-time series, only merge those with the same change direction, that is, the fitting lines of the sub-time series are both increasing or decreasing);

[0136] S32. For the newly merged time series, the least squares first-order linear regression algorithm is used to fit the curve and calculate its overall trend, until any two adjacent sub-time series do not meet the requirements for merging (it should be noted that once the new fitting line and trend are calculated, the minimum adjacent fitting line angle needs to be recalculated);

[0137] S33, when the merging process is completed, new vertices and segmented sub-time series are formed, and the new vertices are fitted to obtain the overall trend characteristics ( Figure 3 F);

[0138] S34, for a non-head-tail vertex, the end point of the previous sub-time series and the start point of the next sub-time series may exist simultaneously on a date, and the specific value of the non-head-tail vertex is determined by calculating the average of the end point and the start point;

[0139] S35. After the segmentation and merging are completed, the definition of the change is changed and the descriptive features of the fitting curve are obtained.

[0140] The descriptive features include the core change period, the maximum trend feature (Maximum trend), the change duration feature (Duration) and the change start time (Start date);

[0141] Specifically, the core change period represents a sub-time series that is consistent with the overall trend change direction and has the largest absolute change value;

[0142] The maximum trend feature represents the absolute value of the change of the fitting curve within the core change period;

[0143] The change duration feature represents the change time of the core change period;

[0144] The change start time indicates the date when the core change period starts.

[0145] S4. By analyzing the time series characteristics of MTS and VTS, the calculation and spatial pattern analysis of MTS and VTS are performed to achieve quantitative analysis of the correlation between MTS and VTS, including the following steps:

[0146] S41. Because the Southeast Asia study area is large, it is not realistic to conduct time series analysis of MTS and VTS in the entire region. Therefore, by analyzing the overall trend characteristics of VTS, marking the areas with large changes in human activities as V-change areas and using them as the target areas for analysis can reduce the computational complexity;

[0147] S42, using a Gaussian fitting method to determine the threshold of the VTS change trend;

[0148] S43, performing a time series analysis of the vegetation coverage MTS in the V-change area;

[0149] S44. Calculate the time series analysis results and features of MTS and VTS, define the difference of Start date features as the time lag feature (Lag_effect); define the difference of Duration features as the duration difference feature (Duration_diff). The length of the 16-day synthetic MTS time series is 143. When calculating the above indicators, it needs to be unified to the VTS time series length of 70. The calculation expressions include:

[0150] Lag_effect=Start_date MODIS -Start_date VIIRS

[0151] Duration_diff = Duration MODIS -Duration VIIRS

[0152] Where, Lag_effect represents the time lag feature;

[0153] Duration_diff represents the duration difference feature;

[0154] Start_date MODIS Indicates the starting time characteristics of the time series change of MTS;

[0155] Start_date VIIRS Indicates the starting time characteristics of the time series change of VTS;

[0156] Duration MODIS Represents the duration characteristics of the time series changes of MTS;

[0157] Duration VIIRS Represents the time series variation duration characteristics of VTS.

[0158] S45, randomly selecting certain sample points manually, and marking the time series related features of the sample points;

[0159] S46. The spatial and temporal pattern of the obtained time series characteristics and results are analyzed by the global Moran index. The larger the positive value is, the more concentrated the spatial distribution is, while the larger the negative value is, the more dispersed the spatial distribution is.

[0160] S47. Calculate the Moran index of a single variable based on time series characteristics such as overall trend, maximum trend, time lag and duration;

[0161] S48. The obtained MTS and VTS maximum change trend characteristics are graded, and the VII_MOD_level indicator is constructed to integrate the two types of characteristics, so as to intuitively display the relationship between MTS and VTS, where the ten digits represent the high, medium and low classification of the MTS maximum change trend, and the ones digits represent the high, medium and low classification of the VTS maximum change trend. This indicator can more intuitively display the relationship between MTS and VTS.

[0162] Among them, the calculation expression for constructing the VII_MOD_level indicator to carry out the comprehensive calculation of the MTS maximum change trend classification indicator and the VTS maximum change trend classification indicator includes:

[0163]

[0164] VII_MOD_level=MOD_index×10+VIIRS_index

[0165] In the formula, MOD_index represents the classification index of the maximum change trend of MTS;

[0166] VIIRS_index represents the grading index of the maximum change trend of VTS;

[0167] VII_MOD_level represents the comprehensive grading index;

[0168] high indicates high classification;

[0169] medium means medium grade;

[0170] Low means low grade.

[0171] The following are specific examples, as well as analysis and verification of the results.

[0172] First, the number and proportion of grids under different VTS overall trends (Overall Trend) were counted. The results showed that the overall distribution showed a Gaussian distribution, that is, the proportion of grids with medium change trends was the highest, and as the change trend increased or decreased, the proportion of grids decreased exponentially. According to the results of the Gaussian fitting curve, for areas where the overall change trend of VIIRS was in the range of [-0.2, 0.4], they were set as areas where human activity changes were not significant; for areas where the overall change trend of VIIRS was (>0.4 or <-0.2), they were set as V-change areas where human activity changes were significant, which was also the core area for subsequent analysis ( Figure 4 Finally, 1.85 million 500m x 500m grids (approximately 8.64% of the total land area of ​​the region) were selected as V-change areas in the Southeast Asian study area, where the statistical significance of the VIIRS maximum change trend was less than 0.05.

[0173] 1. Analysis of changes in VTS human activities

[0174] V-change areas are mainly concentrated in and around large cities, with a Moran index of about 0.535 (p<0.01), showing strong spatial autocorrelation and spatial agglomeration. The overall trend of VTS shows a radial pattern with cities as the center and gradually attenuating outward ( Figure 5). At the same time, there are also areas within some cities where the overall trend of VTS is declining, such as the capitals of Indonesia and Brunei. Combined with Google historical image analysis, these areas are mainly urban renewal or ecological restoration areas. At the same time, an increase in VTS can be observed on roads and along the roads, mainly due to the improvement of road lighting infrastructure. As can be seen from Table 2, the overall trend of VTS in Singapore, Brunei, Vietnam and Malaysia has increased significantly. These countries all have large annual average GDP growth rates and population growth rates (Table 1). At the same time, there are also areas with significantly lower VTS in countries such as Singapore, Indonesia, and Malaysia, mainly because the migration of rural population to towns has led to a significant weakening of human activities in rural areas.

[0175] The maximum trend characteristics of VTS show a similar spatial pattern to the overall trend, with a Moran index of 0.420. The average duration of VTS core changes is 48.04 months, of which 39.95% of the regions last for more than 60 months, indicating that human activities are continuously increasing. At the same time, 74.73% of the regional VTS changes began in the research base period (2012.04), and the start time of the remaining regions showed a relatively uniform spatial distribution.

[0176] Table 2: VTS and overall VTS change trends in Southeast Asia and various countries

[0177]

[0178]

[0179] (II) Analysis of MTS changes in the V-change region

[0180] Based on the extraction of V-change regions, the MTS time series change analysis was carried out, and the statistical significance of the maximum trend characteristics was less than 0.05. The statistical results showed that the overall trend of VTS change was 2.12, and the overall trend of MTS change was 2.27×10 -2. This shows that in Southeast Asia, although human activities are increasing year by year, the vegetation coverage is generally increasing. Among them, Myanmar, Vietnam and Thailand have the highest increase in vegetation coverage. For countries where human activity VTS decreases, their vegetation coverage MTS increases. For areas with increased human activities, the average increase in VTS reached 2.38, while the average decrease in vegetation coverage MTS was 2.15×10-2. Using these two means as thresholds and combining the VTS thresholds in the V-change area, the overall trend change of VTS can be divided into three levels: increase (>2.38), slight increase (0.40-2.38), and decrease (<-0.2); the overall trend change of MTS is divided into three levels: increase (>2.15×10 -2 ), slightly increased (0-2.15×10 -2 ), decrease (<0). Therefore, the trends of VTS and MTS can be synthesized according to step S48 and divided into 9 categories. Figure 6 The grid statistical results under different VTS and MTS change trends are shown, where green indicates an overall increase in vegetation coverage MTS, purple indicates a slight increase in vegetation coverage, and brown indicates an overall decrease in vegetation coverage. The darker the color, the less the increase in VTS due to human activities.

[0181] 3. Analysis of the hysteresis effect of human activities on vegetation cover

[0182] The experimental results show that the lag effect of human activities on vegetation cover in Southeast Asia is more obvious, that is, the change of MTS lags behind the change of VTS by about 10.26 months, while previous studies have shown that the lag effect of non-human activity factors is about 2 years. This lag effect shows the feedback of vegetation cover changes on human activity changes. At the same time, the duration of MTS is about 9 months shorter than that of VTS, which is consistent with the obtained lag time. Although the spatial distribution of the time lag effect does not show an obvious pattern, areas close to the city center still show a shorter time lag effect ( Figure 7 ,Table 3), indicating that vegetation changes in densely populated areas respond faster. Among them, 37.10% of the grids did not show obvious time lag effects, while only 44.79% of the grids had MTS changes that lagged behind VTS changes.

[0183] Table 3: Analysis of time lag of MTS and VTS in Southeast Asia and various countries

[0184]

[0185]

[0186] The results of the difference analysis of VTS and MTS durations have a basically consistent spatial pattern with the results of the time lag analysis ( Figure 8 ), the duration difference is also affected by the study period. The present invention carries out binary kernel density analysis of time lag effect (Lag_effect) and duration difference (Duration_diff), and obtains the vertical line of time lag 0 value and the oblique line of negative linear correlation ( Fig. 9 ). The oblique line indicates that the end time of Lag_effect and Duration_diff are the same, and the change of MTS and VTS is likely to continue. The vertical line means that although the start time of VTS and MTS is the same, the duration of the change is different, and the change of MTS lasts longer.

[0187] The verification and analysis of the results adopts a combination of random sampling and manual assistance. For the MTS and VTS of 400 randomly selected grids, the change direction and change intensity of the time series analysis method proposed in the present invention are basically consistent with those of the time series analysis obtained manually, indicating the reliability of the time series and segmentation method proposed in the present invention. At the same time, through comparison with other surface cover products, the basic direction of changes in construction land and ecological land in the sample grids of these products is also basically consistent with the results obtained by the present invention.

[0188] In summary, with the help of the above technical solution of the present invention, by integrating time series segmentation, time series merging, spatial analysis and statistics and other technologies, a quantitative analysis of the impact of human activities on vegetation cover is achieved. First, taking Southeast Asia as the study area, MODIS / NDVI remote sensing data are used to characterize vegetation cover, and VIIRS / DNB remote sensing data are used to characterize the intensity of human activities, and MTS and VTS long time series are constructed respectively. Secondly, time series segmentation, time series merging and other technologies are used to extract time series features such as overall trend, maximum trend, start time of change, and duration. Finally, the time lag and duration difference of VTS and MTS are calculated, and their spatial patterns are analyzed to achieve quantitative analysis of human activities on vegetation cover changes, and the time series features of manual interpretation are analyzed to verify the effectiveness of the method.

[0189] The technical features of the above-described embodiments may be arbitrarily combined. To make the description concise, not all possible combinations of the technical features in the above-described embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0190] The above-mentioned embodiments only express several implementation methods of the present invention, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the invention patent. It should be pointed out that, for ordinary technicians in this field, several variations and improvements can be made without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the patent of the present invention shall be subject to the attached claims.

Claims

1. A quantitative analysis method of the impact of human activities on vegetation cover change. It is characterized in that The method comprises the following steps: S1. Use MODIS / NDVI remote sensing data to represent vegetation cover and VIIRS / DNB remote sensing data to represent human activity intensity, and construct two time series, MTS and VTS, respectively, through time series preprocessing technology; S2. Build time series segmentation technology to achieve iterative segmentation of time series; S3. Build time series merging and feature extraction technology, iterate through sorting angle method to merge time series and extract time series features; S4. By analyzing the time series characteristics of MTS and VTS, the calculation and spatial pattern analysis of MTS and VTS are carried out to achieve quantitative analysis of the correlation between MTS and VTS; The method of constructing two time series, MTS and VTS, respectively by using the time series preprocessing technology comprises the following steps: S11, construct VTS based on cf_cvg quality control file and perform smoothing; S12, smoothing the MTS curve based on the Reliability layer; S13, further smoothing the MTS and VTS through Savitzky-Golay filter; S14, determine the threshold of the 6-neighborhood standard deviation by randomly selecting some samples and remove high-value outliers in MTS and VTS; S15, obtaining the standard deviation threshold by calculating the standard deviation of random sample points and manually visually interpreting; The calculation expression for determining the 6-neighborhood standard deviation threshold by randomly selecting some samples includes: Where n represents the set of 6 neighborhoods in VTS or MTS; x n Represents the value of each element of the 6-neighborhood set; μ represents the average value of the 6-neighborhood set; σ represents the standard deviation of the 6-neighborhood set; m represents the set of valid elements in the 6-neighborhood; card represents the effective length of set n or set m; x m Represents the value of each element of the valid value set m, that is, |x m -μ|<=5σ; x i represents an outlier in the time series, i.e. |x i -μ|>5σ; Indicates the corrected value of the outlier x i ; The time series merging and feature extraction technology is constructed to merge the time series iteratively through the sorting angle method and extract the time series features, including the following steps: S31, using the sorting angle method, select two sub-time series with the smallest angle between adjacent fitting lines and the same change direction, and merge them preferentially; S32, for the newly merged time series, use the least squares first-order linear regression algorithm to fit the curve and calculate its overall trend, until any two adjacent sub-time series do not meet the requirements for merging; S33, when the merging process is completed, new vertices and segmented sub-time series are formed, and the new vertices are fitted to obtain overall trend characteristics; S34, for a non-head-tail vertex, the end point of the previous sub-time series and the start point of the next sub-time series may exist simultaneously on a date, and the specific value of the non-head-tail vertex is determined by calculating the average of the end point and the start point; S35, after the segmentation and merging are completed, defining and obtaining the descriptive features of the fitting curve; The method of analyzing the time series characteristics of MTS and VTS, calculating MTS and VTS and analyzing the spatial pattern, and realizing the quantitative analysis of the correlation between MTS and VTS comprises the following steps: S41. By analyzing the overall change trend characteristics of VTS, the areas with large changes in human activities are marked as V-change areas, and these areas are used as the target areas for analysis; S42, using a Gaussian fitting method to determine the threshold of the VTS change trend; S43, performing a time series analysis of the vegetation coverage MTS in the V-change area; S44, calculate the time series analysis results and characteristics of MTS and VTS; S45, randomly selecting certain sample points manually, and marking the time series related features of the sample points; S46, analyze the spatiotemporal pattern of the obtained time series characteristics and results through the global Moran index; S47. Calculate the Moran index of a single variable based on time series characteristics such as overall trend, maximum trend, time lag and duration; S48. The obtained maximum change trend characteristics of MTS and VTS are classified, and the VII_MOD_level index is constructed to integrate the two types of characteristics, so as to intuitively display the relationship between MTS and VTS.

2. According to claim 1, a quantitative analysis method for the effect of human activities on vegetation cover changes, It is characterized in that The time series segmentation technology is constructed to achieve iterative segmentation of the time series, including the following steps: S21, fitting the MTS and VTS curves using the least squares first-order linear regression algorithm, and calculating the vertical distance between the time series data points and the fitting curve; S22, setting an inflection point according to the vertical distance between the fitting curve and the time series data point; S23, dividing the current time series into two sub-time series by the inflection point, wherein the inflection point becomes the end point of the first sub-time series and the start point of the second sub-time series; S24. Use the principle of minimizing the BIC index to determine the optimal number of time series segmentation times; S25. Repeat the above segmentation steps to iteratively segment the MTS and VTS.

3. According to claim 2, a quantitative analysis method for the effect of human activities on vegetation cover changes, It is characterized in that The step of setting the inflection point according to the vertical distance between the fitting curve and the time series data point comprises the following steps: S221, setting the time series data point with the largest vertical distance from the fitting curve as the inflection point; S222: When the length of the segmented time series corresponding to the inflection point is less than 3 months, that is, the MTS is greater than 6 data points and the VTS is greater than 3 data points, the time series data point with the second largest vertical distance to the fitting curve is taken as the inflection point.

4. According to claim 3, a quantitative analysis method for the effect of human activities on vegetation cover changes, It is characterized in that The descriptive features include core change period, maximum trend feature, change duration feature and change start time; The core change period represents a sub-time series that is consistent with the overall trend change direction and has the largest absolute change value; The maximum trend feature represents the absolute value of the change of the fitting curve within the core change period; The change duration feature represents the change time of the core change period; The change start time indicates the date when the core change period starts.

5. According to claim 4, a quantitative analysis method for the effect of human activities on vegetation cover changes, It is characterized in that The calculation expressions for calculating the time series analysis results and characteristics of MTS and VTS include: Lageffect=Start_date MODIS -Start_date VIIRS Duration_diff=Duration MODIS -Duration VIIRS Where, Lag_effect represents the time lag feature; Duration_diff represents the duration difference feature; Start_date MODIS Indicates the starting time characteristics of the time series change of MTS; Start_date VIIRS Indicates the starting time characteristics of the time series change of VTS; Duration MODIS Represents the duration characteristics of the time series changes of MTS; Duration VIIRS Represents the time series variation duration characteristics of VTS.

6. A quantitative analysis method of the effect of human activities on vegetation cover changes according to claim 5, It is characterized in that The obtained MTS and VTS maximum change trend characteristics are classified, and the VII_MOD_level indicator is constructed to perform a comprehensive calculation expression of the two types of characteristics, including: VII_MOD_level=MOD_index×10+VIIRS_index In the formula, MOD_index represents the classification index of the maximum change trend of MTS; VIIRS_index represents the grading index of the maximum change trend of VTS; VII_MOD_level represents the comprehensive grading index; high indicates high classification; medium means medium grade; Low means low grade.

Images