Land surface vulnerability background evaluation method based on time sequence remote sensing data

Abstract

The invention discloses a land surface vulnerability background evaluation method based on time sequence remote sensing data, and solves the technical problem that the continuous dynamic change of a land surface process cannot be accurately reflected in the prior art. Comprising the following steps: collecting basic remote sensing data of four quantified water-soil-gas-generation subsystems in a region; quantitative indexes of the four subsystems are calculated year by year; calculating year-by-year distance mean values of the four subsystems based on the reference data; and counting the frequency, the intensity and the cumulant of the positive anomaly or the negative anomaly pixel by pixel. Performing normalization processing on the frequency, the intensity and the cumulant of the sub-system anomaly; and carrying out dimensionality reduction on a normalization result by using a PCA method, respectively selecting first principal components after frequency, strength and cumulant dimensionality reduction, and constructing a land surface vulnerability index after normalization. According to the method, time sequence remote sensing data of four subsystems of'water-soil-gas-generation 'can be combined, and the land surface vulnerability index is constructed by utilizing multi-year dynamic information of the subsystems to evaluate the land surface vulnerability background state.

Description

Technical Field

[0001] The present invention relates to the field of environmental impact assessment, and in particular to a land surface vulnerability background assessment method based on time series remote sensing data. Background Art

[0002] The land surface system is one of the most complex and important subsystems of the earth's surface, and it is also the area most affected by human activities. With global warming, frequent extreme weather events and increasing human activities, key concepts such as vulnerability, resilience and adaptive capacity have become important aspects of in-depth exploration of the science of the land surface system. The land surface vulnerability background is the inherent physical properties of the land surface, including vegetation, soil, water bodies and topography, and their interactions with environmental factors such as climatic conditions, reflecting the overall state of the system under typical environmental conditions. As an important component of terrestrial ecosystems, changes in the land surface directly affect the function and structure of the ecosystem. Therefore, incorporating dynamic information into the vulnerability background assessment is crucial to capture the temporal changes and complexity of the land surface system and its impact on the ecosystem.

[0003] However, most existing environmental impact assessment methods use static baseline references, which capture the land surface state at a single time point or a typical historical period (e.g., averages and percentiles). Under this assumption, baseline conditions at different times may not accurately reflect the continuous dynamic changes in surface processes, leading to a misunderstanding of the natural environmental background conditions. Therefore, it is necessary to consider the dynamic and evolutionary characteristics of the land surface system to conduct a land surface vulnerability background assessment. In view of this, the present invention proposes a novel land surface vulnerability background assessment method based on time series remote sensing data. Summary of the invention

[0004] The purpose of the present invention is to provide a land surface vulnerability background assessment method based on time series remote sensing data to solve the technical problem that the existing technology cannot accurately reflect the continuous dynamic changes of surface processes.

[0005] To achieve the above object, the present invention provides the following technical solutions:

[0006] The present invention provides a land surface vulnerability background assessment method based on time series remote sensing data, which specifically includes the following steps:

[0007] Step 1: Subsystem remote sensing data collection: collect raster remote sensing data for calculating the surface water area ratio (SWAR), soil erosion (SE), drought index (SPEI) and vegetation status (VS) in the selected year in the region, and use the surface water area ratio (SWAR), soil erosion (SE), drought index (SPEI) and vegetation status (VS) to describe the four subsystems of "water-soil-air-life";

