Ecological transition zone instant ecological vulnerability evaluation method and application thereof

Abstract

The present application relates to resource, workflow, personnel or project management technical field, disclose a kind of ecological transition zone instant ecological vulnerability evaluation method and its application.The present application, by the ecological transition zone ecological vulnerability evaluation index factor data of real-time acquisition from cloud platform, the current ecological vulnerability of ecological transition zone is calculated, so as to ensure that ecological vulnerability data can keep up with local environmental changes, with sufficient accuracy.Meanwhile, using the ecological vulnerability evaluation index factor data of years to calculate mutation year, then compare mutation year with reality event, and judge whether the selected evaluation index is reliable by this, so as to ensure that the calculation method itself is reliable.

Description

Technical Field

[0001] This invention relates to the fields of resource, workflow, personnel, or project management technology, and in particular to a method for real-time ecological vulnerability assessment of ecological transition zones and its application. Background Technology

[0002] Ecological vulnerability is a measure of the stability of a regional ecosystem. It refers to the variability of the ecological environment under various natural and human pressures within a certain spatial and temporal scale. It is the result of the combined effects of natural and anthropogenic factors, and this interaction often leads to developments in directions unsuitable for human use. In recent years, research on ecological vulnerability has gradually become a hot topic and focus of current research on ecological environment restoration and sustainable development, attracting considerable attention from scholars both domestically and internationally. Research on regional ecological vulnerability assessment not only allows us to understand the distribution of vulnerability levels at different spatial scales and the current state of the regional ecological environment, but also facilitates the rational allocation of ecological resources and the scientific regulation of human activities to achieve harmonious development between humans and nature, which is of great significance to the sustainable development of regional ecosystems.

[0003] For most regions, ecological vulnerability is relatively uniform and stable. As described in geography textbooks, a specific ecosystem corresponds to a specific level of ecological stability; for example, the tundra ecosystem, repeatedly emphasized in textbooks, is extremely fragile. Therefore, existing methods for calculating ecological vulnerability can accurately assess it and provide guidance. However, for ecological transition zones (commonly found in topographically fragmented areas and / or semi-arid regions), existing methods perform poorly, often yielding results that differ significantly from reality.

[0004] The inventors discovered that the reason existing methods for calculating ecological vulnerability perform poorly in ecological transition zones is not due to inaccuracy, but rather because the local environment has changed, rendering previously calculated ecological vulnerabilities inapplicable. Ecological transition zones contain multiple ecosystems within a small area, exhibiting highly uneven ecological vulnerability and significant volatility. A heavy rain can drastically alter the distribution of ecosystems, consequently changing ecological vulnerability. This makes existing ecological vulnerability calculation methods quite inaccurate for ecological transition zones, as the local ecological environment changes rapidly.

[0005] Taking the Wuliangsuhai Lake basin, which is the subject of this invention, as an example, it is located in a semi-arid region. Many historical events since modern times have had a significant impact on the ecological vulnerability of this area, and even the size of Wuliangsuhai Lake itself is constantly changing. In recent years, with rising temperatures and the northward shift of the precipitation belt, as well as environmental protection efforts, the ecological environment of the lake and its basin has tended to improve.

[0006] The aforementioned problems could be solved if the ecological vulnerability of ecological transition zones could be accurately predicted or the existing level of ecological vulnerability could be estimated in real time. However, the erroneous estimation phenomena found in relevant papers published in the top journal *Science* in recent years (DOI: 10.1126 / science.291.5503.438) indicate that a sufficiently effective method for predicting ecological vulnerability has not yet been mastered. The inventors attempted to program existing ecological vulnerability calculation methods to achieve real-time automatic acquisition of current ecological vulnerability results, but found that this also presents difficulties. It requires the use of some field data that cannot be obtained immediately (such as soil texture, vegetation management, etc.). Without field data, the accuracy of the results needs to be reassessed, but the accuracy of the assessment results also requires some field data that cannot be obtained immediately. Summary of the Invention

[0007] This invention provides a method for real-time ecological vulnerability assessment of ecological transition zones and its application.

[0008] The technical problem to be solved is that existing methods for assessing the ecological vulnerability of ecological transition zones are not adapted to the rapidly changing local ecological environment and have low accuracy.

[0009] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: a method for real-time ecological vulnerability assessment of ecological transition zones, wherein the ecological transition zones are located in topographically fragmented areas and / or semi-arid areas, and the assessment method includes the following steps: Step 1: Select a set of ecological vulnerability assessment index factors that can be obtained online in real time and have historical data, and obtain the real-time data and historical data of the ecological vulnerability assessment index factors; Step 2: Calculate the historical ecological vulnerability grading index of the ecological transition zone; Step 3: Conduct long-term time-series analysis of the changing trends and spatial distribution characteristics of ecological vulnerability, and compare whether the calculated abrupt change years match the actual events. If they do not match, repeat steps one to three until the calculated abrupt change years match the actual events. This is considered to have obtained a set of effective ecological vulnerability assessment index factors. Step 4: Using a set of effective ecological vulnerability assessment index factors, programmatically implement the automatic acquisition of current ecological vulnerability assessment index factor data and the real-time calculation of the ecological vulnerability classification index.

