Assessment and attribution methods for changes in regional ecological security pattern

Abstract

The present invention relates to an evaluation and attribution method for changes in regional ecological security patterns. The present invention is to solve the problem that the existing technology cannot realize quantitative evaluation and attribution analysis of changes in regional ecological security patterns. The method includes: obtaining various types of geographic spatial data of the ecological region; determining the ecological source according to the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region; obtaining a comprehensive ecological resistance surface by weighting various resistance factors; constructing an ecological security pattern; quantitatively evaluating source connectivity and network connectivity according to ecological sources and ecological corridors; and performing attribution analysis on the relative contributions of climate change and human activities to changes in regional ecological security patterns in the ecological region. The present invention realizes quantitative evaluation and attribution analysis of changes in regional ecological security patterns. The present invention belongs to the field of evaluation and attribution of changes in ecological security patterns.

Description

Technical Field

[0001] The invention relates to an evaluation and attribution method for changes in ecological security pattern. Background Art

[0002] Ecological security pattern is a mode or structure for planning and constructing spatial layout to achieve ecological security. It emphasizes the integrity, connectivity and balance of ecological functions to ensure the continuity and stability of ecosystem services. A reasonable and scientific ecological security pattern can effectively guarantee the health and stability of the ecosystem. Therefore, the ecological security pattern is an important means and way to achieve the goal of ecological security. At present, climate change, unsustainable land use change and rapid urbanization have greatly destroyed the balance of natural ecosystems, resulting in the reduction of regional ecosystem service functions, the increase of ecological sensitive areas, the deterioration of landscape pattern, and the significant changes in the ecological security pattern, which in turn threaten ecological security. Therefore, it is urgent to evaluate the changes in the regional ecological security pattern and quantify the contribution rate of its driving factors. At present, there is a lack of evaluation methods for the effectiveness or rationality of the ecological security pattern, and there are few breakthroughs in the research on the driving factors affecting the changes in the ecological security pattern. Scholars focus on comparing the constructed regional ecological security pattern with the existing land use mode and regional development needs, and qualitatively evaluating the advantages and disadvantages of the current ecological security pattern. Quantitative evaluation research is far from enough, which makes the protection and restoration of the regional ecological security pattern lack of reliable theoretical support. In summary, the problems with existing technologies are: they still lack evaluation methods to accurately assess changes in regional ecological security patterns and to attribute and analyze their driving factors. Summary of the invention

[0003] The purpose of the present invention is to solve the technical problem that the existing technology cannot realize the quantitative evaluation and attribution analysis of the changes in the regional ecological security pattern, and to provide an evaluation and attribution method for the changes in the regional ecological security pattern.

[0004] The evaluation and attribution methods of changes in regional ecological security pattern are as follows:

[0005] S1. Obtain various geospatial data of the ecological region;

[0006] S2. According to various types of geographic spatial data of the ecological region, obtain the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region, and determine the ecological source area according to the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region;

[0007] S3. According to the natural conditions of the ecological region and the resistance factors of social and economic activities, the weights are superimposed to obtain the comprehensive ecological resistance surface;

[0008] S4. According to each ecological source and modified resistance surface, determine the ecological corridors, pinch points, and obstacle points between each ecological source, and build an ecological security pattern;

[0009] S5. Quantitatively assess source connectivity and network connectivity based on ecological sources and ecological corridors;

[0010] S6. Conduct attribution analysis on the relative contributions of climate change and human activities to changes in regional ecological security patterns in the ecological regions;

[0011] The various types of geographic spatial data in step S1 include climate, soil, hydrology, land use and remote sensing image data of the same temporal and spatial scales;

[0012] Step S2 evaluates the three ecosystem service functions of habitat quality, water conservation, and soil conservation based on climate and land use data, and superimposes them with equal weights to obtain an ecological importance index;

[0013] (1) Habitat quality: Based on Formula 1, the habitat quality is calculated by combining the suitability of the habitat, the intensity of external threat sources, and the sensitivity. The calculation formula is as follows:

[0014]

[0015] Where: Q xj represents the habitat quality of pixel x on land cover type j; H j is the habitat suitability of type j land cover;

[0016] (2) Water conservation: Water conservation is obtained by subtracting the total runoff of the previous year from the water production. The water production is estimated based on Formula 2. The total runoff of the previous year is estimated by assigning an average runoff coefficient value to each land use type and multiplying it by the annual precipitation of each pixel in that year.

[0017]

[0018] Where: Y xj is the annual water yield m of grid cell x on land use type j 3 ;AET xj is the annual average evapotranspiration of grid cell x on land use type j; P x is the annual precipitation of grid cell x; AET xj / P x is the Zhang coefficient, which is an algorithm simulated based on the Budyko equation combined with the actual situation;

[0019] Runoff ij =P ij ×C j (3)

[0020] WRij =Y xj -Runoff ij (4)

[0021] In the formula, WR ij It is LUCC j The water conservation capacity of each pixel in the year, Y xj It is LUCC j The water output of each pixel in the year, in mm, Runoff ij It is LUCC j The total annual runoff of each pixel surface, in mm / a, P ij LUCC j The annual precipitation of each pixel in the year, in mm, C j represents the jth type of land use;

[0022] (3) Soil retention: Equation 5 is used to estimate the erosion reduction and sediment retention. Finally, the soil retention is calculated by subtracting the sediment retention from the soil loss.

[0023] SEDRET=PKLS x -USLE x (5)

[0024] In the formula: SEDRET x Represents the soil conservation amount of grid x on a certain land use type; PKLS x represents the potential soil loss of grid x on a certain land use type, and USLEx represents the actual erosion of grid x on a certain land use type, that is, the soil erosion under vegetation cover and soil and water conservation measures;

[0025] The three types of ecosystem service function indicators have the same contribution to ecological importance, so habitat quality, water conservation and soil conservation are normalized. The cumulative ecological importance index score of each grid is the sum of the three types of ecosystem service function scores. The calculation formula is as follows:

[0026] ESI=NorHQ(i,t)+NorWR(i,t)+NorSR(i,t)(6)

[0027] Where: ESI is the ecological importance index score, NorHQ(i, t), NorWR(i, t), and NorSR(i, t) are the normalized values ​​of habitat quality, water conservation, and soil conservation of grid i in year t, respectively;

[0028] Step S2 uses climate, land use and remote sensing image data to calculate three sensitivity indexes: sensitivity to nitrogen and phosphorus pollution, sensitivity to cultivated land soil erosion, and sensitivity to human disturbance, and superimposes them with equal weights to obtain the ecological sensitivity index:

[0029] (1) Nitrogen and phosphorus pollution sensitivity: The nitrogen and phosphorus output of each land use type was estimated based on equations 7, 8, and 9;

[0030] ALV x =HSS x ×POL x (7)

[0031]

[0032] λ x =log(∑ U Y U ) (9)

[0033] Where ALV x is the load value adjusted at grid cell x; POL x is the output coefficient of grid cell x; HSS x is the hydrological sensitivity score at grid cell x, is the average runoff index of the catchment area; x is the runoff index at grid cell x, which is used to express the sum of the water yields of all grids in the catchment area;

[0034] (2) Sensitivity of cultivated land soil erosion: The sensitivity of cultivated land soil erosion is characterized by grain supply, and the net primary productivity (NPP) is combined to achieve a reasonable allocation of grain output on cultivated land in each county of the ecological region. The formula is:

[0035]

[0036] Where: Px is the grain output in the x-th pixel, in kg; NPPx is the NPP value of the x-th pixel on cultivated land; Psum is the total grain output of a county; NPPsum is the sum of the NPP of a county on agricultural land;

[0037] (3) Human interference sensitivity: Based on geographic spatial data such as population density, land use, and road network, population concentration pressure, land use intensity, traffic network pressure, and electricity consumption pressure are estimated; and the cumulative human interference sensitivity index is obtained by adding up the equal weights;

[0038] Population concentration pressure calculation: The population density grid data with a resolution of 1km×1km from the LandScan global population database is used; for grids with more than 1,000 people per square kilometer, the interference intensity score is set to 10 points, and for grids with a population density of 1,000 people / km2 The following grid uses the logarithmic method to calibrate the population density pressure score to obtain the ecological region population concentration pressure index:

[0039]

[0040] Among them, PDI s is the population concentration pressure, pd is the population density, PDI s The value range of is [0, 10], the higher the score, the higher the interference intensity;

[0041] Land use intensity calculation: by assigning values ​​to the relative pressures of various land use forms on natural ecosystems, the land use intensity index of the ecological region is obtained;

[0042] The land use intensity index is constructed by calculating the area ratio of different land use types in the grid to describe the relative size of regional land use risk. The formula is as follows:

[0043]

[0044] Among them, LDI s is the land use intensity; n is the total number of land use types contained in the grid; P i is the percentage of land use type i in the grid area; S i Quantify pressure scores for land-use types;

[0045] Traffic network pressure measurement: Based on the Euclidean distance, the pressure of the traffic network with different buffer radii is assigned to obtain the traffic network pressure index of the ecological region;

[0046] Electricity consumption pressure calculation: convert nighttime light data into electricity consumption to characterize the intensity of electricity consumption; fit the relationship between the ecological area's electricity consumption data and the nighttime light intensity value by constructing a power function fitting model, thereby obtaining the ecological area's electricity consumption pressure;

[0047] E ij =DN sij a (13)

[0048] In the formula, E ij is the total electricity consumption of city j in year i; DN sij is the sum of the DN values ​​of the nighttime light data of city j in year i; a is the estimated parameter of the power function fitting model;

[0049] Cumulative human disturbance intensity calculation: The four types of human disturbance indicators contribute equally to the cumulative human pressure, so the population concentration pressure, land use pressure, traffic network pressure and power consumption pressure are normalized. The cumulative human disturbance pressure score of each grid is the sum of the four types of human pressure scores. The formula is:

[0050] HDI(i, t) = NorPDI s (i,t)+NorLDI s (i, t)+NorTDI s (i,t)+NorEDI s (i,t) (14)