[0008] The data used to calculate SWAR include land cover data (LC); the data used to calculate SE include rainfall erosivity (R), soil erodibility (K), topographic factor (LS), vegetation cover and management factor (C), soil and water conservation measures factor (P), soil erosion and wind erosion data (S L ), maximum amount of wind-blown sand transferred (Q max ) data, key plot length data (S), downwind maximum wind erosion distance (z) data, climate factor (WF) data, soil erodible fraction (EF) data, soil crust factor (SCF) data, soil roughness (K') data, and comprehensive vegetation factor (COG) data for calculating SPEI; the data for calculating SPEI include precipitation (PRE) data and potential evapotranspiration (PET) data; the data for calculating VS include GIMSS-3G+DNVI data and China regional MODIS-NDVI data.

[0009] Step 2: Calculation of surface water area ratio (SWAR): Obtain pixels representing water body types from land cover data year by year; then, calculate the proportion of pixels representing water bodies in a specified grid unit in that year, called surface water area ratio, and use 30-meter resolution land cover data to extract surface water pixels, and calculate the proportion of water body pixels within a 1 km × 1 km grid.

[0010] Step 3: Calculation of soil erosion (SE): The sum of hydraulic erosion and wind erosion is used to represent soil erosion. The RUSLE model and RWEQ model are used to calculate hydraulic erosion and wind erosion, respectively.

[0011] Among them, the calculation formula of the RUSLE model is:

[0012] A=R×K×LS×C×P

[0013] Where A is the hydraulic erosion amount (t·hm -2 ·a -1 ), R is the rainfall erosivity factor (MJ·mm·hm -2 ·h -1 ·a -1 ), K is the soil erodibility factor (t·hm 2 ·h·MJ -1 mm -1 ·hm -2 ), LS is the terrain factor (slope length and slope), C is the vegetation cover and management factor, and P is the soil and water conservation measure factor. Among them, LS, C and P factors are dimensionless, and it is assumed that K factor and LS factor do not change with interannual changes.

[0014] The calculation formula of the RWEQ model is:

[0015]

[0016] S=150.71(WF×EF×SCF×K′×COG) -0.3711

[0017] Q max =109.8(WF×EF×SCF×K′×COG)

[0018] In the formula, S L is the amount of wind erosion (t·hm -2 ·a -1 ), Q max is the maximum amount of wind-blown sand transferred (kg·m -1 ), S is the length of the key plot (meters), z is the distance of maximum wind erosion in the downwind direction, which is 50 meters, and WF is the climate factor (kg·m -1 ), EF is the soil erodible fraction, SCF is the soil crust factor, K′ is the soil roughness factor, and COG is the comprehensive vegetation factor.

[0019] Step 4: Calculation of drought index (SPEI): Based on PRE data and PET data, the SPEI is calculated using the Python extension library climate_indices, and the calculation scale of SPEI is selected as three months.

[0020] Step 5. Calculation of vegetation status (VS): The monthly maximum NDVI sequence is obtained by maximum value synthesis from the GIMSS-3G+DNVI monthly dataset, and linear correction is performed using MODIS-NDVI data of the same period each month. Subsequently, the corrected monthly GIMSS-3G+DNVI dataset is generated based on the monthly linear correction coefficient. Finally, the annual vegetation status VS is calculated by averaging the corrected data within the year.

[0021] Step 6. Calculation of annual anomaly data: Select a period when the intensity of human activities increases slowly, and calculate the reference benchmark for each subsystem separately.

[0022] Step 7. Calculate the frequency, intensity and cumulative amount of positive / negative anomalies of the subsystem: For SWAR, SE, SPEI and VS, count the frequency F, intensity I and cumulative amount C of each negative or positive anomaly pixel by pixel. Among them, SWAR, SPEI and VS count negative anomaly pixels, SE counts positive anomaly pixels, and the specific meanings of the anomaly values ​​F, I and C are as follows:

[0023] F: the number of times a single pixel has positive or negative anomalies in the selected year;

[0024] I: average value of positive or negative anomaly pixel values ​​for individual pixels in the selected year;

[0025] C: The sum of the pixel values ​​with positive or negative anomalies for a single pixel in the selected year.

[0026] Step 8. Construct the Land Surface Vulnerability Index (LSVI): Normalize the maximum and minimum values ​​of F, I, and C in each subsystem anomaly. Use the principal component analysis method to reduce the dimension of F, I, and C, and use the first principal component PC1 as the intensity representation of land surface vulnerability. The formula is as follows:

[0027] F=PC1(F a , F b , F c , F d )

[0028] I=PC1(I a , I b , I c , I d )

[0029] C=PC1(C a , C b , C c , C d )

[0030] LSVI=F+I+C

[0031] In the formula, a, b, c and d represent SWAR, SE, SPEI and VS respectively. LSVI represents the land surface vulnerability value, and the range is [0, 3]. Low LSVI values ​​indicate low vulnerability levels, and high LSVI values ​​indicate high vulnerability levels. The natural break classification method is used to divide LSVI into five levels (non-vulnerable, slightly vulnerable, moderately vulnerable, highly vulnerable and severely vulnerable), which can characterize the background state of land surface vulnerability.

[0032] Based on the above technical solution, the embodiments of the present invention can at least produce the following technical effects:

[0033] The present invention provides a land surface vulnerability background assessment method based on time series remote sensing data. The method not only takes into account the four subsystems of "water-soil-air-life", but also uses the long-term dynamic information of the subsystems to evaluate the land surface vulnerability background. The method is mainly suitable for environmental impact assessment in arid and semi-arid areas. It can combine the characteristics of regional climate and natural environment, consider the dynamic characteristics of time series of different subsystems, and construct land surface vulnerability background assessment indicators in a targeted manner, thus providing an extremely effective new method for land surface vulnerability background assessment. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] 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 or the description of the prior art 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 the structures shown in these drawings without paying creative work.