[0010] Furthermore, the effective set of ecological vulnerability assessment index factors obtained in step three includes elevation, topographic relief, landscape disturbance, annual average rainfall, atmospheric temperature, surface temperature, surface humidity, surface dryness, population density, vegetation cover, and economic density.

[0011] Furthermore, in step two, the calculation process for the ecological vulnerability grading index is as follows: Each ecological vulnerability assessment index factor was reprojected into a WGS84 ellipsoid and a UTM plane coordinate system, resampled to 1 km, and then subjected to Min-Max normalization. Spatial principal component analysis was used to normalize the first principal component (PC1) into positive indicators for each ecological vulnerability assessment index factor. The ecological vulnerability index is denoted as EVI, and is restored to the range of 0 to 1, EVI = 1 - PC1; Based on the vulnerability classification according to EVI, the ecological vulnerability classification index is obtained, denoted as EVGI.

[0012] Further, step two involves calculating an ecological vulnerability grading index for at least 10 years.

[0013] Furthermore, in step three, the Mann-Kendall mutation test is first used to perform a Mann-Kendall mutation test on the ecological vulnerability grading index to explore the mutation time points and extract the ecological vulnerability grading for important years. Then, the transformation trajectory method is used to represent the changes in the ecological vulnerability grading for important years, and the spatial distribution of ecological vulnerability in important years is analyzed intuitively. Finally, Sen+Menn-Kendall trend analysis is used to explore the changing trend of ecological vulnerability and its spatial distribution characteristics.

[0014] Furthermore, the ecological vulnerability grading indexes for each region calculated in step four are plotted on a color cloud map, with the same ecological vulnerability grading index marked with the same color on the color cloud map.

[0015] An application of an instant ecological vulnerability assessment method for ecological transition zones is proposed. This method is used to monitor the ecological vulnerability status of ecological transition zones in real time, and to recommend the orderly optimization of industrial and agricultural production activities in areas where the ecological vulnerability grading index reaches or exceeds level three. It also recommends prioritizing ecological restoration in areas where the ecological vulnerability grading index reaches or exceeds level four.

[0016] Furthermore, the factor detector and interaction detector in the geographic detector were used to analyze the factors of each ecological vulnerability assessment index, explore the driving force of each ecological vulnerability assessment index factor on the overall ecological vulnerability, and prioritize intervention on factors with high single-factor driving force and factors with high interaction synergy when withdrawing industrial and agricultural production activities and carrying out ecological restoration.

[0017] The present invention provides a method for real-time ecological vulnerability assessment of ecological transition zones and its application, which, compared with existing technologies, has the following advantages: In this invention, the current ecological vulnerability of the ecological transition zone is calculated by acquiring real-time ecological vulnerability assessment index factor data from a cloud platform, thereby ensuring that the ecological vulnerability data keeps pace with local environmental changes and has sufficient accuracy. Simultaneously, the year of abrupt change is calculated using historical ecological vulnerability assessment index factor data, and then the year of abrupt change is compared with the actual event to determine whether the selected assessment indicators are reliable, thus ensuring the reliability of the calculation method itself. Attached Figure Description

[0018] Figure 1 This is a flowchart of a method for real-time ecological vulnerability assessment of ecological transition zones according to the present invention; Figure 2 The figure shows the Menn-Kendall mutation curve. The two dashed lines represent the statistic Z when α=0.05, with a value of ±1.96. The mutation years in the figure are 2009, 2013, and 2017. Figure 3 Spatial distribution map of ecological vulnerability index in different years; Figure 4 The diagram shows the trajectory changes for different levels of ecological vulnerability. The number of digits in the trajectory code from left to right represents the ecological vulnerability levels in 2000, 2009, 2013, 2017, and 2019, respectively. Figure 5 This is a trend analysis chart of ecological vulnerability. Figure 6 This diagram illustrates the impact of two-factor interactions on ecological vulnerability. The values ​​in the diagram represent the q-values ​​of the two-factor interactions, i.e., q(X1∩X2). Detailed Implementation

[0019] This invention has been applied to the Wuliangsuhai River Basin in Inner Mongolia as an example, and the details are as follows: like Figure 1 As shown, an instant ecological vulnerability assessment method for ecological transition zones, located in topographically fragmented and / or semi-arid regions, is presented. The assessment method includes the following steps: Step 1: Select a set of ecological vulnerability assessment index factors that can be obtained online in real time and have historical data, and obtain the real-time data and historical data of the ecological vulnerability assessment index factors; Step 2: Calculate the historical ecological vulnerability grading index of the ecological transition zone; Step 3: Conduct long-term time-series analysis of the changing trends and spatial distribution characteristics of ecological vulnerability, and compare whether the calculated abrupt change years match the actual events. If they do not match, repeat steps one to three until the calculated abrupt change years match the actual events. This is considered to have obtained a set of effective ecological vulnerability assessment index factors. Step 4: Using a set of effective ecological vulnerability assessment index factors, programmatically implement the automatic acquisition of current ecological vulnerability assessment index factor data and the real-time calculation of the ecological vulnerability classification index.

