[0001] 本发明涉及生态保护修复技术领域，具体为一种基于生态系统服务和生态脆弱性的生态管理分区方法。

[0002] 生态系统服务是指人类直接或间接从生态系统中获得的各种惠益，成为连接自然环境与人类福祉的纽带，因而将生态系统服务评估融入到生态管理中已受到广泛关注，将生态系统服务纳入生态管理决策中，识别重点生态功能区，可以更好地实现生态系统保护与恢复的低成本与高效益目标，生态脆弱性是指易受到来自外部压力和扰动和对外界干扰缺乏适应能力的状态，涉及到暴露性(生态系统暴露于外界干扰的程度)、敏感性(生态系统对环境变化和人类干扰的敏感程度)和适应能力(生态系统为适应外界环境变化而不断发展的能力)三个维度，随着全球气候变化和生态环境问题日益突出，脆弱性评估已经成为可持续发展领域的一个核心议题，在生态保护规划中得到了广泛应用，例如，根据生态脆弱性级别制定生态保护措施，确定生态保护重点区域，在区域尺度上提出土地利用规划与管理策略，但是，仅从生态系统服务或生态脆弱性的角度来制定生态管理策略可能会产生偏差；例如，对于生态系统服务价值低、生态脆弱性高的区域，从生态系统服务的角度这些区域不应该受到保护，但在考虑生态脆弱性时这些区域应该被重点关注；

[0003] 近年来，结合生态系统服务重要性与生态敏感性的国土空间规划(如生态空间划定、生态红线划定、生态功能分区、生态安全格局构建等)在国内的长江流域、黑河流域、京津冀地区、江西省、昆明市、南宁市、大理白族自治州、云浮市等区域以及省、市级行政区的生态空间规划和生态安全格局构建中得到广泛应用，其中，生态敏感性指标主要考虑水土流失、土地沙化、石漠化、地质灾害等，然而，开展生态系统服务重要性与生态脆弱性评估的国土空间规划则相对较少；敏感性作为生态脆弱性的一个维度，是生态系统的内在属性，其难以体现人类活动对生态系统的干扰强度，当前，人类活动对自然生态系统的影响十分广泛和深入，生态脆弱性包含了暴露性、敏感性和适应能力，更能反映出人类活动影响下生态系统的实际状况，因此，结合生态系统服务重要性和生态脆弱性的国土空间规划能够在一定程度上弥补当前不足。

[0022] 为了便于本领域技术人员对本技术方案更加清楚，下面将结合附图1‑图5，以祁连山地区为案例区，详细阐述本发明的技术方案。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。

[0023] 本发明提供了一种基于生态系统服务和生态脆弱性的生态管理分区方法，主要包括生态系统服务评估与价值核算、生态脆弱性指标体系构建与评估、基于生态系统服务重要性和生态脆弱性的空间分区等(如图1)，主要步骤如下：

[0024] 步骤1、生态系统服务评估与价值核算

[0025] 生态系统服务评估与价值核算中，本技术方案针对生态脆弱区，选择水源涵养、水土保持、防风固沙、固碳和生境质量5种典型生态系统服务进行评估与价值核算(如图2)。未来可根据区域生态环境特点，有针对性地剔除或增加不同类型生态系统服务，评估其生态系统服务重要性。

[0026] 1)水源涵养服务评估与价值核算

[0027] 水源涵养服务采用水量平衡方程评估，该方程根据降水量和实际蒸散发的差异计算每个栅格的产水量，输入参数包括潜在蒸散量、降水、土地利用类型、植物有效含水量、根系最大深度和生物物理参数表。主要计算公式如下：

[0028]

[0029] 式中，Yxj为土地利用类型j下给定栅格x的年水源涵养量(mm)；AETxj为土地利用类型j下给定栅格x的年实际蒸散发量(mm)；Px表示给定栅格x的年降水量(mm)。

[0030]

[0031] 式中，Rxj为土地利用类型j下给定栅格x的Budyko系数，计算公式如下：