[0035] Figure 1 It is a technical flow chart of the land surface vulnerability background assessment method based on time series remote sensing data;

[0036] Figure 2 is the geographical location and land cover type of the instance area;

[0037] Figure 3 It is the spatial distribution map of the background vulnerability of the land surface in the example area. DETAILED DESCRIPTION

[0038] The technical solutions in the embodiments of the present invention will be described clearly and completely below. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments in the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. In addition, the technical solutions between the various embodiments can be combined with each other, but they must be based on the ability of ordinary technicians in this field to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0039] Figure 2 The location and land cover type of the example area are given. The natural environment background of the Hohhot-Baotou-Ordos-Yulin urban agglomeration is fragile, which can comprehensively show the evaluation results and applicability of the land surface vulnerability index to the vulnerability background of the region. The flowchart of the land surface vulnerability background evaluation method based on time series remote sensing data proposed in the present invention is as follows: Figure 1 As shown, the following steps are included:

[0040] Step 1: Collect remote sensing data that characterize the four subsystems of "water-soil-air-biology": collect raster remote sensing data for calculating the surface water area ratio (SWAR), soil erosion (SE), drought index (SPEI) and vegetation status (VS) of the Hohhot-Baotou-Ordos-Yulin urban agglomeration from 1990 to 2022. The data for calculating SWAR include land cover data (LC). The data for calculating SE include rainfall erosivity (R), soil erodibility (K), terrain factor (LS), vegetation cover and management factor data (C), soil and water conservation measures factor data (P), soil erosion and wind erosion data (S L ), maximum amount of wind-blown sand transferred (Qmax ) data, key plot length data (S), downwind maximum wind erosion distance (z) data, climate factor (WF) data, soil erodible fraction (EF) data, soil crust factor (SCF) data, soil roughness (K′) data, and comprehensive vegetation factor (COG) data. The data used to calculate SPEI include precipitation (PRE) data and potential evapotranspiration (PET) data. The data used to calculate VS include 8 km GIMSS-3G+DNVI data and 250 m MODIS-NDVI data for China.

[0041] Step 2: Calculation of surface water area ratio (SWAR): Pixels representing water body types from 1990 to 2022 are obtained year by year from the 30-meter resolution land cover data CLCD; then, the average value aggregation method is used year by year to calculate the proportion of pixels represented as water bodies in the 1 km × 1 km grid unit in that year, which is called the surface water area ratio.

[0042] Step 3: Calculation of soil erosion (SE): Soil erosion includes hydraulic erosion and wind erosion, which are calculated year by year using the RUSLE model and the RWEQ model respectively. The calculation formula of the RUSLE model is:

[0043] A=R×K×LS×C×P

[0044] Where A is the hydraulic erosion amount (t·hm -2 ·a -1 ), R is the rainfall erosivity factor (MJ·mm·hm -2 ·h -1 ·a -1 ), K is the soil erodibility factor (t·hm 2 ·h·MJ -1 mm -1 ·hm -2 ), LS is the terrain factor (slope length and slope), C is the vegetation cover and management factor, and P is the soil and water conservation measure factor. Among them, LS, C and P factors are dimensionless, and it is assumed that K factor and LS factor will not change from 1990 to 2022 and are fixed values.

[0045] The calculation formula of the RWEQ model is:

[0046]

[0047] S=150.71(WF×EF×SCF×K′×COG) -0.3711

[0048] Q max =109.8(WF×EF×SCF×K′×COG)