[0020] Step three yielded a set of effective ecological vulnerability assessment indicators, including elevation, topographic relief, landscape disturbance, average annual rainfall, atmospheric temperature, surface temperature, surface humidity, surface dryness, population density, vegetation cover, and economic density. These ecological vulnerability assessment indicators are empirically determined without a standardized selection process, and their reliability in calculating ecological vulnerability is not guaranteed; therefore, verification is essential. However, for ecological transition zones, existing technologies lack effective verification methods. Therefore, historical data was used to calculate abrupt change years for verification.

[0021] The sources of this data are as follows: Raw data are Landsat remote sensing imagery; MODIS products and elevation data are directly referenced from the GEE platform (https: / / developers.google.com / earth-engine / datasets). Land cover and atmospheric temperature data are from the Zendo community (https: / / zenodo.org / ). Population density data are from the NASA Center for Socioeconomic Data Applications (https: / / sedac.ciesin.columbia.edu / data / collection / gpw-v4), and economic density and average annual rainfall data are from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (https: / / www.resdc.cn). See the table below for details: In step two, the calculation process for the ecological vulnerability grading index is as follows: After obtaining the corresponding indices, the 11 evaluation index factors were reprojected onto the WGS84 ellipsoid and UTM plane coordinate system, resampled to 1 km, and subjected to Min-Max normalization. Then, spatial principal component analysis was used to normalize the first principal component (PC1) of the 11 evaluation indicators. The normalized value is related to the degree of ecological vulnerability; a higher normalized value indicates a better ecological environment. To make the results more objective and scientific, EVI = 1 - PC1 was calculated and restored to the range of 0 to 1 to obtain the Ecological Vulnerability Index (EVI) for the Wuliangsuhai Lake Basin. Finally, based on the table below, vulnerability grading was calculated to obtain the Ecological Vulnerability Grading Index (EVGI).

[0022] In the formula, i Vulnerability is categorized into five levels, where S represents the total area of ​​the study region. i Vulnerability level is i The area.

[0023] The reason for grading here is to ensure that the calculated results are meaningful; a single numerical value cannot provide guidance.

[0024] Step two involves calculating an ecological vulnerability grading index covering at least 10 years. If the years are shorter than this value, it becomes difficult to obtain sufficient marker events indicating ecological vulnerability abrupt changes for verification. In this embodiment, steps one through three incorporate data from the period 2000 to 2019 into the calculation.

[0025] In step three, the Mann-Kendall mutation test is first used to perform a Mann-Kendall mutation test on the ecological vulnerability grading index to explore the mutation time points and extract the ecological vulnerability grading for important years. Then, the transformation trajectory method is used to represent the changes in ecological vulnerability grading in important years, and the spatial distribution of ecological vulnerability in important years is analyzed intuitively. Finally, Sen+Menn-Kendall trend analysis is used to explore the changing trend of ecological vulnerability and its spatial distribution characteristics.

[0026] The formula for the Mann-Kendall mutation test is as follows: in style k When =1, S1=0. r i Indicates the first i time pointx i > x j (1≤ j ≤ i The cumulative number of ) x i Indicates time point i EVGI value. UF k It follows a standard normal distribution. UF 1 = 0. For a given significance level α = 0.05, if UF k When the value is greater than 0, the original sequence shows an upward trend. UF k When <0, the original sequence shows a downward trend. Then, an inverse time series is defined. x n , x n-1 ,…, x 1, for UB k Repeat the same operation to satisfy UB k =- UF k , UB n =0. UF k and UB k The time of the intersection within the confidence level is the time when the mutation occurs.

[0027] Menn-Kendall mutation curve Figure 2 As shown, the mutation years are 2009, 2013, and 2017; In reality: In 2009, the lake area repeatedly experienced large-scale outbreaks of yellow algae, a sign of ecological vulnerability and sudden change. At the same time, the area of ​​artificial impermeable layers expanded rapidly, and the area of ​​wetlands, forests and grasslands decreased, leading to a gradual increase in the fluctuation of ecological vulnerability in the watershed.

[0028] In 2013, the development of irrigation areas in the western Ulan Buh Desert led to a sharp decrease in the number of plant species in the grasslands, a decline in vegetation cover and stability, and frequent changes in landscape patches such as bare desert soil base and vegetation. This increased the fluctuation of landscape disturbance and affected the changes in the ecological vulnerability index.

[0029] In 2017, the preliminary ecological restoration project was nearing completion. In 2018, the Wuliangsuhai Basin Ecological Protection and Restoration Project was included in the third batch of national pilot projects for ecological protection and restoration of mountains, rivers, forests, fields, lakes and grasslands. The bare desert soil north of Wulashan has basically completed vegetation restoration, and the ecological vulnerability has been restored to stability.