[0032]

[0033] 式中，kxj土地利用类型j下给定栅格x的蒸散发系数；ETox是栅格x的相对蒸散发(mm)；wx给定栅格x的有效水系数，表征植被对土壤水分利用的响应，计算公式如下：

[0034]

[0035] 式中，Z为季节或气候格局系数(范围是1‑10)，表示降水量的大小和分布；AWCx为某一土地利用类型下给定栅格有效含水量。

[0036] 水源涵养价值的计算一般涉及3个方面：1)未经任何处理的天然水资源，可根据当地水利部门制定的水资源使用价格定价；2)利用修建相应库容的水库成本计算，单位库容造价即为单位水资源价格；3)经过处理后的水资源，可用自来水、工业用水价格表示。采用不同方法计算的水价存在差异，可根据区域实际情况合理选取。本技术方案水源涵养服务价值采用市场价值法计算，公式如下：

[0037] ESVwy＝Q×P

[0038] 式中：ESVwy为水源涵养价值(元)；Q为单位面积产水量(m3)；P为水价(元/m3)。

[0039] 2)水土保持服务评估及价值核算

[0040] 水土保持服务采用通用土壤流失方程(RUSLE)通过量化潜在和实际土壤侵蚀量对水土保持服务进行评估。该模型所需的数据主要包括日降水量、DEM、土地利用类型、土壤类型和管理因子数据。本技术方案首先计算每月的土壤保持量，并将其相加得到年土壤保持量。主要计算公式如下：

[0041] WAEP＝R·K·LS·(1‑C·P)

[0042] 式中，WAEP为年水土保持量(t·hm2/a)；R为降水侵蚀因子(MJ·mm/(hm2·h))；K为2 2
土壤可蚀性因子(t·hm·h/(hm·MJ·mm))；LS为坡长坡度因子(无量纲)；C和P分别为植被覆盖因子和管理因子(无量纲)。

[0043] R因子计算公式如下：

[0044]

[0045]

[0046] 式中，Ri是第i月的R因子(MJ·mm·hm‑2·h‑1)；k是每月的天数；Pj第i月第j天的有效降水量，即日降水量≥12mm的降水量，否则记为0mm。

[0047] 参数α和β定义如下：

[0048]

[0049] α＝21.586·β‑7.1891

[0050] 式中，pd12为≥12mm日均降水量；py12日均降水量≥12mm的年均降水量。

[0051] K因子计算公式如下：

[0052]

[0053] 式中，SA、SI、CL、SOC分别为土壤沙粒、粉粒、黏粒和有机碳含量。

[0054] LS因子代表了坡长和坡度对土壤流失的影响，计算公式如下：

[0055]

[0056] λ＝Di/cosθi

[0057] α＝n/(1+n)

[0058] n＝(sinθ/0.0896)/(3.0sin0.8θ+0.56)

[0059]

[0060] 式中，L为坡长因子；λ为坡长因子；α为坡长指数；Di为各栅格沿径流方向坡度的水平投影距离；θi是每个栅格的坡度；S是坡度因子；θ为坡度。

[0061] 植被覆盖因子C计算公式如下：

[0062]

[0063] 式中，VC为植被覆盖度，可以由归一化植被指数计算。

[0064] 水土保持服务价值通过三个方面核算，包括土壤养分保持的价值、减少废弃土地的价值和减少水库泥沙的价值，其计算方法分别为市场价值法、机会成本法和影子工程法。下面分别简要介绍水土保持服务价值核算的三个方面：土壤养分保持价值的估算采用市场价格法。首先，基于《面向陆面模拟的中国土壤数据集》提取区域土壤速效氮、速效磷、速效钾和有机质含量的空间分布数据。进一步，采用如下公式计算单位质量土壤的养分价值：

[0065]