[0051] Where HDI(i,t) is the cumulative human disturbance pressure score of grid i in year t; NorPDI s (i, t), NorLDI s (i, t), NorTDI s (i, t), NorEDI s (i, t) are the normalized values ​​of population concentration pressure, land use pressure, traffic network pressure and electricity consumption pressure scores of grid i in year t;

[0052] The three types of ecological sensitivity indicators contribute equally to ecological sensitivity. Therefore, the sensitivity to nitrogen and phosphorus pollution, sensitivity to cultivated land soil erosion, and sensitivity to human disturbance are normalized. The cumulative ecological sensitivity index score of each grid is the sum of the three types of ecological sensitivity scores. The formula is:

[0053] ESS=NorNDR(i,t)+NorFP(i,t)+NorHDI(i,t) (15)

[0054] Where: ESS is the cumulative ecological sensitivity index, NorNDR(i, t), NorFP(i, t), and NorHDI(i, t) are the normalized values ​​of the nitrogen and phosphorus pollution sensitivity, cultivated land soil erosion sensitivity, and human disturbance sensitivity scores of grid i in year t, respectively;

[0055] The landscape dominance index in step S2 can highly condense the landscape pattern information and reflect the characteristics of its structural composition and spatial configuration. The landscape pattern index method can better measure the regional landscape pattern and reflect the changes in the characteristics of the landscape pattern in a specific region from different scales. The landscape fragmentation index, landscape complexity index, landscape aggregation connectivity index, and landscape heterogeneity index are estimated using land use data, and the landscape dominance index is obtained by superimposing them with equal weights.

[0056] (1) Landscape fragmentation index: The landscape fragmentation index is obtained by the equal-weighted superposition of patch density and average patch area. Patch density reflects the complexity of the landscape spatial structure; the average patch area represents an average condition and can reflect the degree of landscape fragmentation in landscape structure analysis. The formula is:

[0057]

[0058] Where: is the PD patch density, N1 is the number of patches contained in patch type i in the landscape, which is used to express the total number of patches in the landscape; A1 represents the total area of ​​the landscape in each ecological region;

[0059]

[0060] Where: AREA MN is the average patch area, A2 is the total area of ​​a patch type, and N2 is the total number of patches of this patch type;

[0061] LFI=NorPD(i,t)+NorAREA_MN(i,t) (18)

[0062] Where: LFI is the landscape fragmentation index, NorPD(i, t) and NorAREA_MN(i, t) are the normalized values ​​of the patch density and average patch area index score of grid i in year t, respectively;

[0063] (2) Landscape complexity index: The landscape complexity index is obtained by the weighted superposition of edge density and average patch shape. Edge density is an important indicator of patch shape, indicating the degree to which the type is segmented. MSI is used to reflect the overall complexity of the patch shape of the landscape type.

[0064] ED=E / A3(19)

[0065] Where: ED is the edge density, E is the total length of the patch boundary in the study area, and A3 represents the total area of ​​the landscape. The larger the ED value, the greater the edge density, the higher the degree of segmentation of the landscape type, and the more dispersed the layout.

[0066]

[0067] Where: MSI is used to characterize the overall complexity of the patch shape of landscape types, a ij is the area of ​​landscape patch ij, P ij is the perimeter of landscape patch ij, and n is the number of landscape patches;

[0068] When MSI ≥ 1, all plaques are square, MSI = 1, and the larger the value, the more irregular the plaque shape;

[0069] LCI=NorED(i,t)+NorMSI(i,t)(21)

[0070] Where: LCI is the landscape complexity index, NorED(i, t) and NorMSI(i, t) are the normalized values ​​of the edge density and average patch shape index score of grid i in year t, respectively;

[0071] (3) Landscape aggregation connectivity index: The landscape aggregation connectivity index is obtained by the weighted superposition of aggregation and spread. The larger the aggregation value, the greater the degree of spatial aggregation of patch types. The higher the spread value, the better the connectivity of a certain dominant patch type in the landscape. On the contrary, it reflects that the landscape is a dense pattern with multiple elements and the degree of landscape fragmentation is high. The formula is:

[0072]

[0073] Where: AI is the degree of aggregation, g ii is the number of adjacent pixels of the patch type, maxg ii It is the maximum value of the number of adjacent pixels of each patch type, which is used to reflect the spatial aggregation and connectivity of the patch type;

[0074]

[0075] Where: CONTAG is the degree of contagion, P i is the area proportion of patch type i in the landscape, gik is the number of adjacencies between pixels of patch types i and k based on the double counting method, and m is the number of patch types present in the landscape;

[0076] LAFI=NorAI(i,t)+NorCONTAG(i,t) (24)

[0077] Where: LAFI is the landscape aggregation connectivity index, NorAI(i, t) and NorCONTAG(i, t) are the normalized values ​​of the aggregation and spread scores of grid i in year t, respectively;

[0078] (4) Landscape heterogeneity index: The landscape heterogeneity index is characterized by Shannon diversity, which reflects the degree of landscape heterogeneity. The higher the diversity index, the richer the landscape types. When the difference in the proportion of landscape area of ​​each type increases, the diversity decreases. The formula is:

[0079]

[0080] Where: SHDI is Shannon diversity, the ratio of the number of patches to the area in the landscape type is P i It indicates that it reflects the number of landscape elements and the specific proportion of each landscape element;

[0081] The evaluation results of the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region were superimposed with equal weights and divided into five levels according to the natural breakpoint method. The threshold method was used to screen out the areas with an area greater than 200km in the first two levels. 2 The ecological patches are superimposed on the ecological protection areas to obtain the ecological source areas of the ecological areas;

[0082] Step S3 selects seven resistance factor indicators, namely, land use type, altitude, NDVI, night light, distance to road, distance to railway and distance to water, divides the resistance factor level into five categories, and sets the resistance coefficient to 1-5; based on the weight value, the weights are superimposed to obtain the comprehensive ecological resistance surface of the ecological region;

[0083] (1) Land use type

[0084] Land use type affects regional information exchange, material and energy flow; the closer the land use is to the ecological source, the stronger the ecological function is and the smaller the ecological resistance is; on the contrary, the more human interference the land use type is subject to, the greater the ecological resistance is; forest land is the land use type with the least ecological resistance, followed by grassland, water area, cultivated land, unused land and construction land, with resistance values ​​of 1, 2, 3, 4 and 5 respectively;

[0085] (2) Altitude

[0086] Altitude is used to reflect the topography of the ecological area, that is, the higher the altitude, the higher the resistance value; the altitude is divided into five levels: less than 255m, 255m-491m, 491m-748m, 748m-1079m, and greater than 1079m, with resistances of 1, 2, 3, 4, and 5 respectively;

[0087] (3) NDVI

[0088] NDVI is the "normalized difference vegetation index", which is a common index used to detect vegetation growth, vegetation coverage, and dynamic changes in vegetation in time series. NDVI is proportional to vegetation coverage. The larger the NDVI value, the higher the vegetation coverage, the higher the habitat quality, and the smaller the ecological resistance. The natural breakpoint method is used to divide NDVI into five levels: 0-0.296, 0.296-0.384, 0.384-0.469, 0.469-0.563, 0.563-0.804, and the resistance values ​​are 5, 4, 3, 2, and 1 respectively.

[0089] (4) Night lighting

[0090] Night light data is used to reflect the intensity of human activities in the ecological area. That is, the greater the night light data value, the greater the intensity of human interaction, the worse the ecological environment, and the greater the resistance value. The night light is divided into five levels: 0, 0-15, 15-29, 29-47, 47-63, and their resistances are 1, 2, 3, 4, and 5 respectively.

[0091] (5) Distance from the road

[0092] The development of transportation reflects the intensity of human activities. The more developed the transportation is, the more serious the interference to the ecosystem is, the greater the ecological flow resistance is, and the migration and information exchange of organisms are seriously hindered. National roads, provincial roads and expressways are selected as the factors affecting the ecological source, and the distance from the road is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m, and the resistance values ​​of more than 4000m are 5, 4, 3, 2 and 1 respectively.

[0093] (6) Distance from railway

[0094] The distance from the railway is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m, >4000m, and the resistance values ​​are 5, 4, 3, 2, and 1 respectively;

[0095] (7) Distance from water bodies

[0096] Water areas are important pathways for the migration of material energy and species. Therefore, the closer to the water area, the better the ecological quality and the smaller the ecological resistance value. The distance to the water area is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m. If the distance is greater than 4000m, the resistance is 1, 2, 3, 4, 5 respectively.

[0097] Step S4 is based on each ecological source and the comprehensive ecological resistance surface, based on the circuit theory principle, and determines the corridors, pinch points, and obstacle points between each ecological source through Formula 26, and constructs an ecological security pattern:

[0098]

[0099] Among them, I is the size of the ecological flow passing through the corridor; V is the size of the ecological source measured in the entire corridor; R is the cumulative resistance of the corridor;

[0100] Step S5 quantitatively evaluates source connectivity and network connectivity based on ecological sources and corridors in the ecological region;

[0101] (1) Evaluation of source connectivity Based on the ecological source, the connectivity of the ecological source is analyzed using formulas 27 and 28 to calculate its connectivity probability and connectivity comprehensive index:

[0102]

[0103] Among them, PC is the connectivity probability, IIC is the comprehensive connectivity index, n is the total number of ecological sources in the landscape, and a i and a j is the property of ecological source i and ecological source, nl ij is the number of ecological corridors in the shortest path between patches i and j, is the maximum value of the path probability product between ecological source i and ecological source j, A L is the total area of ​​the ecoregion in question;

[0104] (2) Network connectivity analysis is based on ecological source areas and ecological corridors, which are simplified into network structure indexes of abstract “point-line” graphs. Network closure, network point-line ratio, and network connection ratio are selected as indicators for evaluating the connectivity of the ecological network of the ecological region through formulas 29, 30, and 31:

[0105] Network Closure: Measures the closure of a network, which is the ratio of the actual number of loops in the network to the maximum possible number of loops that exist in the network. It is the degree to which the loops of landscape nodes in the network are connected:

[0106]

[0107] Where: α represents the network closure, L represents the number of corridors, V is the number of nodes, L-V+1 is the actual number of loops, and 2V-5 is the maximum possible number of loops;