[0030] This shows that the years of the mutations coincide with real-world marker events indicating ecological vulnerability. From Figure 3-5 It allows for a more intuitive view of the spatial distribution characteristics of the changing trends in ecological vulnerability in the Wuliangsuhai watershed area.

[0031] The formula for the trajectory transformation method is as follows: T ij For the first in the trajectory layer i Line 1 j The trajectory code of the column pixels has no mathematical meaning; n The number of time points; It is the code for the land use / cover type at each time point for a given pixel.

[0032] The Sen+Menn-Kendall trend analysis formula is as follows: In the formula, Q EVGI Indicates changes in the trend of time series data. EVGI j and EVGI i They are time points respectively j and time point i The ecological vulnerability index value, when Q EVGI The trend is upward when | >0.0005, when | EVGI The trend remains stable when |≤0.0005, when Q EVGI The trend decreases when <-0.0005. n When the value is ≥10, S approximately follows a standard normal distribution, and the trend test is performed using the Z statistic. The significance of Z is tested at a confidence level of α=0.05. A significant change occurs when |Z|>1.96, and a non-significant change occurs when |Z|≤1.96.

[0033] The ecological vulnerability grading indices for each region calculated in step four are plotted on a color cloud map, with the same ecological vulnerability grading index marked by the same color. This provides a clear visual understanding of the local ecological vulnerability status, facilitating guidance for practical production and daily life activities.

[0034] An application of an instant ecological vulnerability assessment method for ecological transition zones is proposed. This method is used to monitor the ecological vulnerability status of ecological transition zones in real time, and to recommend the orderly optimization of industrial and agricultural production activities in areas where the ecological vulnerability grading index reaches or exceeds level three. It also recommends prioritizing ecological restoration in areas where the ecological vulnerability grading index reaches or exceeds level four.

[0035] The factor detector and interaction detector in the geographic detector were used to analyze the factors of each ecological vulnerability assessment index, exploring the driving force of each factor on the overall ecological vulnerability. When withdrawing industrial and agricultural production activities and carrying out ecological restoration, priority was given to intervening in factors with high single-factor driving force and factors with high interaction synergy. The reason for prioritizing intervention in high-ranking factors is to ensure the effectiveness of the intervention.

[0036] Among them, factor vulnerability is mainly used to detect the magnitude of the influence of a single factor on the ecological vulnerability of a watershed. The calculation method of the factor detector is as follows: In the formula: q Impact Factor i Explanatory power of factors related to ecological vulnerability; n For sample size; L It is the number of index factor categories; N h and They are respectively h The variance of layer sample size and ecological vulnerability. q =[0,1], q A higher value indicates a higher impact factor. i The stronger the explanatory power of the factors on the ecological vulnerability of the watershed.

[0037] from Figure 6 It can be seen that the interactive detection results show that the factors have interactive synergistic effects. Among the five factors with strong single-factor explanatory power, surface humidity and landscape disturbance are nonlinearly enhanced, while the other two factors are bilinearly enhanced.