[0049] In the formula, S L is the amount of wind erosion (t·hm -2 ·a -1 ), Q max is the maximum amount of wind-blown sand transferred (kg·m -1 ), S is the length of the key plot (meters), z is the distance of maximum wind erosion in the downwind direction, which is 50 meters, and WF is the climate factor (kg·m -1 ), EF is the soil erodible fraction, SCF is the soil crust factor, K′ is the soil roughness factor, and COG is the comprehensive vegetation factor. It is assumed that the EF factor, SCF factor, and K′ factor will not change from 1990 to 2022 and are fixed values.

[0050] Step 4: Calculation of drought index (SPEI): Based on the PRE data and PET data from 1990 to 2022, the Python extension library climate_indices is used to calculate the SPEI representing the precipitation surplus and deficit on a three-month scale.

[0051] Step 5, construction of annual vegetation status (VS) sequence: First, the maximum value synthesis method is used to obtain the monthly maximum NDVI sequence from the 8-km resolution GIMSS-3G+DNVI monthly dataset from 1990 to 2022, and the sequence is resampled to 1 km using bilinear interpolation. Secondly, the 250-meter MODIS-NDVI in China is resampled to 1 km using bilinear interpolation to ensure pixel alignment of the two datasets. Then, linear modeling is performed using the 1-km resolution GIMSS-3G+DNVI and MODIS-NDVI data for the same period of time (February 2000 to December 2022) to obtain the monthly linear correction coefficients. Next, based on the monthly linear correction coefficients, the corrected 1-km monthly GIMSS-3G+DNVI dataset from 1990 to 2022 is generated. Finally, the interannual vegetation status VS was calculated using the annually averaged 1 km GIMSS-3G+DNVI dataset from January 1990 to January 2000 and the 1 km MODIS-NDVI dataset from February 2000 to December 2022.

[0052] Step 6: Calculation of annual anomaly data: Select the period of slow growth of human activity intensity from 1990 to 1999, and calculate the reference benchmark for each subsystem. Then calculate the annual anomaly value of each subsystem according to the following formula.

[0053] A i =V i -Ave 1990-1999

[0054] In the formula, A i represents the anomaly value of the i-th year (2000-2022), Vi Represents the raster data for year i. 1990-1999 It is the average value of the raster data from 1990 to 1999.

[0055] Step 7: Calculate the frequency, intensity and cumulative amount of positive / negative anomalies of the subsystem: For SWAR, SE, SPEI and VS, count the frequency F, intensity I and cumulative amount C of each negative or positive anomaly pixel by pixel. Among them, SWAR, SPEI and VS count negative anomaly pixels, and SE counts positive anomaly pixels. The specific meanings of the anomaly values ​​F, I and C are as follows:

[0056] F: the number of times a single pixel has positive or negative anomalies during 2000-2022;

[0057] I: average value of positive or negative anomaly pixel values ​​for a single pixel during 2000-2022;

[0058] C: the sum of the pixel values ​​with positive or negative anomalies in a single pixel during 2000-2022;

[0059] Step 8: Construct the Land Surface Vulnerability Index (LSVI): Normalize the maximum and minimum values ​​of F, I, and C in each subsystem anomaly. Use the principal component analysis method to reduce the dimension of F, I, and C, and use the first principal component PC1 as the intensity representation of land surface vulnerability. The formula is as follows:

[0060] F=PC1(F a , F b , F c , F d )

[0061] I=PC1(I a , I b , I c , I d )

[0062] C=PC1(C a , C b , C c , C d )

[0063] LSVI=F+I+C

[0064] In the formula, a, b, c and d represent SWAR, SE, SPEI and VS respectively. LSVI represents the land surface vulnerability value, and the range is [0, 3]. Low LSVI values ​​indicate low vulnerability levels, and high LSVI values ​​indicate high vulnerability levels. The natural break classification method is used to divide LSVI into five levels (non-vulnerable, slightly vulnerable, moderately vulnerable, highly vulnerable and severely vulnerable), which can characterize the background state of land surface vulnerability.

[0065] The present invention is mainly applicable to the evaluation of land surface vulnerability background in arid and semi-arid areas. The technology combines the time series dynamics of the four subsystems of "water-soil-air-biology" and provides a new method for scientifically evaluating the land surface vulnerability background. Figure 3 The spatial distribution characteristics of the land surface vulnerability background assessment of the Hohhot-Baotou-Ordos-Yulin urban agglomeration are given as an example. Figure 3 It can be seen that the areas with high vulnerability and severe vulnerability levels basically overlap with the ecological governance areas, and the spatial distribution characteristics of land surface vulnerability are consistent with the background cognition of the natural environmental conditions in the region, indicating that the land surface vulnerability index proposed in this invention can more comprehensively evaluate the background state of land surface vulnerability in arid and semi-arid areas.