[0108] Network point-line ratio: It measures the average number of links per node in the network and is an indicator of how easy it is for a landscape node to connect with other landscape nodes:

[0109]

[0110] In the formula: β represents the network point-line rate, L represents the number of corridors, and V is the number of nodes;

[0111] Network connectivity: measures the connectivity of a network, describing the degree to which all landscape nodes in a network are connected:

[0112]

[0113] In the formula: γ represents the network connectivity rate, L represents the number of corridors, and V represents the number of nodes;

[0114] Step S6 quantifies the relative contributions of climate change and human activities to changes in regional ecological security patterns through attribution analysis;

[0115] The independent impact of climate change or human activities on the ecological security pattern was estimated by calculating the difference between the baseline scenario and the climate scenario or human activity scenario, and the interactive effect of climate change and human activities was calculated by subtracting the independent impact of climate change and human activities from the combined effect of climate change and human activities. The contribution of the combined effect, climate change effect, human activity effect and interactive effect to the connectivity of the ecological security pattern of the described eco-region was quantified, thereby achieving attribution analysis of the changes in the ecological security pattern.

[0116] Calculate the relative contribution η a ,η b , η c , η d The expressions are:

[0117]

[0118] Among them, η a ,η b , η c , η d Respectively represent the relative contribution rates of interaction effects, human activities, climate change and combined effects, Δr a ,Δr b , Δr c , Δr d They represent the change values ​​of ecological security pattern connectivity caused by interaction effects, human activities, climate change and comprehensive effects respectively, and R is the ecological security pattern connectivity value in the baseline scenario.

[0119] The present invention quantifies the stability of the regional ecological security pattern based on source connectivity and network connectivity, and quantifies the impact of climate change and human activities on changes in the regional ecological security pattern based on different scenarios of climate change and human activities, thus achieving quantitative assessment and attribution analysis of changes in the regional ecological security pattern. BRIEF DESCRIPTION OF THE DRAWINGS

[0120] Figure 1 It is a schematic diagram of the overall process of the method for assessing and attributing changes in the regional ecological security pattern of the present invention;

[0121] Figure 2 This is the trend chart of connectivity probability (PC) and comprehensive connectivity index (IIC) in the source connectivity analysis proposed in Experiment 1;

[0122] Figure 3 This is the trend chart of the network closure (α index), network point-line ratio (β index), and network connection ratio (γ index) in the network connectivity analysis proposed in Experiment 1;

[0123] Figure 4 This is the trend chart of the contribution rate of climate change and human activities to changes in the ecological security pattern proposed in Experiment 1. DETAILED DESCRIPTION

[0124] The technical solution of the present invention is not limited to the specific implementation modes listed below, but also includes any combination of the specific implementation modes.

[0125] Specific implementation method 1: The evaluation and attribution method of the change of regional ecological security pattern in this implementation method is as follows:

[0126] S1. Obtain various geospatial data of the ecological region;

[0127] S2. According to various types of geographic spatial data of the ecological region, obtain the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region, and determine the ecological source area according to the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region;

[0128] S3. According to the natural conditions of the ecological region and the resistance factors of social and economic activities, the weights are superimposed to obtain the comprehensive ecological resistance surface;

[0129] S4. According to each ecological source and modified resistance surface, determine the ecological corridors, pinch points, and obstacle points between each ecological source, and build an ecological security pattern;

[0130] S5. Quantitatively assess source connectivity and network connectivity based on ecological sources and ecological corridors;

[0131] S6. Conduct attribution analysis on the relative contributions of climate change and human activities to changes in regional ecological security patterns in the ecological regions;

[0132] The various types of geographic spatial data in step S1 include climate, soil, hydrology, land use and remote sensing image data of the same temporal and spatial scales;

[0133] Step S2 evaluates the three ecosystem service functions of habitat quality, water conservation, and soil conservation based on climate and land use data, and superimposes them with equal weights to obtain an ecological importance index;

[0134] (1) Habitat quality: Based on Formula 1, the habitat quality is calculated by combining the suitability of the habitat, the intensity of external threat sources, and the sensitivity. The calculation formula is as follows:

[0135]

[0136] Where: Q xj represents the habitat quality of pixel x on land cover type j; H j is the habitat suitability of type j land cover;

[0137] (2) Water conservation: Water conservation is obtained by subtracting the total runoff of the previous year from the water production. The water production is estimated based on Formula 2. The total runoff of the previous year is estimated by assigning an average runoff coefficient value to each land use type and multiplying it by the annual precipitation of each pixel in that year.

[0138]

[0139] Where: Y xj is the annual water yield m of grid cell x on land use type j 3 ;AET xj is the annual average evapotranspiration of grid cell x on land use type j; P x is the annual precipitation of grid cell x; AET xj / P x is the Zhang coefficient, which is an algorithm simulated based on the Budyko equation combined with the actual situation;

[0140] Runoff ij =P ij ×C j (3)

[0141] WR ij =Y xj -Runoff ij (4)

[0142] In the formula, WR ij It is LUCC j The water conservation capacity of each pixel in the year, Y xj It is LUCC j The water output of each pixel in the year, in mm, Runoff ij It is LUCC j The total annual runoff of each pixel surface, in mm / a, P ij LUCC j The annual precipitation of each pixel in the year, in mm, C j represents the jth type of land use;

[0143] (3) Soil retention: Equation 5 is used to estimate the erosion reduction and sediment retention. Finally, the soil retention is calculated by subtracting the sediment retention from the soil loss.

[0144] SEDRET=PKLS x -USLE x (5)

[0145] In the formula: SEDRET x Represents the soil conservation amount of grid x on a certain land use type; PKLS xrepresents the potential soil loss of grid x on a certain land use type, and USLEx represents the actual erosion of grid x on a certain land use type, that is, the soil erosion under vegetation cover and soil and water conservation measures;

[0146] The three types of ecosystem service function indicators have the same contribution to ecological importance, so habitat quality, water conservation and soil conservation are normalized. The cumulative ecological importance index score of each grid is the sum of the three types of ecosystem service function scores. The calculation formula is as follows:

[0147] ESI=NorHQ(i,t)+NorWR(i,t)+NorSR(i,t) (6)

[0148] Where: ESI is the ecological importance index score, NorHQ(i, t), NorWR(i, t), and NorSR(i, t) are the normalized values ​​of habitat quality, water conservation, and soil conservation of grid i in year t, respectively;

[0149] Step S2 uses climate, land use and remote sensing image data to calculate three sensitivity indexes: sensitivity to nitrogen and phosphorus pollution, sensitivity to cultivated land soil erosion, and sensitivity to human disturbance, and superimposes them with equal weights to obtain the ecological sensitivity index:

[0150] (1) Nitrogen and phosphorus pollution sensitivity: The nitrogen and phosphorus output of each land use type was estimated based on equations 7, 8, and 9;

[0151] ALV x =HSS x ×POL x i7)

[0152]

[0153] λ x =log(Z U Y U ) (9)

[0154] Where ALV x is the load value adjusted at grid cell x; POL x is the output coefficient of grid cell x; HSS x is the hydrological sensitivity score at grid cell x, is the average runoff index of the catchment area; x is the runoff index at grid cell x, which is used to express the sum of the water yields of all grids in the catchment area;

[0155] (2) Sensitivity of cultivated land soil erosion: The sensitivity of cultivated land soil erosion is characterized by grain supply, and the net primary productivity (NPP) is combined to achieve a reasonable allocation of grain output on cultivated land in each county of the ecological region. The formula is:

[0156]

[0157] Where: Px is the grain output in the x-th pixel, in kg; NPPx is the NPP value of the x-th pixel on cultivated land; Psum is the total grain output of a county; NPPsum is the sum of the NPP of a county on agricultural land;

[0158] (3) Human interference sensitivity: Based on geographic spatial data such as population density, land use, and road network, population concentration pressure, land use intensity, traffic network pressure, and electricity consumption pressure are estimated; and the cumulative human interference sensitivity index is obtained by adding up the equal weights;

[0159] Population concentration pressure calculation: The population density grid data with a resolution of 1km×1km from the LandScan global population database is used; for grids with more than 1,000 people per square kilometer, the interference intensity score is set to 10 points, and for grids with a population density of 1,000 people / km 2 The following grid uses the logarithmic method to calibrate the population density pressure score to obtain the ecological region population concentration pressure index:

[0160]

[0161] Among them, PDI s is the population concentration pressure, pd is the population density, PDI s The value range of is [0, 10], the higher the score, the higher the interference intensity;

[0162] Land use intensity calculation: by assigning values ​​to the relative pressures of various land use forms on natural ecosystems, the land use intensity index of the ecological region is obtained;

[0163] The land use intensity index is constructed by calculating the area ratio of different land use types in the grid to describe the relative size of regional land use risk. The formula is as follows:

[0164]

[0165] Among them, LDI s is the land use intensity; n is the total number of land use types contained in the grid; P i is the percentage of land use type i in the grid area; S i Quantify pressure scores for land-use types;

[0166] Traffic network pressure measurement: Based on the Euclidean distance, the pressure of the traffic network with different buffer radii is assigned to obtain the traffic network pressure index of the ecological region;

[0167] Electricity consumption pressure calculation: convert nighttime light data into electricity consumption to characterize the intensity of electricity consumption; fit the relationship between the ecological area's electricity consumption data and the nighttime light intensity value by constructing a power function fitting model, thereby obtaining the ecological area's electricity consumption pressure;

[0168] E ij =DN sij a (13)

[0169] In the formula, E ij is the total electricity consumption of city j in year i; DN sij is the sum of the DN values ​​of the nighttime light data of city j in year i; a is the estimated parameter of the power function fitting model;

[0170] Cumulative human disturbance intensity calculation: The four types of human disturbance indicators contribute equally to the cumulative human pressure, so the population concentration pressure, land use pressure, traffic network pressure and power consumption pressure are normalized. The cumulative human disturbance pressure score of each grid is the sum of the four types of human pressure scores. The formula is:

[0171] HDI(i, t) = NorPDI s (i,t)+NorLDI s (i, t)+NorTDI s (i,t)+NorEDI s (i,t) (14)

[0172] Wherein, HDI(i, t) is the cumulative human disturbance pressure score of grid i in year t; NorPDI s (i, t), NorLDI s (i, t), NorTDI s (i, t), NorEDI s (i, t) are the normalized values ​​of population concentration pressure, land use pressure, traffic network pressure and electricity consumption pressure scores of grid i in year t;

[0173] The three types of ecological sensitivity indicators contribute equally to ecological sensitivity. Therefore, the sensitivity to nitrogen and phosphorus pollution, sensitivity to cultivated land soil erosion, and sensitivity to human disturbance are normalized. The cumulative ecological sensitivity index score of each grid is the sum of the three types of ecological sensitivity scores. The formula is:

[0174] ESS=NorNDR(i,t)+NorFP(i,t)+NorHDI(i,t) (15)

[0175] Where: ESS is the cumulative ecological sensitivity index, NorNDR(i, t), NorFP(i, t), and NorHDI(i, t) are the normalized values ​​of the nitrogen and phosphorus pollution sensitivity, cultivated land soil erosion sensitivity, and human disturbance sensitivity scores of grid i in year t, respectively;

[0176] The landscape dominance index in step S2 can highly condense the landscape pattern information and reflect the characteristics of its structural composition and spatial configuration. The landscape pattern index method can better measure the regional landscape pattern and reflect the changes in the characteristics of the landscape pattern in a specific region from different scales. The landscape fragmentation index, landscape complexity index, landscape aggregation connectivity index, and landscape heterogeneity index are estimated using land use data, and the landscape dominance index is obtained by superimposing them with equal weights.

[0177] (1) Landscape fragmentation index: The landscape fragmentation index is obtained by the equal-weighted superposition of patch density and average patch area. Patch density reflects the complexity of the landscape spatial structure; the average patch area represents an average condition and can reflect the degree of landscape fragmentation in landscape structure analysis. The formula is:

[0178]

[0179] Where: is the PD patch density, N1 is the number of patches contained in patch type i in the landscape, which is used to express the total number of patches in the landscape; A1 represents the total area of ​​the landscape in each ecological region;

[0180]

[0181] Where: AREA MN is the average patch area, A2 is the total area of ​​a patch type, and N2 is the total number of patches of this patch type;

[0182] LFI=NorPD(i,t)+NorAREA_MN(i,t) (18)

[0183] Where: LFI is the landscape fragmentation index, NorPD(i, t) and NorAREA_MN(i, t) are the normalized values ​​of the patch density and average patch area index score of grid i in year t, respectively;

[0184] (2) Landscape complexity index: The landscape complexity index is obtained by the weighted superposition of edge density and average patch shape. Edge density is an important indicator of patch shape, indicating the degree to which the type is segmented. MSI is used to reflect the overall complexity of the patch shape of the landscape type.

[0185] ED=E / A3 (19)

[0186] Where: ED is the edge density, E is the total length of the patch boundary in the study area, and A3 represents the total area of ​​the landscape. The larger the ED value, the greater the edge density, the higher the degree of segmentation of the landscape type, and the more dispersed the layout.

[0187]

[0188] Where: MSI is used to characterize the overall complexity of the patch shape of landscape types, a ij is the area of ​​landscape patch ij, P ij is the perimeter of landscape patch ij, and n is the number of landscape patches;

[0189] When MSI ≥ 1, all plaques are square, MSI = 1, and the larger the value, the more irregular the plaque shape;

[0190] LCI=NorED(i,t)+NorMSI(i,t) (21)

[0191] Where: LCI is the landscape complexity index, NorED(i, t) and NorMSI(i, t) are the normalized values ​​of the edge density and average patch shape index score of grid i in year t, respectively;

[0192] (3) Landscape aggregation connectivity index: The landscape aggregation connectivity index is obtained by the weighted superposition of aggregation and spread. The larger the aggregation value, the greater the degree of spatial aggregation of patch types. The higher the spread value, the better the connectivity of a certain dominant patch type in the landscape. On the contrary, it reflects that the landscape is a dense pattern with multiple elements and the degree of landscape fragmentation is high. The formula is:

[0193]

[0194] Where: AI is the aggregation, gii is the number of adjacent pixels of the patch type, maxgii is the maximum value of the number of adjacent pixels of the patch type, which is used to reflect the spatial aggregation and connectivity of the patch type;

[0195]

[0196] Where: CONTAG is the degree of contagion, P i is the area proportion of patch type i in the landscape, gik is the number of adjacencies between pixels of patch types i and k based on the double counting method, and m is the number of patch types present in the landscape;

[0197] LAFI=NorAI(i,t)+NorCONTAG(i,t) (24)

[0198] Where: LAFI is the landscape aggregation connectivity index, NorAI(i,t) and NorCONTAG(i,t) are the normalized values ​​of the aggregation and spread scores of grid i in year t, respectively;

[0199] (4) Landscape heterogeneity index: The landscape heterogeneity index is characterized by Shannon diversity, which reflects the degree of landscape heterogeneity. The higher the diversity index, the richer the landscape types. When the difference in the proportion of landscape area of ​​each type increases, the diversity decreases. The formula is:

[0200]

[0201] Where: SHDI is Shannon diversity, the ratio of the number of patches to the area in the landscape type is P i It indicates that it reflects the number of landscape elements and the specific proportion of each landscape element;

[0202] The evaluation results of the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region were superimposed with equal weights and divided into five levels according to the natural breakpoint method. The threshold method was used to screen out the areas with an area greater than 200km in the first two levels. 2 The ecological patches are superimposed on the ecological protection areas to obtain the ecological source areas of the ecological areas;

[0203] Step S3 selects seven resistance factor indicators, namely, land use type, altitude, NDVI, night light, distance to road, distance to railway and distance to water, divides the resistance factor level into five categories, and sets the resistance coefficient to 1-5; based on the weight value, the weights are superimposed to obtain the comprehensive ecological resistance surface of the ecological region;

[0204] (1) Land use type

[0205] Land use type affects regional information exchange, material and energy flow; the closer the land use is to the ecological source, the stronger the ecological function is and the smaller the ecological resistance is; on the contrary, the more human interference the land use type is subject to, the greater the ecological resistance is; forest land is the land use type with the least ecological resistance, followed by grassland, water area, cultivated land, unused land and construction land, with resistance values ​​of 1, 2, 3, 4 and 5 respectively;

[0206] (2) Altitude

[0207] Altitude is used to reflect the topography of the ecological area, that is, the higher the altitude, the higher the resistance value; the altitude is divided into five levels: less than 255m, 255m-491m, 491m-748m, 748m-1079m, and greater than 1079m, with resistances of 1, 2, 3, 4, and 5 respectively;

[0208] (3) NDVI

[0209] NDVI is the "normalized difference vegetation index", which is a common index used to detect vegetation growth, vegetation coverage, and dynamic changes in vegetation in time series. NDVI is proportional to vegetation coverage. The larger the NDVI value, the higher the vegetation coverage, the higher the habitat quality, and the smaller the ecological resistance. The natural breakpoint method is used to divide NDVI into five levels: 0-0.296, 0.296-0.384, 0.384-0.469, 0.469-0.563, 0.563-0.804, and the resistance values ​​are 5, 4, 3, 2, and 1 respectively.

[0210] (4) Night lighting

[0211] Night light data is used to reflect the intensity of human activities in the ecological area. That is, the greater the night light data value, the greater the intensity of human interaction, the worse the ecological environment, and the greater the resistance value. The night light is divided into five levels: 0, 0-15, 15-29, 29-47, 47-63, and their resistances are 1, 2, 3, 4, and 5 respectively.

[0212] (5) Distance from the road

[0213] The development of transportation reflects the intensity of human activities. The more developed the transportation is, the more serious the interference to the ecosystem is, the greater the ecological flow resistance is, and the migration and information exchange of organisms are seriously hindered. National roads, provincial roads and expressways are selected as the factors affecting the ecological source, and the distance from the road is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m, and the resistance values ​​of more than 4000m are 5, 4, 3, 2 and 1 respectively.

[0214] (6) Distance from railway

[0215] The distance from the railway is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m, >4000m, and the resistance values ​​are 5, 4, 3, 2, and 1 respectively;

[0216] (7) Distance from water bodies

[0217] Water areas are important pathways for the migration of material energy and species. Therefore, the closer to the water area, the better the ecological quality and the smaller the ecological resistance value. The distance to the water area is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m. If the distance is greater than 4000m, the resistance is 1, 2, 3, 4, 5 respectively.

[0218] Step S4 is based on each ecological source and the comprehensive ecological resistance surface, based on the circuit theory principle, and determines the corridors, pinch points, and obstacle points between each ecological source through Formula 26, and constructs an ecological security pattern:

[0219]

[0220] Among them, I is the size of the ecological flow passing through the corridor; V is the size of the ecological source measured in the entire corridor; R is the cumulative resistance of the corridor;

[0221] Step S5 quantitatively evaluates source connectivity and network connectivity based on ecological sources and corridors in the ecological region;

[0222] (1) Evaluation of source connectivity Based on the ecological source, the connectivity of the ecological source is analyzed using formulas 27 and 28 to calculate its connectivity probability and connectivity comprehensive index:

[0223]

[0224] Among them, PC is the connectivity probability, IIC is the comprehensive connectivity index, n is the total number of ecological sources in the landscape, and a i and a j is the property of ecological source i and ecological source, nl ij is the number of ecological corridors in the shortest path between patches i and j, is the maximum value of the path probability product between ecological source i and ecological source j, A L is the total area of ​​the ecoregion in question;

[0225] (2) Network connectivity analysis is based on ecological source areas and ecological corridors, which are simplified into network structure indexes of abstract “point-line” graphs. Network closure, network point-line ratio, and network connection ratio are selected as indicators for evaluating the connectivity of the ecological network of the ecological region through formulas 29, 30, and 31:

[0226] Network Closure: Measures the closure of a network, which is the ratio of the actual number of loops in the network to the maximum possible number of loops that exist in the network. It is the degree to which the loops of landscape nodes in the network are connected:

[0227]

[0228] Where: α represents the network closure, L represents the number of corridors, V is the number of nodes, L-V+1 is the actual number of loops, and 2V-5 is the maximum possible number of loops;

[0229] Network point-line ratio: It measures the average number of links per node in the network and is an indicator of how easy it is for a landscape node to connect with other landscape nodes:

[0230]

[0231] In the formula: β represents the network point-line rate, L represents the number of corridors, and V is the number of nodes;

[0232] Network connectivity: measures the connectivity of a network, describing the degree to which all landscape nodes in a network are connected:

[0233]

[0234] In the formula: γ represents the network connectivity rate, L represents the number of corridors, and V represents the number of nodes;

[0235] Step S6 quantifies the relative contributions of climate change and human activities to changes in regional ecological security patterns through attribution analysis;

[0236] The independent impact of climate change or human activities on the ecological security pattern was estimated by calculating the difference between the baseline scenario and the climate scenario or human activity scenario, and the interactive effect of climate change and human activities was calculated by subtracting the independent impact of climate change and human activities from the combined effect of climate change and human activities. The contribution of the combined effect, climate change effect, human activity effect and interactive effect to the connectivity of the ecological security pattern of the described eco-region was quantified, thereby achieving attribution analysis of the changes in the ecological security pattern.

[0237] Calculate the relative contribution η a ,η b , η c , η d The expressions are:

[0238]

[0239] Among them, η a ,η b , η c , η d Respectively represent the relative contribution rates of interaction effects, human activities, climate change and combined effects, Δr a ,Δr b , Δr c , Δr d They represent the change values ​​of ecological security pattern connectivity caused by interaction effects, human activities, climate change and comprehensive effects respectively, and R is the ecological security pattern connectivity value in the baseline scenario.

[0240] The following experiments were used to verify the effect of the present invention:

[0241] Experiment 1:

[0242] The evaluation and attribution methods of changes in regional ecological security pattern are as follows:

[0243] S1. Obtain various geospatial data of the ecological region;

[0244] S2. Based on various types of geographic spatial data of the ecological region, obtain the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region, and determine the ecological source area based on the three types of indexes;

[0245] S3. According to the natural conditions of the ecological region and the resistance factors of social and economic activities, the weights are superimposed to obtain the comprehensive ecological resistance surface;

[0246] S4. According to each ecological source and modified resistance surface, determine the ecological corridors, pinch points, and obstacle points between each ecological source, and build an ecological security pattern;

[0247] S5. Quantitatively assess source connectivity and network connectivity based on ecological sources and ecological corridors;

[0248] S6. Conduct attribution analysis on the relative contributions of climate change and human activities to changes in regional ecological security patterns in the ecological regions;

[0249] The various types of geographic spatial data in step S1 include climate, soil, hydrology, land use and remote sensing image data of the same temporal and spatial scales;

[0250] Step S2 evaluates the three ecosystem service functions of habitat quality, water conservation, and soil conservation based on climate and land use data, and superimposes them with equal weights to obtain an ecological importance index;

[0251] (1) Habitat quality (HQ): Based on formula 1, the habitat quality is calculated by combining the suitability of the habitat, the intensity of external threat sources, and the sensitivity. The calculation formula is as follows:

[0252]

[0253] Where: Q xj represents the habitat quality of pixel x on land cover type j; H j is the habitat suitability of type j land cover;

[0254] (2) Water conservation (WR): Water conservation is obtained by subtracting the total runoff of the previous year from the water production. The water production is estimated based on Formula 2. The total runoff is estimated by assigning an average runoff coefficient value to each land use type and multiplying it by the annual precipitation of each pixel in that year.

[0255]

[0256] Where: Y xj is the annual water yield of grid cell x on land use type j (m 3 ); AET xjis the annual average evapotranspiration of grid cell x on land use type j; P x is the annual precipitation of grid cell x; AET xj / P x is the Zhang coefficient, which is an algorithm simulated based on the Budyko equation combined with the actual situation;

[0257] Runoff ij =P ij ×C ij (3)

[0258] WR ij =Y xj -Runoff ij (4)

[0259] In the formula, WR ij It is LUCC j The water conservation capacity of each pixel in the year, Y xj It is LUCC j The water output of each pixel in the year, in mm, Runoff ij It is LUCC j The total annual runoff of each pixel surface, in mm / a, P ij LUCC j The annual precipitation of each pixel in the year, in mm, C j It represents the surface runoff coefficient of the jth land use type as shown in Table 1;

[0260] Table 1 Average runoff coefficient data

[0261]

[0262]

[0263] (3) Soil retention (SR): Equation 5 is used to estimate the erosion reduction and sediment retention. Finally, the soil retention is calculated by subtracting the sediment retention from the soil loss.

[0264] SEDRET=PKLS x -USLE x (5)

[0265] In the formula: SEDRET x Represents the soil conservation amount of grid x on a certain land use type; PKLS x represents the potential soil loss of grid x on a certain land use type, and USLEx represents the actual erosion of grid x on a certain land use type, that is, the soil erosion under vegetation cover and soil and water conservation measures;

[0266] The three types of ecosystem service function indicators have the same contribution to ecological importance, so habitat quality, water conservation and soil conservation are normalized. The cumulative ecological importance index score (ESI) of each grid is the sum of the three types of ecosystem service function scores. The calculation formula is as follows:

[0267] ESI=NorHQ(i,t)+NorWR(i,t)+NorSR(i,t) (6)

[0268] Where: ESI is the ecological importance index score, NorHQ(i, t), NorWR(i, t), and NorSR(i, t) are the normalized values ​​of habitat quality, water conservation, and soil conservation of grid i in year t, respectively;

[0269] Step S2 uses climate, land use and remote sensing image data to calculate three sensitivity indexes: sensitivity to nitrogen and phosphorus pollution, sensitivity to cultivated land soil erosion, and sensitivity to human disturbance, and superimposes them with equal weights to obtain the ecological sensitivity index:

[0270] (1) Nitrogen and phosphorus pollution sensitivity (NDR): The nitrogen and phosphorus output of each land use type was estimated based on equations 7, 8, and 9;

[0271] ALV x =HSS x ×POL x (7)

[0272]

[0273] λ x =log(Σ U Y U ) (9)

[0274] Where ALV x is the load value adjusted at grid cell x; POL x is the output coefficient of grid cell x; HSS x is the hydrological sensitivity score at grid cell x, is the average runoff index of the catchment area; x is the runoff index at grid cell x, which is used to express the sum of the water yields of all grids in the catchment area;

[0275] (2) Sensitivity of cultivated land to soil erosion (FP): The sensitivity of cultivated land to soil erosion is characterized by grain supply, and the net primary productivity (NPP) is combined to achieve a reasonable distribution of grain output on cultivated land in each county of the ecological region. The formula is:

[0276]

[0277] Where: P x is the grain yield in the x-th pixel, in kg; NPP x is the NPP value of the x-th pixel on cultivated land; Psum is the total grain output of a certain county; NPPsum is the sum of the NPP of a certain county on agricultural land;

[0278] (3) Human Disturbance Sensitivity (HDI): Based on geographic spatial data such as population density, land use, and road network, population concentration pressure, land use intensity, traffic network pressure, and electricity consumption pressure are estimated; and their weighted superposition is used to obtain the cumulative human disturbance sensitivity index;

[0279] Population concentration pressure calculation: The population density grid data with a resolution of 1km×1km from the LandScan global population database is used; for grids with more than 1,000 people per square kilometer, the interference intensity score is set to 10 points, and for grids with a population density of 1,000 people / km 2 The following grid uses the logarithmic method to calibrate the population density pressure score to obtain the ecological region population concentration pressure index:

[0280]

[0281] Among them, PDI s is the population concentration pressure, pd is the population density, PDI s The value range of is [0, 10], the higher the score, the higher the interference intensity;

[0282] Land use intensity calculation: By assigning values ​​to the relative pressure of each land use form on the natural ecosystem, the land use intensity index of the ecological region is obtained (Table 2);

[0283] Table 2 Quantitative assignment of pressure on land use types

[0284]

[0285] The land use intensity index is constructed by calculating the area ratio of different land use types in the grid to describe the relative size of regional land use risk. The formula is as follows:

[0286]

[0287] Among them, LDI s is the land use intensity; n is the total number of land use types contained in the grid; P i is the percentage of land use type i in the grid area; S i Quantify pressure scores for land-use types;

[0288] Traffic network pressure measurement: Based on the Euclidean distance, the pressure of the traffic network with different buffer radii is assigned to obtain the traffic network pressure index of the ecological region (Table 3);

[0289] Table 3 Traffic pressure assignment based on Euclidean distance

[0290]

[0291] Electricity consumption pressure calculation: convert nighttime light data into electricity consumption to characterize the intensity of electricity consumption; fit the relationship between the ecological area's electricity consumption data and the nighttime light intensity value by constructing a power function fitting model, thereby obtaining the ecological area's electricity consumption pressure;

[0292] E ij =DN sij a (13)

[0293] In the formula, E ij is the total electricity consumption of city j in year i; DN sij ; is the sum of the DN values ​​of the nighttime light data of city j in year i; a is the estimated parameter of the power function fitting model;

[0294] Cumulative human disturbance intensity calculation: The four types of human disturbance indicators contribute equally to the cumulative human pressure, so the population concentration pressure, land use pressure, traffic network pressure and power consumption pressure are normalized. The cumulative human disturbance pressure score of each grid is the sum of the four types of human pressure scores. The formula is:

[0295] HDI(i, t) = NorPDI s (i,t)+NorLDI s (i, t)+NorTDI s (i,t)+NorEDI s (i,t) (14)

[0296] Wherein, HDI(i, t) is the cumulative human disturbance pressure score of grid i in year t; NorPDI s (i,t), NorLDI s (i, t), NorTDI s (i, t), NorEDI s (i, t) are the normalized values ​​of population concentration pressure, land use pressure, traffic network pressure and electricity consumption pressure scores of grid i in year t;

