区域生态安全格局变化的评估与归因方法

本发明涉及一种生态安全格局变化的评估与归因方法。

背景技术：

1、生态安全格局是为实现生态安全，对空间布局进行规划和构建的一种模式或结构，它强调生态功能的完整性、连通性和均衡性，以确保生态系统服务的连续性和稳定性。一个合理、科学的生态安全格局可以有效保障生态系统的健康和稳定，因此，生态安全格局是实现生态安全目标的一种重要手段和途径。目前，气候变化、不可持续的土地利用变化和快速城市化等过程极大地破坏了自然生态系统的平衡，使得区域生态系统服务功能降低、生态敏感性区域增多、景观格局变差，生态安全格局发生了明显的变化，进而威胁到了生态安全。所以，急需对区域生态安全格局的变化进行评估并量化其变化驱动因子的贡献率。当前对生态安全格局的有效性或合理性缺少评价方法，针对影响生态安全格局变化的驱动力因素研究更是鲜少有突破性的成果，学者们的关注点集中在将已构建的区域生态安全格局与现有的土地利用方式及区域发展需求进行对比，定性地评价当前生态安全格局的优劣，定量评价的研究远远不够，这使得区域生态安全格局的保护与修复工作开展缺乏可靠的理论支撑。综上所述，现有技术存在的问题是：现技术还缺乏准确评估区域生态安全格局变化及归因分析其驱动因素的评价方法。

技术实现思路

1、本发明的目的是为了解决现有技术无法实现对区域生态安全格局变化的定量评估及归因分析的技术问题，提供了一种区域生态安全格局变化的评估与归因方法。

2、1、区域生态安全格局变化的评估与归因方法，其特征在于所述区域生态安全格局变化的评估与归因方法如下：

3、s1.获取所述生态区域的各类地理空间数据；

4、s2.根据生态区域各类地理空间数据，得到所述生态区域生态重要性指数、生态敏感性指数及景观优势度指数，并根据所述三类指数确定生态源地；

5、s3.根据生态区域自然条件和社会经济活动的各个阻力因子，权重叠加得到综合生态阻力面；

6、s4.根据每个生态源地和修正阻力面，确定各个生态源地之间的生态廊道、夹点、障碍点，并构建生态安全格局；

7、s5.根据生态源地及生态廊道，定量评估源地连通性和网络连通性；

8、s6.对所述生态区域，根据气候变化和人类活动对区域生态安全格局变化相对贡献进行归因分析；

9、步骤s1所述各类地理空间数据包含时空尺度相同的气候、土壤、水文、土地利用及遥感影像数据；

10、步骤s2基于气候和土地利用数据，评估生境质量、水源涵养、土壤保持这三种生态系统服务功能，并将其等权叠加得到生态重要性指数；

11、(1)生境质量:基于公式1，结合生境的适宜度、外界威胁源的强度以及敏感度来计算生境质量，计算公式如下：

12、

13、式中：qxj代表j型土地覆盖上x像元的生境质量；hj为j型土地覆被的生境适宜性；

14、(2)水源涵养：水源涵养由产水量减去年径流总量所得，基于公式2估算产水量；通过对各土地利用类型进行平均径流系数值赋值，乘以各个像元当年的年降水量，估算出年径流总量；

15、

16、式中：yxj为土地利用类型j上的栅格单元x的年产水量(m3)；aetxj为土地利用类型j上栅格单元x年平均蒸散发量；px为栅格单元x的年降水量；aetxj/px为zhang系数，是根据budyko方程结合实际情况模拟的算法；

17、runoffij＝pij×cij    (3)

18、wrij＝yxj-runoffij    (4)

19、式中，wrij是luccj上各个像元当年的水源涵养量，yxj是luccj上各个像元当年的产出水量，单位mm，runoff是luccj上各个像元地表的年径流总量，单位mm/a，pij表示luccj上各个像元当年的年降水量，单位mm，cj表示第j类土地使用方式；

20、(3)土壤保持量:通过公式5，估算侵蚀减少量和泥沙截留量，最终通过土壤流失量减去泥沙截留量求得土壤保持量；

21、sedret=pklsx-uslex    (5)

22、式中：sedretx表示某类土地利用类型上栅格x的土壤保持量；pklsx表示某类土地利用类型上栅格x潜在土壤流失量，uslex分别表示某类土地利用类型上栅格x的实际侵蚀量，即有植被覆盖和水土保持措施下的土壤侵蚀量；

23、三类生态系统服务功能指标对生态重要性的贡献度相同，因此对生境质量、水源涵养及土壤保持进行归一化处理，各格网的累计生态重要性指数得分是三类生态系统服务功能得分的总和，计算公式如下:

24、esi＝norhq(i，t)+norwr(i，t)+norsr(i，t)    (6)