[0038] The programming mentioned in step four of this embodiment has been completed using JavaScript and can run in Google Earth Engine. The detailed code is as follows: / **** Start of imports. If edited, may not auto-convert in theplayground. **** / var RFi = yes.Image("users / young0801 / DEM_RFi"), wuliangsu2 = yes.FeatureCollection("users / yufanli0801 / wuliangsuROIs / wuliangsu"), wuliangsu = yes.FeatureCollection("users / yunfanli0801 / WuliangsuROIs / LandNoWater"); FVC2001 = yes.Image("users / young0801 / FVC2001_2010a / FVC2013"), GDP2001 = yes.Image("users / child0801 / GDP00_10 / GDP13"), a2001 = yes.Image("users / youngli0801 / Landscape2000_2010a / Landscape13"), Pop01 = yes.Image("users / young0801 / Population / Pop13_49N"), b2001 = yes.Image("users / young0801 / Rainfall2000_2015a / 2013"), Temp2001 = yes.Image("users / young0801 / Temperature00_10 / Mean2013"); / ***** End of imports. If edited, may not auto-convert in theplayground. ***** / var roi = geometry().bounds(); Map . centerObject ( row , 9 ) ; function removeCloud(image){ var qa = image.select('BQA') var cloudMask = qa.bitwiseAnd(1<<4).eq(0) var cloudShadowMask = qa.bitwiseAnd(1<<8).eq(0) var valid = cloudMask.and(cloudShadowMask) return image.updateMask(valid)} var L8 = ee.ImageCollection("LANDSAT / LC08 / C02 / T1_L2") .filterBounds(roi) .filterDate('2013-01-01', '2020-12-31') .filterMetadata('CLOUD_COVER', 'less_than',50) .map(function(image){ return image.set('year', ee.Image(image).date().get('year'))}); .map(removeCloud); var L8imgList = ee.List([]); for(var a = 2013; a<2020; a++){ var img = L8.filterMetadata('year', 'equals', a).median().clip(roi); var L8img = img.set('year', a); L8imgList = L8imgList.add(L8img); } var L8imgCol = ee.ImageCollection(L8imgList) .map(function(img){ return img.clip(roi); }); L8imgCol = L8imgCol.map(function(img){ / *DEM* / var dem = ee.Image('CGIAR / SRTM90_V4'); dem = dem.clip(wuliangsu); function CmpDEM(Band,region,scale){ var num = Band.reduceRegion({ reducer:ee.Reducer.minMax(), geometry:region, scale:scale, bestEffort:true, }); var min = ee.Number(num.get("elevation_min")); var max = ee.Number(num.get("elevation_max")); var Forward=Band.subtract(min).divide(max.subtract(min)); return Forward.unitScale(0,1).rename('Forward'); } var DEM=CmpDEM(dem,wuliangsu,9); print(DEM); Map.addLayer(DEM); img = img.addBands(DEM.rename('DEM')); / *RFi* / function CmpRFi(Band,region,scale){ var num = Band.reduceRegion({ reducer:ee.Reducer.minMax(), geometry:region, scale:scale, bestEffort:true, }); var min = ee.Number(num.get("b1_min")); var max = ee.Number(num.get("b1_max")); var Forward=Band.subtract(min).divide(max.subtract(min)); return Forward.unitScale(0,1).rename('Forward'); } var RFI=CmpRFi(RFi,wuliangsu,9); print(RFI); Map.addLayer(RFI); img = img.addBands(RFI.rename('RFI')); / *FVC* / var CmpNDVI =ee.ImageCollection("LANDSAT / LT05 / C02 / T1_L2") ; CmpNDVI = CmpNDVI.filterBounds(wuliangsu) .filterDate('2001-01-01', '2001-12-31') .map(function(image){ varimage_clip =image.clip(wuliangsu); returnimage_clip;}) .map(function(image){ var ndvi = image.normalizedDifference(["SR_B4","SR_B3"]); returnimage.addBands(ndvi.rename("NDVI"));}) .qualityMosaic('NDVI'); var bestndvi =CmpNDVI.select('NDVI').rename('bestNDVI'); Map.addLayer(bestndvi,{palette: ['red', 'green', 'blue']},'bestndvi'); functionCmpFVC(BestVI,region,scale){ var num = BestVI.reduceRegion({ reducer:ee.Reducer.percentile([5,95]), geometry:region, scale:scale, maxPixels:1e13 }); var min = ee.Number(num.get("bestNDVI_p5")); var max = ee.Number(num.get("bestNDVI_p95")); print(min); print(max); var greaterPart = BestVI.gt(max); var lessPart = BestVI.lt(min); var middlePart =ee.Image(1).subtract(greaterPart).subtract(lessPart); var tempf1=BestVI.subtract(min).divide(max.subtract(min)); var FVC=ee.Image(1).multiply(greaterPart).add(ee.Image(0).multiply(lessPart)) .add(tempf1.multiply(middlePart)); return FVC.rename('FVC'); } var FVC=CmpFVC(bestndvi,wuliangsu,9); Map.addLayer(FVC2001,{palette: ['red', 'green', 'blue']},'FVC'); img = img.addBands(FVC2001.rename('FVC')); / *Land2001* / a2001 = a2001.clip(wuliangsu2); function CmpLand(Band,region,scale){ var num = Band.reduceRegion({ reducer:ee.Reducer.minMax(), geometry:region, scale:scale, bestEffort:true, }); var min = ee.Number(num.get("b1_min")); var max = ee.Number(num.get("b1_max")); var Forward=Band.subtract(min).divide(max.subtract(min)); return Forward.unitScale(0,1).rename('Forward'); } var Land2001=CmpLand(a2001,wuliangsu2,9); print(Land2001); Map.addLayer(Land2001); img = img.addBands(Land2001.rename('Land')); / *Population* / Pop01 = Pop01.clip(wuliangsu); function CmpPop(Band,region,scale){ var num = Band.reduceRegion({ reducer:ee.Reducer.minMax(), geometry:region, scale:scale, bestEffort:true, }); var min = ee.Number(num.get("b1_min")); var max = ee.Number(num.get("b1_max")); var Forward=Band.subtract(min).divide(max.subtract(min)); return Forward.unitScale(0,1).rename('Forward'); } var Pop=CmpPop(Pop01,wuliangsu,9); print(Pop); Map.addLayer(Pop); img = img.addBands(Pop.rename('Pop')); / *GDP* / GDP2001 = GDP2001.clip(wuliangsu); function CmpBackward(Band,region,scale){ var num = Band.reduceRegion({ reducer:ee.Reducer.minMax(), geometry:region, scale:scale, bestEffort:true, }); var min = ee.Number(num.get("b1_min")); var max = ee.Number(num.get("b1_max")); var Forward=Band.subtract(min).divide(max.subtract(min)); var Backward = Forward.expression('1-Forward', {'Forward':Forward}); return Backward.unitScale(0,1).rename('Backward'); } var GDP=CmpBackward(GDP2001,wuliangsu,9); print(GDP); Map.addLayer(GDP); img = img.addBands(GDP.rename('GDP')); / *Temperature* / Tem2001 = Tem2001.clip(wuliangsu); function CmpTem(Band,region,scale){ var num = Band.reduceRegion({ reducer:ee.Reducer.minMax(), geometry:region, scale:scale, bestEffort:true, }); var min = ee.Number(num.get("b1_min")); var max = ee.Number(num.get("b1_max")); var Forward=Band.subtract(min).divide(max.subtract(min)); return Forward.unitScale(0,1).rename('Forward'); } var Temperature=CmpRainfall(Tem2001,wuliangsu,9); print(Temperature); Map.addLayer(Temperature); img = img.addBands(Temperature.rename('Tem')); / *Rainfall2001* / b2001 = b2001.clip(wuliangsu); function CmpRainfall(Band,region,scale){ var num = Band.reduceRegion({ reducer:ee.Reducer.minMax(), geometry:region, scale:scale, bestEffort:true, }); var min = ee.Number(num.get("b1_min")); var max = ee.Number(num.get("b1_max")); var Forward=Band.subtract(min).divide(max.subtract(min)); return Forward.unitScale(0,1).rename('Forward'); } var Rainfall2001=CmpRainfall(b2001,wuliangsu,9); print(Rainfall2001); Map.addLayer(Rainfall2001); img = img.addBands(Rainfall2001.rename('Rainfall')); / *Wet* / var Wet = img.expression('B*(0.1511) + G*(0.1973) + R*(0.3283)+ NIR*(0.3407) +SWIR1*(-0.7117) + SWIR2*(-0.4559)',{ 'B': img.select(['SR_B2']), 'G': img.select(['SR_B3']), 'R': img.select(['SR_B4']), 'NIR': img.select(['SR_B5']), 'SWIR1': img.select(['SR_B6']), 'SWIR2': img.select(['SR_B7']) }); img = img.addBands(Wet.rename('WET')); / *NDVI* / var ndvi = img.normalizedDifference(["SR_B4","SR_B3"]); img = img.addBands(ndvi.rename('NDVI')); / *LST_MODIS* / var lst_MODIS = ee.ImageCollection('MODIS / 006 / MOD11A1') .filterBounds(roi) .map(function(image){ var image_clip =image.clip(roi); return image_clip;}) .filter(ee.Filter.date('2013-01-01', '2013-12-31')) .select(['LST_Day_1km', 'LST_Night_1km']) .mean(); var lst_mean = lst_MODIS.expression('((Day + Night) / 2)',{ 'Day': lst_MODIS.select(['LST_Day_1km']), 'Night': lst_MODIS.select(['LST_Night_1km'])}) .rename("LST"); function CmpForward(Band,region,scale){ var num = Band.reduceRegion({ reducer:ee.Reducer.minMax(), geometry:region, scale:scale, bestEffort:true, }); var min = ee.Number(num.get("LST_min")); var max = ee.Number(num.get("LST_max")); var Forward=Band.subtract(min).divide(max.subtract(min)); return Forward.unitScale(0,1).rename('Forward'); } var LST=CmpForward(lst_mean,roi,9); print(LST); Map.addLayer(LST); img = img.addBands(LST.rename('LST')); ndbsi = ( ibi + si ) / 2 var ibi = img.expression('(2 * SWIR1 / (SWIR1 + NIR) - (NIR / (NIR+ RED) + GREEN / (GREEN + SWIR1))) / (2 * SWIR1 / (SWIR1 + NIR) + (NIR / (NIR+ RED) + GREEN / (GREEN + SWIR1)))', { 'SWIR1': img.select('SR_B6'), 'NIR': img.select('SR_B5'), 'NETWORK': img.select('SR_B4'), 'GREEN': img.select('SR_B3') }) var si = img.expression('((SWIR1 + RED) - (NIR + BLUE)) / ((SWIR1 +RED) + (NIR +BLUE))', { 'SWIR1': img.select('SR_B6'), 'NIR': img.select('SR_B5'), 'NETWORK': img.select('SR_B4'), 'BLUE': img.select('SR_B2') }) var ndbsi = (ibi.add(si)).divide(2) return img . addBands ( ndbsi . rename ( 'NDBSI ' ))} ) ; var bandNames = ["NDBSI","WET","LST","DEM","POP","GDP","RFI","FVC","Temperature","Rainfall","Land"]; L8imgCol = L8imgCol . select ( bandNames ) / / ingredient(substitute ingredient ingredient) var img_normalize = function(img){ var minMax = img.reduceRegion({ reducer:ee.Reducer.minMax(), geometry: roi, scale: 1000, maxPixels: 10e13, }) var year = img.get('year') var normalize = ee.ImageCollection.fromImages( img.bandNames().map(function(name){ name = ee.String(name); var band = img.select(name); return band.unitScale(ee.Number(minMax.get(name.cat('_min'))),ee.Number(minMax.get(name.cat('_max'))));}) ).toBands().rename(img.bandNames()).set('year', year); return normalize; } var imgNorcol = L8imgCol.map(img_normalize); / / Principal Component var pca = function(img){ var bandNames = img.bandNames(); var region = roi; var year = img.get('year') / / Mean center the data to enable a faster covariance reducer / / and an SD stretch of the principal components. var meanDict = img.reduceRegion({ reducer: ee.Reducer.mean(), geometry: region, scale: 1000, maxPixels: 10e13 }); var means = ee.Image.constant(meanDict.values(bandNames)); var centered = img.subtract(means).set('year', year); / / This helper function returns a list of new band names. var getNewBandNames = function(prefix, bandNames){ var seq = ee.List.sequence(1, 11); / / var seq = ee.List.sequence(1, bandNames.length()); return seq.map(function(n){ return ee.String(prefix).cat(ee.Number(n).int()); }); }; / / This function accepts mean centered imagery, a scale and / / a region in which to perform the analysis. It returns the / / Principal Components (PC) in the region as a new image. var getPrincipalComponents = function(centered, scale, region){ var year = centered.get('year') var arrays = centered.toArray(); / / Compute the covariance of the bands within the region. var covar = arrays.reduceRegion({ reducer: ee.Reducer.centeredCovariance(), geometry: region, scale: scale, bestEffort:true, maxPixels: 10e13 }); / / Get the 'array' covariance result and cast to anarray. / / This represents the band-to-band covariance within theregion. var covarArray = ee.Array(covar.get('array')); / / Perform an eigen analysis and slice apart the valuesand vectors. var eigens = covarArray.eigen(); / / This is a P-length vector of Eigenvalues. var eigenValues = eigens.slice(1, 0, 1); / / This is a PxP matrix with eigenvectors in rows. var eigenVectors = eigens.slice(1, 1); / / Convert the array image to 2D arrays for matrixcomputations. var arrayImage = arrays.toArray(1) / / Left multiply the image array by the matrixofeigenvectors. var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage); / / Turn the square roots of the Eigenvalues into a P-bandimage. var sdImage = ee.Image(eigenValues.sqrt()) .arrayProject([0]).arrayFlatten([getNewBandNames('SD',bandNames)]); / / Turn the PCs into a P-band image, normalized by SD. return principalComponents / / Throw out an an unneeded dimension, [[]]->[]. .arrayProject([0]) / / Make the one band array image a multi-band image, []->image. .arrayFlatten([getNewBandNames('PC', bandNames)]) / / Normalize the PCs by their SDs. .divide(sdImage) .set('year', year); } / / Get the PCs at the specified scale and in the specified region img = getPrincipalComponents(centered, 1000, region); return img; }; var PCA_imgcol = imgNorcol.map(pca); print(PCA_imgcol); Map.addLayer(PCA_imgcol.first(), {"bands":["PC1"]}, 'pc1'); / / Calculate RSEI using PC1 and perform normalization var RSEI_imgcol = PCA_imgcol.map(function(img){ img = img.addBands(ee.Image(1).rename('constant')) var rsei = img.expression('constant - pc1' , { constant: img.select('constant'), pc1: img.select('PC1') ) rsei = img_normalize(rsei) return img.addBands(rsei.rename('rsei')) ) print(RSEI_imgcol) var visParam = { palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' + '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301' }; Map.addLayer(RSEI_imgcol.first().select('rsei'), visParam, 'rsei') var sign = function(i, j) { / / i and j are images return ee.Image(j).neq(i) / / Zero case .multiply(ee.Image(j).subtract(i).clamp(-1, 1)).int(); }; var kendall = ee.ImageCollection(joined.map(function(current) { var afterCollection = ee.ImageCollection.fromImages(current.get('after')); return afterCollection.map(function(image) { return ee.Image(sign(current, image)).unmask(0); }); }).flatten()).reduce('sum', 2); var palette = ['red', 'white', 'green']; Map.addLayer(kendall, {palette: palette}, 'kendall'); / / Mann-Kendall Test var slope = function(i, j) { return ee.Image(j).subtract(i) .divide(ee.Image(j).date().difference(ee.Image(i).date(), 'days')) .rename('slope') .float(); }; var slopes = ee.ImageCollection(joined.map(function(current) { var afterCollection = ee.ImageCollection.fromImages(current.get('after')); return afterCollection.map(function(image) { return ee.Image(slope(current, image)); }); }).flatten()); var sensSlope = slopes.reduce(ee.Reducer.median(), 2); Map.addLayer(sensSlope, {palette: palette}, 'sensSlope'); / / Sen mutation test / / Parameter extraction and analysis var groups = coll.map(function(i) { var matches = coll.map(function(j) { return i.eq(j); / / i and j are images. }).sum(); return i.multiply(matches.gt(1)); }); / / Compute tie group sizes in a sequence. The first group isdiscarded. var group = function(array) { var length = array.arrayLength(0); / / Array of indices. These are 1-indexed. var indices = ee.Image([1]) .arrayRepeat(0, length) .arrayAccum(0, ee.Reducer.sum()) .toArray(1); var sorted = array.arraySort(); var left = sorted.arraySlice(0, 1); var right = sorted.arraySlice(0, 0, -1); / / Indices of the end of runs. var mask = left.neq(right) / / Always keep the last index, the end of the sequence. .arrayCat(ee.Image(ee.Array([[1]])), 0); var runIndices = indices.arrayMask(mask); / / Subtract the indices to get run lengths. var groupSizes = runIndices.arraySlice(0, 1) .subtract(runIndices.arraySlice(0, 0, -1)); return groupSizes; }; / / See equation 2.6 in Sen (1968). var factors = function(image) { return image.expression('b() * (b() - 1) * (b() * 2 + 5)'); }; var groupSizes = group(groups.toArray()); var groupFactors = factors(groupSizes); var groupFactorSum = groupFactors.arrayReduce('sum', [0]) .arrayGet([0, 0]); var count = joined.count(); var kendallVariance = factors(count) .subtract(groupFactorSum) .divide(18) .float(); Map.addLayer(kendallVariance, {}, 'kendallVariance'); Export.image.toDrive({ image: RSEI_imgcol.first().select('rsei').clip(wuliangsu), / / Cropping the output image using the study area region: wuliangsu, scale: 1000, / / Resolution (meters / pixel) crs: "EPSG:4326", / / Set projection maxPixels: 1e13 / / Limits the number of pixels exported }); / / Export to drive.

[0039] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method for real-time ecological vulnerability assessment of ecological transition zones, wherein the ecological transition zones are located in topographically fragmented areas and / or semi-arid areas, characterized in that: The evaluation method includes the following steps: Step 1: Select a set of ecological vulnerability assessment index factors that can be obtained online in real time and have historical data, and obtain the real-time data and historical data of the ecological vulnerability assessment index factors; Step 2: Calculate the historical ecological vulnerability grading index of the ecological transition zone; Step 3: Conduct long-term time-series analysis of the changing trends and spatial distribution characteristics of ecological vulnerability, and compare whether the calculated abrupt change years match the actual events. If they do not match, repeat steps one to three until the calculated abrupt change years match the actual events. This is considered to have obtained a set of effective ecological vulnerability assessment index factors. Step 4: Using a set of effective ecological vulnerability assessment index factors, programmatically implement the automatic acquisition of current ecological vulnerability assessment index factor data and the real-time calculation of the ecological vulnerability classification index; In step three, the Mann-Kendall mutation test is first used to perform a Mann-Kendall mutation test on the ecological vulnerability grading index to explore the mutation time points and extract the ecological vulnerability grading for important years. Then, the transformation trajectory method is used to represent the changes in ecological vulnerability grading in important years, and the spatial distribution of ecological vulnerability in important years is analyzed intuitively. Finally, Sen+Mann-Kendall trend analysis is used to explore the changing trend of ecological vulnerability and its spatial distribution characteristics.

2. The method for real-time ecological vulnerability assessment of ecological transition zones according to claim 1, characterized in that: Step 3 yields a set of effective ecological vulnerability assessment indicators, including elevation, topographic relief, landscape disturbance, annual average rainfall, atmospheric temperature, surface temperature, surface humidity, surface dryness, population density, vegetation cover, and economic density.

3. The method for real-time ecological vulnerability assessment of ecological transition zones according to claim 1, characterized in that: In step two, the calculation process for the ecological vulnerability grading index is as follows: Each ecological vulnerability assessment index factor was reprojected into a WGS84 ellipsoid and a UTM plane coordinate system, resampled to 1 km, and then subjected to Min-Max normalization. Spatial principal component analysis was used to normalize the first principal component (PC1) into positive indicators for each ecological vulnerability assessment index factor. The ecological vulnerability index is denoted as EVI, and is restored to the range of 0 to 1, EVI = 1 - PC1; Based on the vulnerability classification according to EVI, the ecological vulnerability classification index is obtained, denoted as EVGI.

4. The method for real-time ecological vulnerability assessment of ecological transition zones according to claim 1, characterized in that: Step 2: Calculate the ecological vulnerability grading index for at least 10 years.

5. The method for real-time ecological vulnerability assessment of ecological transition zones according to claim 1, characterized in that: The ecological vulnerability grading indexes for each region calculated in step four are plotted on a color cloud map, with the same ecological vulnerability grading index marked with the same color on the color cloud map.

6. An application of a method for real-time ecological vulnerability assessment in ecological transition zones, characterized in that: Using the real-time ecological vulnerability assessment method for ecological transition zones as described in claim 1, the ecological vulnerability status of ecological transition zones can be monitored in real time. It is recommended to optimize industrial and agricultural production activities in areas where the ecological vulnerability grading index reaches or exceeds level three, and to prioritize ecological restoration in areas where the ecological vulnerability grading index reaches or exceeds level four.

7. The application of the instant ecological vulnerability assessment method for ecological transition zones according to claim 6, characterized in that: The factor detector and interaction detector in the geographic detector were used to analyze the factors of each ecological vulnerability assessment index, explore the driving force of each ecological vulnerability assessment index factor on the overall ecological vulnerability, and prioritize intervention on factors with high single-factor driving force and high interaction synergy when withdrawing industrial and agricultural production activities and carrying out ecological restoration.

Images