[0297] The three types of ecological sensitivity indicators contribute equally to ecological sensitivity. Therefore, the sensitivity to nitrogen and phosphorus pollution, sensitivity to cultivated land soil erosion, and sensitivity to human disturbance are normalized. The cumulative ecological sensitivity index score of each grid is the sum of the three types of ecological sensitivity scores. The formula is:

[0298] ESS=NorNDR(i,t)+NorFP(i,t)+NorHDI(i,t) (15)

[0299] Where: ESS is the cumulative ecological sensitivity index, NorNDR(i, t), NorFP(i, t), and NorHDI(i, t) are the normalized values ​​of the nitrogen and phosphorus pollution sensitivity, cultivated land soil erosion sensitivity, and human disturbance sensitivity scores of grid i in year t, respectively;

[0300] The landscape dominance index in step S2 can highly condense the landscape pattern information and reflect the characteristics of its structural composition and spatial configuration. The landscape pattern index method can better measure the regional landscape pattern and reflect the changes in the characteristics of the landscape pattern in a specific region from different scales. The landscape fragmentation index, landscape complexity, landscape aggregation connectivity index, and landscape heterogeneity index are estimated using land use data, and the landscape dominance index is obtained by superimposing them with equal weights.

[0301] (1) Landscape fragmentation index (LFI): The landscape fragmentation index is composed of patch density (PD) and average patch area (AREA MN ) are obtained by weighted superposition; patch density reflects the complexity of landscape spatial structure; the average patch area represents an average state, which can reflect the degree of landscape fragmentation in landscape structure analysis. The formula is:

[0302]

[0303] Where: is the PD patch density, N1 is the number of patches contained in patch type i in the landscape, which is used to express the total number of patches in the landscape; A1 represents the total area of ​​the landscape in each ecological region;

[0304]

[0305] Where: AREA MN is the average patch area, A2 is the total area of ​​a patch type, and N2 is the total number of patches of this patch type;

[0306] LFI=NorPD(i,t)+NorAREA_MN(i,t) (18)

[0307] Where: LFI is the landscape fragmentation index, NorPD(i, t) and NorAREA_MN(i, t) are the normalized values ​​of the patch density and average patch area index score of grid i in year t, respectively;

[0308] (2) Landscape complexity index (LCI): The landscape complexity index is obtained by weighted superposition of edge density (ED) and mean patch shape (MSI). Edge density is an important indicator of patch shape, indicating the degree of segmentation of the type. MSI is used to reflect the overall complexity of the patch shape of the landscape type.

[0309] ED=E / A3 (19)

[0310] Where ED is the edge density, E is the total length of the patch boundary in the study area, and A3 represents the total landscape area.

[0311] The larger the ED value, the greater the edge density, the higher the degree of segmentation of the landscape type, and the more dispersed the layout;

[0312]

[0313] Where: MSI is used to characterize the overall complexity of the patch shape of landscape types, a ij is the area of ​​landscape patch ij, P ij is the perimeter of landscape patch ij, and n is the number of landscape patches. When MSI ≥ 1, all patches are square, MSI = 1, and the larger the value, the more irregular the patch shape;

[0314] LCI=NorED(i,t)+NorMSI(i,t) (21)

[0315] Where: LCI is the landscape complexity index, NoRED(i, t) and NorMSI(i, t) are the normalized values ​​of the edge density and average patch shape index score of grid i in year t, respectively;

[0316] (3) Landscape aggregation connectivity index (LAFI): The landscape aggregation connectivity index is obtained by the weighted superposition of aggregation (AI) and contagion (CONTAG). The larger the aggregation value, the greater the degree of spatial aggregation of patch types. The higher the contagion value, the better the connectivity of a certain dominant patch type in the landscape. On the contrary, it reflects that the landscape is a dense pattern with multiple elements and the degree of landscape fragmentation is high. The formula is:

[0317]

[0318] Where: AI is the degree of aggregation, g ii is the number of adjacent pixels of the patch type, maxg iiIt is the maximum value of the number of adjacent pixels of each patch type, which is used to reflect the spatial aggregation and connectivity of the patch type;

[0319]

[0320] Where: CONTAG is the degree of contagion, P i is the area proportion of patch type i in the landscape, gik is the number of adjacencies between pixels of patch types i and k based on the double counting method, and m is the number of patch types present in the landscape;

[0321] LAFI=NorAI(i,t)+NorCONTAG(i,t) (24)

[0322] Where: LAFI is the landscape aggregation connectivity index, NorAI(i, t) and NorCONTAG(i, t) are the normalized values ​​of the aggregation and spread scores of grid i in year t, respectively;

[0323] (4) Landscape heterogeneity index (LHI): The landscape heterogeneity index is represented by the Shannon diversity index (SHDI), which reflects the degree of landscape heterogeneity. The higher the diversity index, the richer the landscape types. When the difference in the proportion of landscape area of ​​each type increases, the diversity decreases. The formula is:

[0324]

[0325] Where: SHDI is Shannon diversity, the ratio of the number of patches to the area in the landscape type is P i It indicates that it reflects the number of landscape elements and the specific proportion of each landscape element;

[0326] The evaluation results of the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region were superimposed with equal weights and divided into five levels according to the natural breakpoint method. The threshold method was used to screen out the areas with an area greater than 200km in the first two levels. 2 The ecological patches are superimposed on the ecological protection areas to obtain the ecological source areas of the ecological areas;

[0327] Step S3 selects seven resistance factor indicators, namely, land use type, altitude, NDVI, night light, distance to road, distance to railway and distance to water, divides the resistance factor level into five categories, and sets the resistance coefficient to 1-5 (Table 4); based on the weight value, the weights are superimposed to obtain the comprehensive ecological resistance surface of the ecological region;

[0328] (1) Land use type

[0329] Land use type affects regional information exchange, material and energy flow; the closer the land use is to the ecological source, the stronger the ecological function is and the smaller the ecological resistance is; on the contrary, the more human interference the land use type is subject to, the greater the ecological resistance is; forest land is the land use type with the least ecological resistance, followed by grassland, water area, cultivated land, unused land and construction land, with resistance values ​​of 1, 2, 3, 4 and 5 respectively;

[0330] (2) Altitude

[0331] Altitude is used to reflect the topography of the ecological area, that is, the higher the altitude, the higher the resistance value; the altitude is divided into five levels: less than 255m, 255m-491m, 491m-748m, 748m-1079m, and greater than 1079m, with resistances of 1, 2, 3, 4, and 5 respectively;

[0332] (3) NDVI

[0333] NDVI is the "normalized difference vegetation index", which is a common index used to detect vegetation growth, vegetation coverage, and dynamic changes in vegetation in time series. NDVI is proportional to vegetation coverage. The larger the NDVI value, the higher the vegetation coverage, the higher the habitat quality, and the smaller the ecological resistance. The natural breakpoint method is used to divide NDVI into five levels: 0-0.296, 0.296-0.384, 0.384-0.469, 0.469-0.563, 0.563-0.804, and the resistance values ​​are 5, 4, 3, 2, and 1 respectively.

[0334] (4) Night lighting

[0335] Night light data is used to reflect the intensity of human activities in the ecological area. That is, the greater the night light data value, the greater the intensity of human interaction, the worse the ecological environment, and the greater the resistance value. The night light is divided into five levels: 0, 0-15, 15-29, 29-47, 47-63, and their resistances are 1, 2, 3, 4, and 5 respectively.

[0336] (5) Distance from the road

[0337] The development of transportation reflects the intensity of human activities. The more developed the transportation is, the more serious the interference to the ecosystem is, the greater the ecological flow resistance is, and the migration and information exchange of organisms are seriously hindered. National roads, provincial roads and expressways are selected as the factors affecting the ecological source, and the distance from the road is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m, and the resistance values ​​of more than 4000m are 5, 4, 3, 2 and 1 respectively.

[0338] (6) Distance from railway

[0339] The distance from the railway is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m, >4000m, and the resistance values ​​are 5, 4, 3, 2, and 1 respectively;

[0340] (7) Distance from water bodies

[0341] Water areas are important pathways for the migration of material energy and species. Therefore, the closer to the water area, the better the ecological quality and the smaller the ecological resistance value. The distance to the water area is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m. If the distance is greater than 4000m, the resistance is 1, 2, 3, 4, 5 respectively.

[0342] Table 4 Determination of single factor resistance coefficient and weight

[0343]

[0344]

[0345] Step S4 is based on each ecological source and the comprehensive ecological resistance surface, based on the circuit theory principle, and determines the corridors, pinch points, and obstacle points between each ecological source through Formula 26, and constructs an ecological security pattern:

[0346]

[0347] Among them, I is the size of the ecological flow passing through the corridor; V is the size of the ecological source measured in the entire corridor; R is the cumulative resistance of the corridor;

[0348] Step S5 quantitatively evaluates source connectivity and network connectivity based on ecological sources and corridors in the ecological region;

[0349] (1) Evaluation of source connectivity Based on the ecological source, the connectivity of the ecological source was analyzed using formulas 27 and 28 to calculate its connectivity probability (PC) and comprehensive connectivity index (IIC):

[0350]

[0351] Among them, PC is the connectivity probability, IIC is the comprehensive connectivity index, n is the total number of ecological sources in the landscape, and a i and a j is the property of ecological source i and ecological source, nl ij is the number of ecological corridors in the shortest path (topological distance) between patches i and j, is the maximum value of the path probability product between ecological source i and ecological source j, A Lis the total area of ​​the ecoregion in question;

[0352] (2) Network connectivity analysis is based on ecological sources and ecological corridors. With reference to landscape ecology, it is simplified into a network structure index of abstract “point-line” graphs. The network closure (α index), network point-line ratio (β index), and network connection ratio (γ index) are selected as indicators for evaluating the connectivity of the ecological network of the ecological region through formulas 29, 30, and 31:

[0353] Network closure (α index): measures the closure of the network, that is, the ratio of the actual number of loops in the network to the maximum possible number of loops in the network, and is the degree of connecting the landscape node loops in the network:

[0354]