[0066] 式中，ESVwa1为土壤养分保持的单位价值(元/t)；Ri为单位质量土壤中速效氮、速效磷、速效钾和有机质的平均质量(无量纲)；ni为磷酸二铵和氯化钾中速效氮、速效磷和速效钾含量，分别为14％、15％和50％；Pi为磷酸二铵、氯化钾和有机质的价格。

[0067] 减少废弃土地的价值

[0068] 减少废弃土地价值的估算采用机会成本法，计算公式如下：

[0069]

[0070] 式中，ESVwa2为减少废弃土地的单位价值(元/t)；P单位面积土地的年均收益(元/2 3
km)；ρ为土壤容重(t/m)；h为土壤层厚度(m)。本技术方案以食物生产服务价值作为各土地利用的年均收益，即机会成本。

[0071] 减少水库泥沙的价值采用影子工程法估算减少水库泥沙淤积的价值，其公式如下：

[0072]

[0073] 式中：ESVwa3为水库每年减少泥沙淤积的价值(元/t)；C为单位库容造工程造价3 3
(元/m)；ρ为土壤容重(t/m)。

[0074] 将以上三者进行空间叠加，乘以土壤保持量，即得到水土保持服务价值(ESVwa)：

[0075] ESVwa＝(ESVwa1+ESVwa2+ESVwa3)×WAEC

[0076] 3)防风固沙服务评估及价值核算

[0077] 防风固沙服务采用修正风蚀方程(RWEQ)评估，该模型考虑了多种因素(如气候、土壤性质、积雪、地形等)，可以相对准确地预测风蚀量。首先在月尺度上计算风蚀量，然后对其进行累加得到年风蚀量。主要计算公式如下：

[0078]

[0079] Qpmax＝109.8(WF×EF×SCF×K′)

[0080] Sp＝150.71(WF×EF×SCF×K′)‑0.3711

[0081] Qamax＝109.8(WF×EF×SCF×K′×COG)

[0082] Sa＝150.71(WF×EF×SCF×K′×COG)‑0.3711

[0083] 式中，WIEP为防风固沙量；WEp和WEA分别代表潜在和实际风蚀量(kg/m2)；Qpmax和Qamax分别为潜在和实际最大输沙能力(kg/m)；Sp和Sa分别为潜在临界长度(m)和实际临界长度(m)；z为最大下风向侵蚀距离(m)；WF为气象因子(kg/m)；EF、SCF、K’和COG分别为土壤可蚀性因子、土壤结皮因子、地表粗糙度因子和植被因子(无量纲)。

[0084] WF因子计算公式如下：

[0085]

[0086] Wf＝u2(u2‑u1)2×Nd

[0087]

[0088]

[0089]

[0090] SD＝1‑P(snow depth＞25.4mm)

[0091] 式中，Wf为风力因子(m3·s‑3)；ρ为空气密度(kg·m‑3)；g为重力加速度(m·s‑2)；SW‑1为土壤湿度因子(无量纲)；SD为雪盖因子(无量纲)；u1为2m高处的临界风速(m·s )，设定‑1 ‑1
为6m·s ；u2为2m高处的实际风速(m·s )；Nd为每个月内风速超过临界风速的天数；EL为海拔(km)；T是绝对温度(K)；ETp为潜在蒸散发(mm)；R为降水量(mm)；Rd是每月降水天数；SR‑2
太阳辐射(cal·cm )；DT为平均气温(C)；P为降雪深度≥25.4mm的概率。

[0092] EF和SCF计算公式如下：

[0093]

[0094]

[0095] 式中，SA、SI、CL、SOM和CC分别代表沙粒、粉粒、黏粒、土壤有机质和碳酸钙含量(％)。本技术方案假定EF和SCF不随时间改变。

[0096] 植被因子(C)计算公式如下：

[0097] C＝e‑0.0483(VC)

[0098]

[0099] 式中，VC为植被覆盖度(％)；NDVI为归一化植被指数；NDVIsoil为土壤的归一化植被指数；NDVIveg为植被的最大归一化植被指数。在本技术方案中，NDVIveg和NDVIsoil分别代表了累积频率为95％和5％的NDVI值。