25、式中：esi是生态重要性指数得分，norhq(i，t)、norwr(i，t)、norsr(i，t)分别为t年网格i的生境质量、水源涵养、土壤保持量的归一化值；

26、步骤s2利用气候、土地利用和遥感影像数据测算氮磷污染敏感性、耕地土壤侵蚀敏感性、人类干扰敏感性三种敏感性指数，并将其等权叠加得到生态敏感性指数：

27、(1)氮磷污染敏感性：基于公式7、8、9估算出各土地利用类型的氮、磷输出量；

28、alvx=hssx×polx    (7)

29、

30、λx＝log(∑uyu)    (9)

31、式中，alvx为栅格单元x处调节的载荷值；polx为栅格单元x的输出系数；hssx为栅格单元x处的水文敏感性得分值，为集水区平均径流指数；λx为栅格单元x处的径流指数，为集水区所有栅格的产水量之和；

32、(2)耕地土壤侵蚀敏感性：使用粮食供给表征耕地土壤侵蚀的敏感性程度，结合净初级生产力npp实现所述生态区域各县域粮食产量在耕地上的合理分配，公式为：

33、

34、式中:px为第x个像元上的粮食产量，单位kg；npp为第x像元在耕地上的npp值；p为某县区粮食总产量；npp为某县区npp在农业用地上的总和；

35、(3)人类干扰敏感性：基于人口密度、土地利用、路网这些地理空间数据估算人口聚集压力、土地利用强度、交通路网压力、电力消耗压力；并将其等权叠加得到累计人类干扰敏感性指数；

36、人口聚集压力测算：使用landscan全球人口数据库1km×1km分辨率的人口密度栅格数据；对于每平方千米超过1000人的格网，将干扰强度评分定为10分，对人口密度在1000人/km2以下的格网，采用取对数的方法标定人口密度的压力得分，得到生态区域人口聚集压力指数：

37、

38、其中，pdis为人口聚集压力，pd为人口密度，pdis的值域为[0，10]，得分越高，干扰强度越高；

39、土地利用强度测算：通过对各土地利用形式对自然生态系统产生的相对压力大小进行赋值，得到所述生态区域土地利用强度指数；

40、通过计算格网内不同土地利用类型面积比例来构建土地利用强度指数，用以描述区域土地利用风险的相对大小，公式如下：

41、

42、其中，ldis为土地利用强度；n为格网中包含的土地利用类型总数；pi为土地利用类型i占栅格面积的百分比；si为土地利用类型的压力量化分值；

43、交通路网压力测算：基于欧式距离对不同缓冲半径的交通路网进行压力赋值，得到所述生态区域交通路网压力指数；

44、电力消耗压力测算：将夜间灯光数据转换为电力消耗量，用以表征电力消耗强度；通过构建幂函数拟合模型对所述生态区域电力耗能数据和夜间灯光强度值之间的关系进行拟合，从而得到所述生态区域电力消耗压力；

45、eij＝dnsija    (13)

46、式中，eij为第i年j市总电力消耗量；dnsij；为第i年j市的内夜间灯光数据dn值的总和；a为幂函数拟合模型的预估参数；

47、累计人为干扰强度测算：四类人为干扰指标对累计人为压力贡献度相同，因此对人口聚集压力、土地利用压力、交通路网压力及电力消耗压力进行归一化处理，各格网的累计人为干扰压力得分是四类人类压力得分的总和，公式为:

48、hdi(i，t)＝norpdis(i，t)+norldis(i，t)+nortdis(i，t)+noredis(i，t)    (14)

49、其中，hdi(i，t)为t年网格i的累计人为干扰压力得分；norpdis(i，t)、norldis(i，t)、nortdis(i，t)、noredis(i，t)分别为t年网格i的人口聚集压力、土地利用压力、交通路网压力及电力消耗压力得分的归一化值；

50、三类生态敏感性指标对生态敏感性的贡献度相同，因此对氮磷污染敏感性、耕地土壤侵蚀敏感性及人类干扰敏感性进行归一化处理，各格网的累计生态敏感性指数得分是三类生态敏感性得分的总和，公式为:

51、ess＝norndr(i，t)+norfp(i，t)+norhdi(i，t)    (15)

52、式中：ess是累计生态敏感性指数，norndr(i，t)、norfp(i，t)、norhdi(i，t)分别为t年网格i的氮磷污染敏感性、耕地土壤侵蚀敏感性、人类干扰敏感性得分的归一化值；

53、步骤s2所述景观优势度指数能够高度浓缩景观格局信息，反映其结构组成和空间配置方面的特征，景观格局指数法能够较好的度量区域景观格局，并从不同尺度上反映特定区域景观格局特征变化；利用土地利用数据估算景观破碎化指数、景观复杂程度、景观聚集连通性指数、景观异质性指数，并将其等权叠加得到景观优势度指数；

