A land surface vulnerability background assessment method based on time series remote sensing data

Abstract

The present invention discloses a method for evaluating the background of land surface vulnerability based on time series remote sensing data, which solves the technical problem that the existing technology cannot accurately reflect the continuous dynamic changes of surface processes. It includes: collecting basic remote sensing data for quantifying the four subsystems of "water, soil, air, and ecology" in the region; calculating the quantitative indicators of the four subsystems year by year; calculating the annual anomaly values ​​of the four subsystems based on the baseline data; and counting the frequency, intensity, and cumulative amount of positive or negative anomaly values ​​pixel by pixel. The frequency, intensity, and cumulative amount of the subsystem anomaly values ​​are normalized; the normalized result is reduced in dimension using the PCA method, and the first principal component of the frequency, intensity, and cumulative amount after dimensionality reduction is selected respectively, and the land surface vulnerability index is constructed after normalization. The present invention can combine the time series remote sensing data of the four subsystems of "water, soil, air, and ecology" and use the multi-year dynamic information of the subsystems to construct a land surface vulnerability index to evaluate the background state of land surface vulnerability.

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 research on the science of land surface systems. The land surface vulnerability background is the inherent physical properties of the land surface, including vegetation, soil, water bodies, and topography, as well as 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 capturing the temporal changes and complexity of the land surface system and its impact on the ecosystem.

[0003] However, existing environmental impact assessment methods often use static baseline references, capturing the land surface state at a single point in time or a typical historical period (e.g., averages and percentiles). Under this assumption, baseline conditions at different times may not accurately reflect the ongoing dynamic changes in surface processes, leading to a misunderstanding of the natural environmental background conditions. Therefore, it is necessary to consider the dynamic and evolving characteristics of the land surface system when conducting background assessments of land surface vulnerability. In light of this, the present invention proposes a novel background assessment method for land surface vulnerability 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 objectives, 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: raster remote sensing data were collected to calculate the surface water area ratio (SWAR), soil erosion (SE), drought index (SPEI), and vegetation status (VS) in the selected year within the region. The surface water area ratio (SWAR), soil erosion (SE), drought index (SPEI), and vegetation status (VS) were used to characterize the four subsystems of "water-soil-air-organism";

[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 data (z), 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 for calculating SPEI include precipitation (PRE) data and potential evapotranspiration (PET) data; the data for calculating VS include GIMSS-3G+DNVI data and MODIS-NDVI data for the Chinese region.

[0009] Step 2: Calculate the surface water area ratio (SWAR): Obtain pixels representing water body types from the land cover data year by year. Then, calculate the proportion of pixels representing water bodies in a specified grid unit year by year, which is called the surface water area ratio. Use 30-meter resolution land cover data to extract surface water pixels and calculate the proportion of water 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] The calculation formula of the RUSLE model is:

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

[0013] Where A is the amount of hydraulic erosion (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 topographic factor (slope length and slope), C is the vegetation cover and management factor, and P is the soil and water conservation measures factor. The LS, C, and P factors are dimensionless, and it is assumed that the K and LS factors do not change with interannual variation.

[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] Where 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 (m), 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 series is obtained by maximal synthesis from the monthly GIMSS-3G+DNVI dataset, and linear correction is performed using the MODIS-NDVI data of the same monthly period. Subsequently, the corrected monthly GIMSS-3G+DNVI dataset is generated based on the monthly linear correction coefficient. Finally, the inter-annual corrected data are averaged to calculate the annual vegetation status VS.

[0021] Step 6. Calculation of annual anomaly data: Select a period of slow growth in human activity intensity and calculate the reference benchmark for each subsystem separately.

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

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

[0024] I: the average value of positive or negative anomaly pixel values ​​for a single pixel 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 within each subsystem anomaly. Use principal component analysis to reduce the dimensionality 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] Where a, b, c, and d represent SWAR, SE, SPEI, and VS, respectively. LSVI stands for the Land Surface Vulnerability Index (LSVI), which ranges from 0 to 3. Low LSVI values ​​indicate low vulnerability, while high LSVI values ​​indicate high vulnerability. Using the natural breaks classification method, the LSVI is categorized into five levels (not vulnerable, slightly vulnerable, moderately vulnerable, highly vulnerable, and severely vulnerable) to characterize the baseline vulnerability of the land surface.

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

[0033] The present invention provides a land surface vulnerability background assessment method based on time series remote sensing data. This method not only takes into account the four subsystems of "water-soil-air-organism", but also utilizes the long-term dynamic information of the subsystems to evaluate the land surface vulnerability background. It is mainly suitable for environmental impact assessment in arid and semi-arid areas. It can combine regional climate and natural environment characteristics, consider the dynamic characteristics of time series of different subsystems, and construct land surface vulnerability background assessment indicators in a targeted manner, 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 following briefly introduces the drawings required for use in the embodiments or the description of the prior art. 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 any creative work.