[0100] K′计算公式如下：

[0101] K′＝cosα

[0102] 式中，α代表坡度，由DEM数据计算得到。

[0103] 防风固沙服务价值通过两方面核算，即土壤养分保持的价值和减少废弃土地的价值，计算方法参照水土保持服务价值核算方法。

[0104] 4)固碳服务评估及价值核算

[0105] 固碳服务是指植被通过光合作用固存碳，而减少大气中二氧化碳含量的生态系统服务。本技术方案以净初级生产力(Net Primary Productivity，NPP)作为评估固碳服务的代理指标。NPP通过Carnegie‑Ames‑Stanford Approach(CASA)模型进行估算，即通过植被吸收的光合有效辐射和植物将太阳辐射转化为有机物的效率估算。目前，常用固碳服务价值核算方法主要有碳税法和造林成本法。本技术方案采用碳税法核算固碳服务价值，以瑞典碳税率计算，即1200元/t。

[0106] 5)生境维持服务评估及价值核算

[0107] 生境维持服务采用InVEST模型中的生境质量模块(Habitat Quality)进行评估，该模块假设生境质量与人类土地利用强度呈负相关关系，因此利用土地利用数据对每种土地利用类型进行生境适宜性评分。该模型还考虑了人类威胁因素及其对生境质量的负面影响。由于每种威胁对栖息地都有其独特的影响，因此该模型需要相关人员根据实际情况来分配各种威胁因子对不同栖息地类型的影响。总的来说，以下四个方面决定了威胁因子的影响：

[0108] ①各威胁因子的相对影响强度；

[0109] ②威胁因子的最大影响范围和影响方式(线性或指数)。威胁因子的最大影响距离和线性或指数距离衰减函数表示如下：

[0110] 如果为线性

[0111] 如果为指数

[0112] 式中，irxy为源自栅格y的威胁因子r对栅格x的影响；dxy为栅格x和y的线性距离；drmax威胁因素r的最大有效距离。

[0113] ③对干扰的保护水平计算如下：

[0114]

[0115] 式中，Dxj为土地利用类型j下栅格x的总威胁因子水平；R为威胁源的总数量；Yr代表威胁因子r的所有栅格；Wr为威胁因素r的权重，用以表示威胁因子r对各土地利用类型的破坏性；ry为源自栅格y的威胁因子r的相对影响；βx为栅格x的可达性；Sjr为土地利用类型j对威胁因子r的敏感性，变化范围为0‑1，1代表敏感性最高。

[0116] ④各生境类型对每种威胁因子的相对敏感性计算如下：

[0117]

[0118] 式中，Qxj为土地利用类型j下栅格x的生境质量；Hj为土地利用类型j的生境适宜性；z和k是比例常数。

[0119] 本技术方案采用二值法(1表示存在威胁，0表示不存在威胁)提取土地利用数据中的威胁源。通过专家知识和文献综述，获得了每种威胁因素的相对影响、每种生境类型的适宜性以及每种生境类型对每种威胁的相对敏感性(表1)，以及各威胁因子的相对影响强度(表2)。常量z默认设置为2.5，常量k设置为栅格分辨率的一半。

[0120] 表1各生境类型对每种威胁的适宜性和敏感性

[0121]

[0122] 表2各威胁因素的相对影响和距离衰减函数

[0123]

[0124] 在ArcGIS中将基于土地利用数据的各生境类型的生境适宜性赋值数据空间化，与生境维持服务评估结果相除，可以得到各威胁因素对生境的影响系数图层。参考谢高地等(2003)的成果，将区域各土地利用类型生物多样性价值应用到生境维持价值计算中，并将其转化为各土地利用类型生物多样性价值(表3)。进而，通过土地利用数据将各生境类型的生境维持价值空间化。最后，将空间化的价值乘以影响系数图层即得到区域生境维持服务价值空间分布。

[0125] 表3各生境类型的生境维持价值

[0126]