54、(1)景观破碎化指数：景观破碎化指数由斑块密度和平均斑块面积等权叠加得到；斑块密度反映了景观空间结构的复杂性；平均斑块面积代表一种平均状况，在景观结构分析中可以反映景观破碎化程度，公式为：

55、

56、式中：为pd斑块密度，n为景观中斑块类型i所包含的斑块数量；n为景观中的斑块总数，a表示各生态区域的景观的总面积；

57、

58、式中：areamn为平均斑块面积，a为某斑块类型总面积，n为此斑块类型的斑块总数；

59、lfi＝norpd(i，t)+norarea_mn(i，t)    (18)

60、式中：lfi为景观破碎化指数，norpd(i，t)、norarea_mn(i，t)分别为t年网格i的斑块密度和平均斑块面积指数得分的归一化值；

61、(2)景观复杂度指数：景观复杂度指数由边缘密度和平均斑块形状等权叠加得到；边缘密度用来反映斑块形状的重要指标，表示类型被分割的程度；msi用来反映景观类型斑块形状的总体复杂程度；

62、ed＝e/a    (19)

63、式中：ed为边缘密度，e为研究区内斑块边界总长度，a表示景观总面积；ed值越大，则表示边缘密度越大，景观类型被分割程度越高，布局越分散；

64、

65、式中：msi用来表征景观类型斑块形状的总体复杂程度，aij为景观斑块ij的面积，pij为景观斑块ij的周长，n为景观斑块数目；

66、msi≥1，所有斑块均为正方形时，msi＝1，值越大，说明斑块形状越不规则；

67、lci＝nored(i，t)+normsi(i，t)    (21)

68、式中：lci为景观复杂度指数，nored(i，t)、normsi(i，t)分别为t年网格i的边缘密度和平均斑块形状指数得分的归一化值；

69、(3)景观聚集连通性指数：景观聚集连通性指数由聚集度和蔓延度等权叠加得到；聚集度值越大，反映斑块类型在空间上的聚集程度越大；蔓延度值越高反映景观中的某种优势斑块类型形成了良好的连接性；反之则反映景观是具有多种要素的密集格局，景观的破碎化程度较高，公式为：

70、

71、式中：ai为聚集度，gii为斑块类型的两两相邻像素的数量，maxgii为斑块类型两两相邻像素数量的最大值；它反映斑块类型的空间聚集度和连接程度；

72、

73、式中：contag为蔓延度，pi为斑块类型i在景观中的面积比重，gik为基于双计数方法的斑块类型i和k的像素之间的邻接的数量，m为景观中存在的斑块类型的数量；

74、lafi＝norai(i，t)+norcontag(i，t)    (24)

75、式中：lafi为景观聚集连通性指数，norai(i，t)、norcontag(i，t)分别为t年网格i的聚集度和蔓延度得分的归一化值；

76、(4)景观异质性指数：景观异质性指数采用香农多样性表征，反映景观异质性程度，多样性指数最高，景观类型越丰富；当各类型景观面积比例差异增大时，多样性下降，公式为：

77、

78、式中：shdi为香农多样性，景观类型中斑块数与面积之比用pi表示，反映的是景观要素多少与各景观要素所占的具体比例；

79、通过等权叠加所述生态区域生态重要性指数、生态敏感性指数及景观优势度指数的评估结果，按照自然断点法将其划分为五个等级，使用阈值法筛选出前两个等级中面积大于200km2的生态斑块；并叠加生态保护区得到所述生态区域的生态源地；

80、步骤s3选择土地利用类型、海拔、ndvi、夜间灯光、距道路距离、距铁路距离及距水域距离七个阻力因子指标，将阻力因子等级分为五类，阻力系数设置为1-5；基于的权重值，权重叠加得到所述生态区域的综合生态阻力面；

81、(1)土地利用类型

82、土地利用类型影响区域信息交流、物质和能量流动；土地利用方式越接近生态源地，生态功能越强，生态阻力就越小，相反，土地利用类型受到人为干扰越多，生态阻力就越大；将林地作为生态阻力最小的土地利用类型，其次是草地、水域，耕地，未利用地，建筑用地，阻力值分别是1、2、3、4、5；

83、(2)海拔

84、海拔用于反映所述生态区域的地形地貌，即海拔越高阻力值越高；将海拔划分为五个等级：小于255m、255m-491m、491m-748m、748m-1079m、大于1079m，其阻力分别是1、2、3、4、5；

85、(3)ndvi