[0355] Where: α represents the network closure, L represents the number of corridors, V is the number of nodes, L-V+1 is the actual number of loops, and 2V-5 is the maximum possible number of loops;

[0356] Network point-line ratio (β index): measures the average number of links for each node in the network and is an indicator of the ease with which a landscape node can connect with other landscape nodes:

[0357]

[0358] In the formula: β represents the network point-line rate, L represents the number of corridors, and V is the number of nodes;

[0359] Network connectivity rate (γ index): measures the connectivity of the network and describes the degree to which all landscape nodes in a network are connected:

[0360]

[0361] In the formula: γ represents the network connectivity rate, L represents the number of corridors, and V represents the number of nodes;

[0362] Step S6: This experiment designs a factor control analysis experiment on climate and land use during the observation period to evaluate the contribution of climate change and human activities to ESN changes. This experiment considers scenarios with the ecological security pattern in the 1990s as the baseline, climate change only as the climate scenario, land use change as the human activity scenario, and the ecological security pattern in the 2020s as the climate change and human activity scenario. Attribution analysis is used to quantify the relative contributions of climate change and human activities to the changes in the regional ecological security pattern.

[0363] Table 5 Experimental design of attribution analysis

[0364]

[0365] This experiment estimates the independent impact of climate change or human activities on the ecological security pattern by calculating the difference between the baseline scenario and the climate scenario or human activity scenario, and calculates the interactive effect of climate change and human activities by subtracting the independent impact of climate change and human activities from the combined effect of climate change and human activities. This experiment quantifies the contribution of the combined effect, climate change effect, human activity effect and interactive effect to the connectivity of the ecological security pattern of the described eco-region, thereby achieving attribution analysis of the changes in the ecological security pattern.

[0366] Calculate the relative contribution η a ,η b , η c , η d The expressions are:

[0367]

[0368]

[0369] Among them, η a ,η b , η c , η d Respectively represent the relative contribution rates of interaction effects, human activities, climate change and combined effects, Δr a ,Δr b , Δr c , Δr d They represent the change values ​​of ecological security pattern connectivity caused by interaction effects, human activities, climate change and comprehensive effects respectively, and R is the ecological security pattern connectivity value in the baseline scenario.

Claims