[0035] Figure 1 This 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 land surface vulnerability background 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 of the present invention, all other embodiments obtained by ordinary technicians in this field without making 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 deemed 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 of the Hohhot-Baotou-Erdos-Yulin urban agglomeration is fragile, which can comprehensively demonstrate the results and applicability of the land surface vulnerability index to the vulnerability background of the region. The flow chart of the land surface vulnerability background evaluation method based on time series remote sensing data proposed in this invention is as follows: Figure 1 As shown, the following steps are included:

[0040] Step 1: Collection of remote sensing data representing 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-Erdos-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) data, soil erodibility (K) data, terrain factor (LS) data, 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 ), key plot length data (S), downwind distance of maximum wind erosion (z), climate factor (WF), soil erodible fraction (EF), soil crust factor (SCF), soil roughness (K′), and comprehensive vegetation factor (COG) data. The SPEI is calculated using precipitation (PRE) and potential evapotranspiration (PET). VS is calculated using 8 km GIMSS-3G+DNVI data and 250 m MODIS-NDVI data for the China region.

[0041] Step 2: Calculate the surface water area ratio (SWAR): Pixels representing water body types are obtained annually from 1990 to 2022 from the 30-meter resolution land cover data CLCD. Subsequently, the average value aggregation method is used to calculate the proportion of pixels representing water bodies in each 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 water erosion and wind erosion, which are calculated annually using the RUSLE model and the RWEQ model respectively. The RUSLE model calculation formula is:

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

[0044] Where A is the amount of hydraulic erosion (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 topographic factor (slope length and slope), C is the vegetation cover and management factor, and P is the soil and water conservation measures factor. The LS, C, and P factors are dimensionless, and the K and LS factors are assumed to remain constant from 1990 to 2022.

[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] Where 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 (m), 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. EF, SCF, and K′ are assumed to remain constant from 1990 to 2022.

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

[0051] Step 5. Construction of annual vegetation status (VS) sequence: First, the maximum 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 bilinear interpolation is used to resample to 1 km. Secondly, the 250-meter MODIS-NDVI in the Chinese region 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 monthly period (February 2000 to December 2022) to obtain the monthly linear correction coefficient. Then, based on the monthly linear correction coefficient, the corrected 1-km monthly GIMSS-3G+DNVI dataset from 1990 to 2022 is generated. Finally, the interannual vegetation state VS was calculated using the annual average of the calibrated 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: Calculate annual anomaly data: Select the period of slow growth in human activity intensity from 1990 to 1999 and calculate the reference baseline for each subsystem. Then, calculate the annual anomaly value for each subsystem using the following formula.

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

[0054] Where A i represents the anomaly value of year i (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 for each subsystem: For SWAR, SE, SPEI, and VS, calculate the frequency F, intensity I, and cumulative amount C of each negative or positive anomaly pixel by pixel. SWAR, SPEI, and VS count negative anomaly pixels, while 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 experiences 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 for 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 within each subsystem anomaly. Use principal component analysis to reduce the dimensionality 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] Where a, b, c, and d represent SWAR, SE, SPEI, and VS, respectively. LSVI stands for the Land Surface Vulnerability Index (LSVI), which ranges from 0 to 3. Low LSVI values ​​indicate low vulnerability, while high LSVI values ​​indicate high vulnerability. Using the natural breaks classification method, the LSVI is categorized into five levels (not vulnerable, slightly vulnerable, moderately vulnerable, highly vulnerable, and severely vulnerable) to characterize the baseline vulnerability of the land surface.

[0065] This technology is primarily applicable to assessing land surface vulnerability in arid and semi-arid regions. It combines the time series dynamics of the four subsystems of water, soil, air, and ecosystems, providing a new approach for scientifically assessing land surface vulnerability. Figure 3 The spatial distribution characteristics of the land surface vulnerability background assessment of the example area Hohhot-Baotou-Erdos-Yulin urban agglomeration are given. 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 understanding 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 status of land surface vulnerability in arid and semi-arid areas.

[0066] The basic principles, main features, and advantages of the present invention are shown and described above. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions are merely illustrative of the principles of the present invention. Various changes and modifications may be made to the present invention without departing from the spirit and scope of the present invention. Such changes and modifications are intended to fall within the scope of the present invention. The scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A land surface vulnerability background assessment method based on time series remote sensing data, characterized by: 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 within 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 proportion of pixels representing water bodies 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 inter-annual data to calculate the annual vegetation status; Step 6. Calculation of annual anomaly data: Select a period of slow growth in human activity intensity, calculate the reference baseline for each subsystem, and then calculate the annual anomaly value of each subsystem based on this reference baseline; 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 the surface water area ratio, soil erosion, drought index, and vegetation status pixel by pixel; Step 8. Construct the land surface vulnerability index: perform maximum and minimum normalization on F, I, and C within each subsystem anomaly, use principal component analysis 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, surface water pixels were extracted using 30-meter resolution land cover data, 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 three-month precipitation surplus and deficit, 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 calibrated 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 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