86、ndvi为“归一化植被指数”，是用来检测植被生长状况、植被覆盖度、时间序列上植被动态变化的常用指数；ndvi与植被覆盖度成正比，ndvi值越大，植被覆盖度就越高，生境质量就越高，生态阻力就越小；采用自然断点法将ndvi分为五个等级：0-0.296、0.296-0.384、0.384-0.469、0.469-0.563、0.563-0.804，其阻力值分别为5、4、3、2、1；

87、(4)夜间灯光

88、夜间灯光数据用于反映所述生态区域的人类活动强度，即夜间灯光数据值越大人类互动强度越大，生态环境越差，阻力值越大；将夜间灯光划分为五个等级：0、0-15、15-29、29-47、47-63，其阻力分别是1、2、3、4、5；

89、(5)距道路距离

90、交通的发展状况反映了人类活动的剧烈程度，交通越发达对生态系统的干扰越严重，生态流动阻力就越大，严重阻碍了生物的迁徙和信息交流；选择国道、省道、高速公路作为道路对生态源地的影响因素，距道路距离划分为五个等级：分别为0-1000m，1000m-2000m，2000m-3000m，3000m-4000m，大于4000m，其阻力值分别为5、4、3、2、1；

91、(6)距铁路距离

92、将距铁路距离划分为5个等级：分别为0-1000m，1000m-2000m，2000m-3000m，3000m-4000m，>4000m，其阻力值分别为5、4、3、2、1；

93、(7)距水域距离

94、水域是物质能量和物种迁移的重要路径，因此离水域越近，生态质量越好，生态阻力值越小；将距水域距离划分为5个等级：分别为0-1000m，1000m-2000m，2000m-3000m，3000m-4000m，大于4000m，其阻力分别是1、2、3、4、5；

95、步骤s4以每个生态源地和综合生态阻力面为基础，基于电路理论原理，通过第26公式确定各个生态源地之间的廊道、夹点、障碍点，并构建生态安全格局：

96、

97、其中，i为穿过廊道的生态流大小；v为整个廊道测量的生态源大小；r为廊道的累积阻力；

98、步骤s5根据所述生态区域的生态源地和廊道定量评估源地连通性和网络连通性；

99、(1)源地连通性的评估以生态源地为基础通过公式27、28对生态源地连通性进行分析，计算其连通性概率、连通性综合指数：

100、

101、

102、其中，pc为连通性概率，iic为连通性综合指数，n是景观中生态源的总数，ai和aj是生态源i和生态源j的属性，nlij是斑块i和j之间最短路径中的生态走廊数量，是生态源i和生态源j之间路径概率乘积的最大值，al是最大横向属性，al是所述生态区域的总面积；

103、(2)网络连通性分析以生态源地和生态廊道为基础，将其简化为抽象“点—线”图形的网络结构指数，选择网络闭合度、网络点线率、网络连接率通过公式29，30，31作为评估所述生态区域生态网络连通性的指标：

104、网络闭合度：测度网络的闭合度，即网络中实际回路数与网络中存在的最大可能回路数之比，是连接网络中景观节点环路的程度：

105、

106、式中：α代表网络闭合度，l代表廊道数，v为节点数，l-v+1为实际环路数，2v-5为最大可能环路数；

107、网络点线率：测度网络中每个节点的平均连线数，是度量一个景观节点与其他景观节点联系难易程度的指标：

108、

109、式中：β代表网络点线率，l代表廊道数，v为节点数；

110、网络连接率：测度网络的连接度，描述一个网络中所有景观节点被连接的程度：

111、

112、式中：γ代表网络连接率，l代表廊道数，v为节点数；

113、步骤s6通过归因分析来量化气候变化和人类活动对区域生态安全格局变化相对贡献；

114、通过计算基线情景与气候情景或人类活动情景之间的差异来估算气候变化或人类活动对生态安全格局的独立影响，并通过从气候变化和人类活动的综合效应中减去气候变化和人类活动的独立影响效应，计算了气候变化和人类活动的交互效应；通过量化了综合效应、气候变化效应、人类活动效应及交互效应其所述生态区域生态安全格局的连通性的贡献作用，从而实现对生态安全格局的变化进行归因分析；

115、计算相对贡献ηa,ηb，ηc，ηd的表达式分别为：

116、

117、

118、

119、

120、其中，ηa,ηb，ηc，ηd分别表示交互效应、人类活动、气候变化及综合效应的相对贡献率，δra,δrb，δrc，δrd分别表示交互效应、人类活动、气候变化及综合效应导致的生态安全格局连通性变化值，r为基线情景中的生态安全格局连通性值。

121、本发明根据源地连通性和网络连通性对区域生态安全格局的稳定性进行量化，并基于气候变化和人类活动的不同情景来量化气候变化和人类活动对区域生态安全格局变化的影响。实现对区域生态安全格局变化的定量评估及归因分析。