[0066] The above shows and describes the basic principles and main features of the present invention and the advantages of the present invention. It should be understood by those skilled in the art that the present invention is not limited to the above embodiments, and the above embodiments and descriptions are only for explaining the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, and these changes and improvements fall within the scope of the present invention to be protected. The scope of the present invention to be protected is defined by the attached claims and their equivalents.

Claims

1. A land surface vulnerability background assessment method based on time series remote sensing data, characterized in that: The following steps are involved: Step 1: Data collection to characterize the four subsystems of water, soil, air and ecology: collect raster remote sensing data for calculating the surface water area ratio SWAR, soil erosion SE, drought index SPEI and vegetation status VS in the selected year in the region; Step 2: Calculation of surface water area ratio (SWAR): Obtain pixels representing water body types from land cover data year by year; then calculate the ratio of pixels representing water body in a given grid unit year by year; Step 3: Calculation of soil erosion SE: Soil erosion includes hydraulic erosion and wind erosion, which are calculated using the RUSLE model and the RWEQ model respectively; Step 4: Calculation of drought index SPEI: Based on precipitation data and potential evapotranspiration data, the drought index is calculated using the Python extension library climate_indices; Step 5: Calculation of vegetation status VS: Use monthly high-resolution NDVI data and low-resolution NDVI data to construct long-term NDVI data, and average the annual data to calculate the annual vegetation status; Step 6: Calculation of annual anomaly data: Select the period of slow growth of human activity intensity, calculate the reference benchmark of each subsystem respectively, and calculate the annual anomaly value of each subsystem accordingly; Step 7: Calculate the frequency, intensity and cumulative amount of positive / negative anomalies of the subsystem: Count the frequency F, intensity I and cumulative amount C of negative or positive anomalies of surface water area ratio, soil erosion, drought index and vegetation status pixel by pixel; Step 8. Construct the land surface vulnerability index: normalize the maximum and minimum values ​​of F, I and C in each subsystem anomaly, use the principal component analysis method to reduce the dimension of F, I and C, and use the first principal component PC1 as the intensity representation of land surface vulnerability. The formula is as follows: F=PC1(F a ,F b ,F c ,F d ) I=PC1(I a ,I b ,I c ,I d ) C=PC1(C a ,C b ,C c ,C d ) LSVI=F+I+C Where a, b, c and d represent the proportion of surface water area, soil erosion, drought index and vegetation status, respectively, and LSVI represents the land surface vulnerability value.

2. The land surface vulnerability background assessment method based on time series remote sensing data according to claim 1 is characterized in that: In step 2, 30-m-resolution land cover data were used to extract surface water pixels, and the proportion of water pixels within a 1 km × 1 km grid was calculated.

3. The land surface vulnerability background assessment method based on time series remote sensing data according to claim 1 is characterized in that: In step 4, the SPEI, which represents the precipitation surplus and deficit on a three-month scale, is used to represent the drought status.

4. The land surface vulnerability background assessment method based on time series remote sensing data according to claim 1 is characterized in that: In step five, the high-resolution NDVI data and the low-resolution NDVI data need to be corrected using a linear method to ensure data consistency.

5. The land surface vulnerability background assessment method based on time series remote sensing data according to claim 1 is characterized in that: In step six, the benchmark of each subsystem is quantified by the average value of data from 1990 to 1999.

6. The land surface vulnerability background assessment method based on time series remote sensing data according to claim 1 is characterized in that: In step seven, the surface water area ratio, drought index and vegetation status are counted as negative anomaly pixels, and the soil erosion amount is counted as positive anomaly pixels.

7. The land surface vulnerability background assessment method based on time series remote sensing data according to claim 1 is characterized in that: In step eight, the value range of the LSVI is [0,3], and the LSVI is divided into five levels using the natural break classification method, namely, non-vulnerable, slightly vulnerable, moderately vulnerable, highly vulnerable and severely vulnerable.

Images