1. The assessment and attribution method of regional ecological security pattern changes is characterized by The evaluation and attribution methods for the changes in the regional ecological security pattern are as follows: S1. Obtain various geospatial data of the ecological region; S2. According to various types of geographic spatial data of the ecological region, obtain the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region, and determine the ecological source area according to the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region; S3. According to the natural conditions of the ecological region and the resistance factors of social and economic activities, the weights are superimposed to obtain the comprehensive ecological resistance surface; S4. According to each ecological source and modified resistance surface, determine the ecological corridors, pinch points, and obstacle points between each ecological source, and build an ecological security pattern; S5. Quantitatively assess source connectivity and network connectivity based on ecological sources and ecological corridors; S6. Conduct attribution analysis on the relative contributions of climate change and human activities to changes in regional ecological security patterns in the ecological regions; The various types of geographic spatial data in step S1 include climate, soil, hydrology, land use and remote sensing image data of the same temporal and spatial scales; Step S2 evaluates the three ecosystem service functions of habitat quality, water conservation, and soil conservation based on climate and land use data, and superimposes them with equal weights to obtain an ecological importance index; (1) Habitat quality: Based on Formula 1, the habitat quality is calculated by combining the suitability of the habitat, the intensity of external threat sources, and the sensitivity. The calculation formula is as follows: Where: Q xj represents the habitat quality of pixel x on land cover type j; H j is the habitat suitability of land cover type j; (2) Water conservation: Water conservation is obtained by subtracting the total runoff of the previous year from the water production. The water production is estimated based on Formula 2; the total annual runoff is estimated by assigning an average runoff coefficient value to each land use type and multiplying it by the annual precipitation of each pixel in that year; Where: Y xj is the annual water yield m of grid cell x on land use type j 3 ;AET xj is the annual average evapotranspiration of grid cell x on land use type j; P x is the annual precipitation of grid cell x; AET xj / P x is the Zhang coefficient, which is an algorithm simulated based on the Budyko equation combined with the actual situation; Runoff ij =P ij ×C j (3) WR ij =Y xj -Runoff ij (4) In the formula, WR ij It is LUCC j The water conservation capacity of each pixel in the year, Y xj It is LUCC j The water output of each pixel in the year, in mm, Runoff ij It is LUCC j The total annual runoff of each pixel surface, in mm / a, P ij LUCC j The annual precipitation of each pixel in the year, in mm, C j represents the jth type of land use; (3) Soil retention: Equation 5 is used to estimate the erosion reduction and sediment retention. Finally, the soil retention is calculated by subtracting the sediment retention from the soil loss. CEDAR x =PKLS x -USLE x (5) In the formula: SEDRET x Represents the soil conservation amount of grid x on a certain land use type; PKLS x Represents the potential soil loss of grid x for a certain land use type, USLE x represents the actual erosion amount of grid x on a certain land use type, that is, the soil erosion amount under vegetation cover and soil and water conservation measures; The three types of ecosystem service function indicators have the same contribution to ecological importance, so habitat quality, water conservation and soil conservation are normalized. The cumulative ecological importance index score of each grid is the sum of the three types of ecosystem service function scores. The calculation formula is as follows: ESI=NorHQ(i,t)+NorWR(i,t)+NorSR(i,t) (6) Where: ESI is the ecological importance index score, NorHQ(i,t), NorWR(i,t), and NorSR(i,t) are the normalized values ​​of habitat quality, water conservation, and soil conservation of grid i in year t, respectively; Step S2 uses climate, land use and remote sensing image data to calculate three sensitivity indexes: sensitivity to nitrogen and phosphorus pollution, sensitivity to cultivated land soil erosion, and sensitivity to human disturbance, and superimposes them with equal weights to obtain the ecological sensitivity index: (1) Nitrogen and phosphorus pollution sensitivity: The nitrogen and phosphorus output of each land use type was estimated based on equations 7, 8, and 9; ALV x =HSS x ×POL x (7) l x =log (∑ U Y U ) (9) Where ALV x is the load value adjusted at grid cell x; POL x is the output coefficient of grid cell x; HSS x is the hydrological sensitivity score at grid cell x, is the average runoff index of the catchment area; λ x is the runoff index at grid cell x, which is used to express the sum of the water yields of all grids in the catchment area; (2) Sensitivity of cultivated land soil erosion: The sensitivity of cultivated land soil erosion is characterized by grain supply, and the net primary productivity (NPP) is combined to achieve a reasonable distribution of grain output on cultivated land in each county of the ecological region. The formula is: Where: Px is the grain yield in the xth pixel, in kg; NPP x is the NPP value of the x-th pixel on cultivated land; P sum is the total grain output of a county; NPP sum It is the sum of NPP on agricultural land in a county; (3) Human interference sensitivity: Based on geographic spatial data such as population density, land use, and road network, population concentration pressure, land use intensity, traffic network pressure, and power consumption pressure are estimated; and the cumulative human interference sensitivity index is obtained by adding up the equal weights; Population concentration pressure calculation: The population density grid data with a resolution of 1km×1km from the LandScan global population database is used; for grids with more than 1,000 people per square kilometer, the interference intensity score is set to 10 points, and for grids with a population density of 1,000 people / km 2 The following grid uses the logarithmic method to calibrate the population density pressure score to obtain the ecological region population concentration pressure index: Among them, PDI s is the population concentration pressure, pd is the population density, PDI s The value range of is [0, 10], the higher the score, the higher the interference intensity; Land use intensity calculation: by assigning values ​​to the relative pressures of various land use forms on natural ecosystems, the land use intensity index of the ecological region is obtained; The land use intensity index is constructed by calculating the area ratio of different land use types in the grid to describe the relative size of regional land use risk. The formula is as follows: Among them, LDI s is the land use intensity; n is the total number of land use types contained in the grid; P i is the percentage of land use type i in the grid area; S i Quantify pressure scores for land-use types; Traffic network pressure measurement: Based on the Euclidean distance, the pressure of the traffic network with different buffer radii is assigned to obtain the traffic network pressure index of the ecological region; Electricity consumption pressure calculation: convert nighttime light data into electricity consumption to characterize the intensity of electricity consumption; fit the relationship between the ecological area's electricity consumption data and the nighttime light intensity value by constructing a power function fitting model, thereby obtaining the ecological area's electricity consumption pressure; And ij =DN sij a (13) In the formula, E ij is the total electricity consumption of city j in year i; DN sij is the sum of the DN values ​​of the nighttime light data of city j in year i; a is the estimated parameter of the power function fitting model; Cumulative human disturbance intensity calculation: The four types of human disturbance indicators contribute equally to the cumulative human pressure, so the population concentration pressure, land use pressure, traffic network pressure and power consumption pressure are normalized. The cumulative human disturbance pressure score of each grid is the sum of the four types of human pressure scores. The formula is: HDI(i,t)=NorPDI s (i,t)+NorLDI s (i,t)+NorTDI s (i,t)+NorEDI s (i,t) (14) Where HDI(i,t) is the cumulative human disturbance pressure score of grid i in year t; NorPDI s (i,t), NorLDI s (i,t), NorTDI s (i,t), NorEDI s (i, t) are the normalized values ​​of population concentration pressure, land use pressure, traffic network pressure and electricity consumption pressure scores of grid i in year t; The three types of ecological sensitivity indicators contribute equally to ecological sensitivity. Therefore, the sensitivity to nitrogen and phosphorus pollution, sensitivity to cultivated land soil erosion, and sensitivity to human disturbance are normalized. The cumulative ecological sensitivity index score of each grid is the sum of the three types of ecological sensitivity scores. The formula is: ESS=NorNDR(i,t)+NorFP(i,t)+NorHDI(i,t) (15) Where: ESS is the cumulative ecological sensitivity index, NorNDR(i,t), NorFP(i,t), and NorHDI(i,t) are the normalized values ​​of the nitrogen and phosphorus pollution sensitivity, cultivated land soil erosion sensitivity, and human disturbance sensitivity scores of grid i in year t, respectively; The landscape dominance index in step S2 can highly condense the landscape pattern information and reflect the characteristics of its structural composition and spatial configuration. The landscape pattern index method can better measure the regional landscape pattern and reflect the changes in the characteristics of the landscape pattern in a specific region from different scales. The landscape fragmentation index, landscape complexity index, landscape aggregation connectivity index, and landscape heterogeneity index are estimated using land use data, and the landscape dominance index is obtained by superimposing them with equal weights. (1) Landscape fragmentation index: The landscape fragmentation index is obtained by the equal-weighted superposition of patch density and average patch area. Patch density reflects the complexity of the landscape spatial structure; the average patch area represents an average condition and can reflect the degree of landscape fragmentation in landscape structure analysis. The formula is: Where: is the PD patch density, N1 is the number of patches contained in patch type i in the landscape, which is used to express the total number of patches in the landscape; A1 represents the total area of ​​the landscape in each ecological region; Where: AREA MN is the average patch area, A2 is the total area of ​​a patch type, and N2 is the total number of patches of this patch type; LFI=NorPD(i,t)+NorAREA_MN(i,t) (18) Where: LFI is the landscape fragmentation index, NorPD(i,t) and NorAREA_MN(i,t) are the normalized values ​​of the patch density and average patch area index score of grid i in year t, respectively; (2) Landscape complexity index: The landscape complexity index is obtained by the weighted superposition of edge density and average patch shape. Edge density is an important indicator of patch shape, indicating the degree to which the type is segmented. MSI is used to reflect the overall complexity of the patch shape of the landscape type. ED=E / A3 (19) Where: ED is the edge density, E is the total length of the patch boundary in the study area, and A3 represents the total area of ​​the landscape. The larger the ED value, the greater the edge density, the higher the degree of segmentation of the landscape type, and the more dispersed the layout. Where: MSI is used to characterize the overall complexity of the patch shape of landscape types, a ij is the area of ​​landscape patch ij, P ij is the perimeter of landscape patch ij, and n is the number of landscape patches; When MSI ≥ 1, all plaques are square, MSI = 1, and the larger the value, the more irregular the plaque shape; LCI=NorED(i,t)+NorMSI(i,t) (21) Where: LCI is the landscape complexity index, NorED(i,t) and NorMSI(i,t) are the normalized values ​​of the edge density and average patch shape index score of grid i in year t, respectively; (3) Landscape aggregation connectivity index: The landscape aggregation connectivity index is obtained by the weighted superposition of aggregation and spread. The larger the aggregation value, the greater the degree of spatial aggregation of patch types. The higher the spread value, the better the connectivity of a certain dominant patch type in the landscape. On the contrary, it reflects that the landscape is a dense pattern with multiple elements and the degree of landscape fragmentation is high. The formula is: Where: AI is the degree of aggregation, g ii is the number of adjacent pixels of the patch type, maxg ii It is the maximum value of the number of adjacent pixels of each patch type, which is used to reflect the spatial aggregation and connectivity of the patch type; Where: CONTAG is the degree of contagion, P i is the area proportion of patch type i in the landscape, g ik is the number of adjacencies between pixels of patch types i and k based on the double counting method, and m is the number of patch types present in the landscape; LAFI=NorAI(i,t)+NorCONTAG(i,t) (24) Where: LAFI is the landscape aggregation connectivity index, NorAI(i,t) and NorCONTAG(i,t) are the normalized values ​​of the aggregation and spread scores of grid i in year t, respectively; (4) Landscape heterogeneity index: The landscape heterogeneity index is characterized by Shannon diversity, which reflects the degree of landscape heterogeneity. The higher the diversity index, the richer the landscape types. When the difference in the proportion of landscape area of ​​each type increases, the diversity decreases. The formula is: Where: SHDI is Shannon diversity, the ratio of the number of patches to the area in the landscape type is P i It indicates that it reflects the number of landscape elements and the specific proportion of each landscape element; The evaluation results of the ecological importance index, ecological sensitivity index and landscape dominance index of the ecological region were superimposed with equal weights and divided into five levels according to the natural breakpoint method. The threshold method was used to screen out the areas with an area greater than 200km in the first two levels. 2 The ecological patches are superimposed on the ecological protection areas to obtain the ecological source areas of the ecological areas; Step S3 selects seven resistance factor indicators, namely, land use type, altitude, NDVI, night light, distance to road, distance to railway and distance to water, divides the resistance factor level into five categories, sets the resistance coefficient to 1-5, and obtains the comprehensive ecological resistance surface of the ecological region based on the weight superposition of each resistance factor indicator; (1) Land use type Land use type affects regional information exchange, material and energy flow; the closer the land use is to the ecological source, the stronger the ecological function and the smaller the ecological resistance; on the contrary, the more human interference the land use type is subject to, the greater the ecological resistance; forest land is the land use type with the least ecological resistance, followed by grassland, water area, cultivated land, unused land and construction land, with resistance values ​​of 1, 2, 3, 4 and 5 respectively; (2) Altitude Altitude is used to reflect the topography of the ecological area, that is, the higher the altitude, the higher the resistance value; the altitude is divided into five levels: less than 255m, 255m-491m, 491m-748m, 748m-1079m, and greater than 1079m, with resistances of 1, 2, 3, 4, and 5 respectively; (3) NDVI NDVI is the "normalized difference vegetation index", which is a common index used to detect vegetation growth, vegetation coverage, and dynamic changes in vegetation in time series. NDVI is proportional to vegetation coverage. The larger the NDVI value, the higher the vegetation coverage, the higher the habitat quality, and the smaller the ecological resistance. The natural breakpoint method is used to divide NDVI into five levels: 0-0.296, 0.296-0.384, 0.384-0.469, 0.469-0.563, and 0.563-0.804, and the resistance values ​​are 5, 4, 3, 2, and 1 respectively. (4) Night lighting Night light data is used to reflect the intensity of human activities in the ecological area. That is, the greater the night light data value, the greater the intensity of human interaction, the worse the ecological environment, and the greater the resistance value. The night light is divided into five levels: 0, 0-15, 15-29, 29-47, 47-63, and their resistances are 1, 2, 3, 4, and 5 respectively. (5) Distance from the road The development of transportation reflects the intensity of human activities. The more developed the transportation is, the more serious the interference to the ecosystem is, the greater the ecological flow resistance is, and the migration and information exchange of organisms are seriously hindered. National roads, provincial roads and expressways are selected as the factors affecting the ecological source, and the distance from the road is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m, and the resistance values ​​of more than 4000m are 5, 4, 3, 2 and 1 respectively. (6) Distance from railway The distance from the railway is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m, >4000m, and the resistance values ​​are 5, 4, 3, 2, and 1 respectively; (7) Distance from water bodies Water areas are important pathways for the migration of material energy and species. Therefore, the closer to the water area, the better the ecological quality and the smaller the ecological resistance value. The distance to the water area is divided into five levels: 0-1000m, 1000m-2000m, 2000m-3000m, 3000m-4000m. If the distance is greater than 4000m, the resistance is 1, 2, 3, 4, 5 respectively. Step S4 is based on each ecological source and the comprehensive ecological resistance surface, based on the circuit theory principle, and determines the corridors, pinch points, and obstacle points between each ecological source through Formula 26, and constructs an ecological security pattern: Among them, I is the size of the ecological flow passing through the corridor; V is the size of the ecological source measured in the entire corridor; R is the cumulative resistance of the corridor; Step S5 quantitatively evaluates source connectivity and network connectivity based on ecological sources and corridors in the ecological region; (1) Evaluation of source connectivity Based on the ecological source, the connectivity of the ecological source is analyzed through formulas 27 and 28 to calculate its connectivity probability and connectivity comprehensive index: Among them, PC is the connectivity probability, IIC is the comprehensive connectivity index, n is the total number of ecological sources in the landscape, and a i and a j is the property of ecological source i and ecological source j, nl ij is the number of ecological corridors in the shortest path between patches i and j, is the maximum value of the path probability product between ecological source i and ecological source j, A L is the total area of ​​the ecoregion in question; (2) Network connectivity analysis is based on ecological source areas and ecological corridors, which are simplified into network structure indexes of abstract "point-line" graphs. Network closure, network point-line ratio, and network connection ratio are selected as indicators for evaluating the connectivity of the ecological network of the ecological region through formulas 29, 30, and 31: Network Closure: Measures the closure of a network, which is the ratio of the actual number of loops in the network to the maximum possible number of loops that exist in the network. It is the degree to which the loops of landscape nodes in the network are connected: Where: α represents the network closure, L represents the number of corridors, V is the number of nodes, L-V+1 is the actual number of loops, and 2V-5 is the maximum possible number of loops; Network point-line ratio: It measures the average number of links per node in the network and is an indicator of how easy it is for a landscape node to connect with other landscape nodes: In the formula: β represents the network point-line rate, L represents the number of corridors, and V is the number of nodes; Network connectivity: measures the connectivity of a network, describing the degree to which all landscape nodes in a network are connected: In the formula: γ represents the network connectivity rate, L represents the number of corridors, and V represents the number of nodes; Step S6 quantifies the relative contributions of climate change and human activities to changes in regional ecological security patterns through attribution analysis; The independent impact of climate change or human activities on the ecological security pattern was estimated by calculating the difference between the baseline scenario and the climate scenario or human activity scenario, and the interactive effect of climate change and human activities was calculated by subtracting the independent impact of climate change and human activities from the combined effect of climate change and human activities. The contribution of the connectivity of the ecological security pattern of the ecological region was indicated by quantifying the combined effect, climate change effect, human activity effect and interactive effect, thus realizing the attribution analysis of the change of the ecological security pattern. Calculate the relative contribution η a ,η b , η c , η d The expressions are: Among them, η a ,η b , η c , η d Respectively represent the relative contribution rates of interaction effects, human activities, climate change and combined effects, Δr a ,Δr b , Δr c , Δr d They represent the change values ​​of ecological security pattern connectivity caused by interaction effects, human activities, climate change and comprehensive effects respectively, and R is the ecological security pattern connectivity value in the baseline scenario.

Images