[0127] 6)生态系统服务重要性评估

[0128] 在ArcGIS软件中将各生态系统服务价值进行空间叠加，得到区域生态系统服务总价值，将总价值归一化处理后，即得到区域生态系统服务重要性空间分布。

[0129] 步骤2、生态脆弱性指标体系构建与评估

[0130] 1)生态脆弱性指标体系构建

[0131] 生态脆弱性指标体系构建中，结合对当地生态环境状况分析，采用“暴露性‑敏感性‑适应性”框架，构建了生态脆弱区生态脆弱性评价指标体系(表4)。

[0132] 表4生态脆弱性评估指标及权重

[0133]

[0134] 注：+为正向指标，值越高，脆弱性越高；‑为负向指标，值越高，脆弱性越低。

[0135] 2)生态脆弱性评估

[0136] 采用yaahp(10.3版)软件建立层次分析模型，并基于相关领域专家的打分构建判断矩阵，数据通过一致性检验，最终得到各指标权重(表4)。在进行脆弱性评估之前，采用极差标准化法将所有指标标准化到0‑1之间，以消除指标之间的维度差异。基于生态脆弱性指标数据及其权重结果，在ArcGIS中将指标图层乘以相应权重进行叠加计算，即得到区域暴露性、敏感性、适应性和生态脆弱性空间分布(如图3)，主要计算公式如下：

[0137]

[0138] 式中，EV为生态脆弱性指数；n为指标个数；xi和wi分别为各指标的归一化值和权重。脆弱性的取值范围在0‑1之间，值越高则脆弱性程度越高。

[0139] 步骤3、基于生态系统服务重要性和生态脆弱性的生态管理分区

[0140] 1)生态管理分区划定

[0141] 采用ArcGIS中的优化热点分析工具识别区域生态系统服务重要性和生态脆弱性冷热点。显著热点区(置信度>95％)被归类为高生态系统服务或高生态脆弱区；不显著地区被归类为中等生态系统服务或中等生态脆弱区；显著冷点区(置信度>95％)被归类为低生态系统服务或低生态脆弱区。在此基础上，将生态系统服务重要性和生态脆弱性冷热点图叠加，形成低服务‑低脆弱、低服务‑中脆弱、低服务‑高脆弱、中服务‑低脆弱、中服务‑中脆弱、中服务‑高脆弱、高服务‑低脆弱、高服务‑中脆弱和高服务‑高脆弱共9种类型区(如图4)。

[0142] 2)分区生态管理策略制定

[0143] 区域I人类活动密集，主要土地利用类型为耕地，除粮食生产和固碳服务以外，其他生态系统服务处于较低水平，且良好的气候条件使这一区域生态脆弱性较低。因此，该地区应进行综合开发利用。区域II主要由温性荒漠覆盖，几乎没有人类活动，生态系统服务重要性低，干旱气候导致生态系统非常脆弱。因此，该地区需要采取重点保护和自然恢复的策略。区域III处于干旱和湿润的过渡区，由稀疏草原和高寒荒漠覆盖，生态系统服务重要性适中，但生态脆弱性较高。因此，该区域应受到严格保护，并对主要退化区域进行人工恢复。区域IV主要土地利用类型为草地和森林，为当地居民的生活提供重要的生态系统服务，且适宜的气候条件使该地区具有较低的生态脆弱性。因此，该地区需适度开发以支持当地居民的生产和生活，同时对重点功能区和退化区分别采取重点保护和人工恢复。区域V主要是高寒荒漠，提供重要的水源涵养服务，但高寒气候导致该地区生态脆弱性较高，一旦遭到破坏就难以恢复，且该地区还容易受到气候变暖的影响。因此，该地区应实施严格保护措施，且对于退化严重的地区应辅以适当人工恢复(如图5)。

[0144] 上述技术方案仅体现了本发明技术方案的优选技术方案，本技术领域的技术人员对其中某些部分所可能做出的一些变动均体现了本发明的原理，属于本发明的保护范围之内。