自然保护地生态功能空间识别是优化国土空间开发保护格局、开展自然资源针对性保护、落实差异化生态修复策略和采取精细化分区管控的基础性工作^[1—4]。在国家公园情景中, 生态功能分区则是国家公园综合规划和科学管理的关键环节^[5], 是协调自然资源严格保护和合理利用的前提性工作, 对于打破生态环境保护与保护区周边社区民生福祉之间的对立冲突具有现实意义^[6]。我国10个国家公园体制试点中, 有2个国家公园(三江源国家公园、祁连山国家公园)位于青藏高原, 目前体制试点目标任务基本完成。2021年10月12日在《生物多样性公约》第十五次缔约方大会上, 我国正式宣布设立三江源、大熊猫、东北虎豹、海南热带雨林、武夷山等第一批国家公园^[7], 标志着建立以国家公园为主体的自然保护地体系迈出坚实一步。从发展沿革看, 我国自然保护地体系建设发展大致经历了纯自然保护阶段、自然保护与区域经济协调发展阶段, 目前正朝着向可持续发展阶段优化升级。在战略规划层面, 推动自然保护地体系顶层设计, 优化自然保护地空间格局, 加快生态安全战略格局关键区保护修复等逐渐达成共识^[8—9]。在技术操作层面, 现有自然保护地与生态功能区还存在空间分布不匹配问题^[10], 开展大范围、深层次的生态系统服务实证研究是生态功能关键区识别和科学管理生态资产的重要依据^[11]。因此, 有必要对国家公园生态资源功能价值、生态环境演化风险阈值进行系统、全面评估, 促进生态空间的科学、完整识别^[12], 进而建立一套精准化、可视化、靶向性的生态空间分区建设方案, 以此支撑国家公园建设兼具操作性且符合科学规律。可见, 开展国家公园生态空间分区识别具有重要的理论实践意义, 理应成为地理学、生态学、管理学、经济学等多学科维度、多元视角关注的时代课题。

国家公园承担着维护国家生态安全和提供公共生态产品的艰巨任务, 提供了各项生态系统服务, 也面临着生态系统服务功能退化甚至丧失的潜在风险, 过度汲取可能导致生态系统服务价值盈余潜力下降^[13], 生态价值溢出效应减退^[14], 生态环境演化逐渐逼近风险预警阈值, 进而影响到人类福祉的有序提升。可以看出, 在时空上生态系统服务价值和景观生态风险之间存在相互作用^[15], 是生态环境质量正反两个方面的体现, 也是生态安全评估的两个维度, 为国家公园生态空间分区识别提供了新思路。现阶段, 生态系统服务价值和景观生态风险评价研究正由两个独立的科学命题逐步走向二者融合研究, 并且出现了两种截然不同的研究范式和思考维度。一方面, 以人为本的生态文明理念环境议题开始回溯人类福祉本源, 尝试将生态系统服务纳入生态风险评价框架^[16—19]开展研究工作并逐渐发展和完善, 成为生态安全研究的热点话题和重要分支, 目前定量实证案例较少；另一方面, 伴随地理学研究生态化、生态学研究宏观化的影响, 基于生态系统服务价值和景观生态风险时空变化及关联特征研究不断创新^[20—23], 以谋划生态空间分区治理和制定自然资源管理方案为目标, 从生态系统服务价值和景观生态风险正反两个方面来表征和评估区域生态安全状况^{[15, 24]}, 开创了生态空间功能分区识别研究新视角, 相关研究处于不断探索阶段, 尚为形成统一的标准范式。

祁连山是我国西北干旱半干旱地区重要战略资源储备基地, 肩负着为河西走廊、柴达木盆地、青海湖流域、河湟谷地等提供生态产品的艰巨重任, 生态系统服务功能外溢效应显著。在生态学意义上祁连山是一个具有明显高原地带性和地理地缘特征的绿色半岛, 在气象学意义上祁连山是我国东部湿润地区伸入西北荒漠中的一个湿岛。2017年中办、国办印发《祁连山国家公园体制试点方案》, 祁连山国家公园体制试点工作全面启动。至此, 如何实现祁连山国家公园生态系统服务稳定供给？如何界定祁连山国家公园生态环境质量风险预警阈值？如何科学制定祁连山国家公园管控分区方案？成为祁连山国家公园建设绕不开的现实命题。目前, 祁连山国家公园国土空间规划与用途管制中还存在“顶层设计”和“基层匹配”的“剥离现象”亟待解决。将生态系统服务价值和景观生态风险(价值—风险)综合纳入同一研究范畴, 可实现土地资源空间信息表达的定量化、生态系统服务价值的可视化、生态风险预警预估的精细化, 有效揭示国土空间生态安全的关联图式及空间规律。

祁连山国家公园依托国家林业和草原局驻西安森林资源监督专员办事处管理, 加挂祁连山国家公园管理局的牌子^[25]。在体制试点过程中, 为更好地服务国家公园管理目标^[26], 保证管理机构行使管理权的有效性和可行性, 甘肃省、青海省分别组建了大熊猫祁连山国家公园甘肃省管理局和祁连山国家公园青海省管理局两个实体机构, 初步探索形成了统一事权、分级行使所有权的协同管理机制。在现有行政区组织管理模式下, 祁连山国家公园管理形成了甘肃片区、青海片区两个稳定的功能组织及内部联系, 地方事权职责分工清晰。依据因“区”制宜, 提高分区治理针对性的原则^[27], 开展祁连山国家公园生态分区识别诊断时需考虑行政区划的独立性, 以便提高分区管控方案的行政执行力。基于此, 本研究以祁连山国家公园青海片区(以下简称为“园区”)为例, 集成遥感技术、地理信息模型方法、景观生态学方法、GIS网格法构建价值(风险)评价模型, 分析土地利用、生态系统服务价值及景观生态风险的时空演变特征, 提出园区生态空间分区差异化管控方案, 以期为科学制定国家公园生态资产管理制度、精准实施国土空间生态修复工程项目和保障青藏高原生态安全等提供科技支撑和研究思路。

1 数据与方法 1.1 研究区概况

祁连山国家公园位于青藏高原东北部, 横跨甘肃和青海两省, 总面积为5.02万km², 其中青海片区面积1.58万km², 包括核心保护区(0.94万km²)和一般控制区(0.64万km²), 涉及德令哈市、祁连县、天峻县和门源县4县(市), 包括青海省祁连山省级自然保护区、仙米国家森林公园、祁连黑河源国家湿地公园等。祁连山国家公园地处中纬度北温带, 深居内陆, 远离海洋, 属大陆性高寒半湿润山地气候, 冬季长而寒冷干燥, 夏季短而温凉湿润, 园区内生态系统独特, 自然景观多样。祁连山国家公园蕴含丰富的草地、林地、湿地、冰川等生态资源, 是极其重要的水源涵养生态功能区, 对维护青藏高原生态平衡、维持河西走廊绿洲稳定、保障北方地区生态安全具有重要意义, 是国家重要的生态安全屏障(图 1)。

1.2 数据来源与处理

本研究中园区1998年、2008年、2018年土地利用数据以landsat-5/TM、landsat-8/OLI遥感影像数据为数据源。根据研究区土地资源特征和研究目标, 结合野外实际考察, 参照土地利用现状分类标准(GB/T 21010—2017)^[28]和土地利用/土地覆盖遥感监测数据分类体系(http://www.resdc.cn/), 把研究区的土地利用/覆被类型划分为旱地、有林地、灌木林、高覆盖度草地、中覆盖度草地、低覆盖度草地、人工牧草地、水域、永久性冰川雪地、建设用地、沼泽地、其他用地等共12种类型。从生态系统类型来看, 旱地属于农田生态系统, 有林地、灌木林属于森林生态系统, 高覆盖度草地、中覆盖度草地、低覆盖度草地、人工牧草地属于草地生态系统, 水域、永久性冰川雪地、沼泽地属于水体与湿地生态系统, 其他用地是以裸岩石质地、裸土地、戈壁等为主的荒漠生态系统。遥感影像土地利用/覆被现状解译工作, 采用计算机智能与人工目视解译互补的方式进行。通过建立遥感影像解译标志库、数据内业编辑处理、外业调查与核查等工作, 最终形成土地利用空间信息数据库, 经Kappa系数检验和野外实地考察验证, 精度达87.60%以上。研究区土地利用类型数据均采用CGCS2000坐标系, 投影为Albers Conical Equal Area, 中央经线为110°E, 双标准纬线为25°N和47°N, 投影起始纬度12°, 中央经线偏差和起始点偏差都为0。

1.3 研究方法 1.3.1 评价单元划分与数据重构

网格尺度强大的空间信息挖掘力和空间图谱表达力是有效揭示地理事物或现象关联图式及空间规律的重要手段和路径方法^[13]。本研究结合园区地势起伏、面积大小、环境特征及后续生态功能分区识别等因素, 综合考虑国土空间规划、生态管控措施、生态补偿实施的易操作性及实践指导价值, 最终确定1 km×1 km网格为本研究的评价单元, 共计17394个。借助ArcGIS 10.6软件Greate Fishnet、Dissolve、Clip等工具进行格网重采样和计算, 最终实现格网尺度下土地利用类型空间数据的微观重构。

1.3.2 生态系统服务价值估算

本研究选取基于货币量的价值当量法^[29]评价生态系统服务价值(Ecosystem Services Value, ESV), 以谢高地等^[30]提出的生态系统服务价值当量为基准, 结合研究区实际情况构修正得到园区单位面积ESV当量。具体来说, 旱地、有林地、灌木林、水域、永久性冰川雪地、沼泽地、其他用地分别选用旱地、针叶、灌木、水系、冰川积雪、湿地、荒漠的生态系统服务价值当量作为基准当量值；高覆盖度草地、中覆盖度草地、低覆盖度草地则参照草地(草原、草甸、灌草丛)基准当量, 按照比例权重适当修正, 且进一步修正人工牧草地的基准当量；建设用地则不考虑生态系统服务价值^[14]。以网格化土地利用数据为基础, 计算网格单元内各土地利用类型的ESV。

每个网格单一土地利用生态系统服务价值公式表达为：

每个网格生态系统服务价值公式表达为：

式中, ESV_j为单元网格中第j类土地利用类型的生态系统服务价值(10⁸元)；A_ij为第j个网格第i种土地利用类型的分布面积(hm²)；VC_i为第i种土地利用类型的生态系统服务价值系数；n为土地利用类型。

单位当量的经济价值量等于当年全国平均粮食价格的1/7^{[14, 31]}, 本研究粮食产量、粮食价格数据来源于青海省国民经济和社会发展统计公报、《青海省统计年鉴》和国家粮食和物资储备局(lswz.gov.cn), 最终确定标准当量基准值的经济价值为1262.66元/hm²。1998—2018年园区单位面积生态服务价值如表 1所示。

1.3.3 景观生态风险指数构建

参考前人研究成果^[32—34], 结合土地利用景观格局和区域生态风险之间的经验联系, 本研究引入景观干扰度指数、景观破碎度指数、景观分离度指数、景观优势度指数、景观脆弱度指数构建了研究区景观生态风险指数(Ecological Risk Index, ERI), 定量评价风险小区的生态风险大小。通过Fragatats 4.2软件计算景观指数, 最终得到每个风险单元ERI, 以此作为每个风险小区中心点的生态风险值。计算公式为：

式中, n为景观类型数量, A_ki为风险小区k中第i类景观类型面积, A_k为风险小区k的总面积, E_i为景观干扰度指数；V_i为景观脆弱度指数；C_i为景观破碎度指数；S_i为景观分离度指数；DO_i为景观优势度指数；a、b、c分别表示相应景观指数的权重, a+b+c=1, a、b、c权重分别是0.5、0.3和0.2。

结合研究区实际对12种土地利用类型脆弱性分级, 进行归一化处理得到各景观类型的脆弱度指数。永久性冰川雪地、其他用地、人工牧草地、沼泽地、水域、旱地、低覆盖度草地、中覆盖度草地、高覆盖度草地、灌木林、有林地、建设用地的景观脆弱度指数分别为0.6565、0.6209、0.5838、0.5449、0.5039、0.4601、0.4129、0.3613、0.3035、0.2369、0.1562、0.0500。

1.3.4 价值—风险生态分区象限划分

基于生态系统服务价值、景观生态风险指数两个维度综合测度, 对祁连山国家公园青海片区17394个网格单元进行Z-score标准化, 据此划分出4个象限, x轴表示标准化后的生态系统服务价值, y轴表示标准化后的景观生态风险值, 第一象限(I)表示高生态系统服务价值—高生态风险, 第二象限(II)低生态系统服务价值—高生态风险, 第三象限(III)表示低生态系统服务价值—低生态风险, 第四象限(IV)表示高生态系统服务价值—低生态风险, 落入不同象限的网格单元隶属于不同生态分区^[35]。价值—风险生态分区象限划分Z-score标准化计算公式如下：

式中, x为网格单元标准化后的生态系统服务价值和景观生态风险值, x_i为第i个网格单元的生态系统服务价值和景观生态风险值, x为研究区平均值, S为研究区标准差, n为研究区网格单元总数。

2 结果与分析 2.1 土地利用类型变化特征

祁连山国家公园青海片区总面积为15830.67 km², 土地利用包括林地、草地、水域、冰川积雪、裸岩石质地等类型, 呈现“山水林田草冰沙”镶嵌交错的分布格局, 具有明显的高原地带性和生态地理地缘特征。具体来看(图 2), 受地理位置、海拔高度、气候条件的影响, 草地是最主要的土地覆被类型, 构成了高寒植被生态系统的主体, 占园区面积的55.00%以上；2018年高覆盖度草地、中覆盖度草地、低覆盖度草地分别为4098.71 km²、2466.51 km²和2181.47 km², 分别占园区面积的25.89%、15.58%和13.78%；林地资源为1347.79 km², 占园区面积的8.51%；水体与湿地生态系统也是园区重要的生态系统类型, 包括水域、永久性冰川雪地和沼泽地, 2018年面积分别为461.39 km²、527.53 km²和615.01 km², 占园区面积的10.13%；其他用地主要包括裸岩石质地、裸土地和戈壁, 为4063.48 km², 占园区面积的25.67%。

通过土地利用变化转移矩阵进一步分析各土地利用类型的相互转化情况(图 3), 30年间(1998—2018年)转移总面积为102.49 km²。从转出来看, 中覆盖度草地转出面积最大(34.83 km²), 转为人工牧草地、高覆盖度草地的面积分别为17.52 km²和15.31 km²；其次是低覆盖度草地, 转为中覆盖度草地、人工牧草地的面积分别为13.14 km²、7.09 km²；沼泽地、其他用地、永久性冰川雪地、高覆盖度草地的转出面积也较大, 分别为11.88 km²、10.42 km²、9.27 km²和8.31 km²；高覆盖度草地主要是转为建设用地(6.33 km²), 永久性冰川雪地主要是转为以高山流石坡为主的其他用地(9.27 km²)。从转入来看, 人工牧草地转入面积最大(24.61 km²), 主要来自于中覆盖度草地和低覆盖度草地, 分别为17.52 km²和7.09 km²；其次是中覆盖度草地、低覆盖度草地和高覆盖度草地, 转入面积分别为20.30 km²、17.63 km²和16.93 km²；其中, 中覆盖度草地的主要转入来源为低覆盖度草地(13.14 km²)和沼泽地(6.09 km²), 转入高覆盖度草地的主要土地利用类型是中覆盖度草地(15.31 km²), 低覆盖度草地的主要转入来源为其他用地(9.45 km²)和沼泽地(5.55 km²)。总之, 园区土地资源结构和景观格局稳定, 土地利用开发强度低, 属于自然生态系统集中分布区。林地、草地、湿地等生态系统是祁连山重要生物生态源地, 构筑起祁连山“山水林田湖草沙冰”生态安全格局, 维系着祁连山山地生态系统稳定。

2.2 生态系统服务价值时空特征 2.2.1 生态系统服务价值估算

根据表 2可知, 1998年、2008年、2018年园区生态系统服务价值分别为274.36×10⁸元、273.52×10⁸元和273.44×10⁸元, 总体变化幅度不大。从土地资源价值来看, 水体与湿地生态系统的ESV贡献量最大, 占ESV比例在44.00%以上, 分别为水域＞沼泽地＞永久性冰川雪地；1998年、2008年、2018年水域ESV分别为73.50 ×10⁸元、73.29×10⁸元和73.18×10⁸元, 分别占ESV的26.79%、26.79%和26.76%；沼泽地和永久性冰川雪地也是水体与湿地生态系统ESV的重要组成部分, ESV贡献量分别为14.76%—14.97%和2.50%—2.54%, 其中永久性冰川雪地ESV呈现持续下降趋势, 从1998年的6.95×10⁸元减少到2018年的6.84×10⁸元, 减少了1.64%。草地生态系统的ESV贡献量在43.00%以上, 其中高覆盖度草地ESV最高, 3个时期分别为65.40×10⁸元、65.36×10⁸元和65.53×10⁸元；中覆盖度草地ESV在31.93×10⁸元至32.22×10⁸元之间, 低覆盖度草地ESV在21.59×10⁸元至21.69×10⁸元之间, 分别占ESV的11.68%—11.78%和7.89%—7.93%；2018年人工牧草地ESV为0.29×10⁸元, 为黑土滩退化草地治理修复工程产出的最直接生态效益。森林生态系统也是园区生态系统的重要组成部分, 其ESV贡献量在10.18%—10.21%, 其中1998年、2008年、2018年有林地ESV分别为15.42×10⁸元、15.41×10⁸元和15.43×10⁸元；灌木林ESV在12.50×10⁸元至12.51×10⁸元之间, 占ESV的4.56%—4.57%。旱地、其他用地的ESV贡献量较小, 分别为0.11×10⁸元和5.64×10⁸元。总之, 湿地、草地、林地构成了园区ESV的主体, 不同覆盖度草地结构优化对于提升祁连山ESV潜力较大。

从单项ESV来看(表 3), 单项ESV可以归结为三个梯级, 水文调节、气候调节属于第一梯级, 占ESV比例高于57.00%, 1998年、2008年和2018年水文调节价值分别为113.68×10⁸元、113.18×10⁸元和113.06×10⁸元, 气候调节价值分别为44.52×10⁸元、44.50×10⁸元和44.52×10⁸元；生物多样性、水土保持、环境净化、气体调节、美学景观、水资源供给属于单项ESV的第二梯级, 占ESV比例高于38.00%, 依次为生物多样性＞水土保持＞环境净化＞气体调节＞美学景观＞水资源供给；原料生产、食物生产、养分循环是单项ESV的第三梯级, ESV在10.00×10⁸元左右, 占比不到4.00%。祁连山作为“中国湿岛”, 水文调节、气候调节、水资源供给等生态服务功能重要, 为河西走廊、柴达木盆地、河湟谷地等周边地区提供了稳定的生态产品, 生态系统服务溢出价值巨大。

从单位面积ESV来看(图 4), 1998—2018年园区单位面积ESV为172.94万元/km²。水体与湿地生态系统单位面积ESV最高, 其中水域为1586.02万元/km², 沼泽地为656.83万元/km², 永久性冰川雪地为129.68万元/km²；森林生态系统单位面积ESV较高, 有林地为221.34万元/km², 灌木林为192.18万元/km²；高覆盖度草地、中覆盖度草地、低覆盖度草地、人工牧草地单位面积ESV分别为159.89万元/km²、129.44万元/km²、98.98万元/km²和116.61万元/km²；农田生态系统(旱地)单位面积ESV为50.63万元/km²；荒漠生态系统(裸岩石质地、裸土地、戈壁等)单位面积ESV为13.89万元/km²。总之, 河流水系、沼泽湿地、高寒草甸、山地森林是发挥祁连山“气候湿岛”、“绿色半岛”生态服务功能的重要载体, 在水源涵养、气候调节、水资源供给等方面生态服务功能价值高。人工牧草地、山区旱地等人工生态系统在生产有机质、调节空气质量等方面的价值较高, 其脆弱性、敏感性等级高, 适应气候变化能力弱, 土壤保持、生物多样性等方面的生态服务价值低。园区西北部的高原荒漠生态系统, 以裸岩石质地、裸土地、戈壁为主, 与柴达木盆地东缘相连, 干燥少雨、植被稀疏、生产力较低, 生态系统服务能力较低。

2.2.2 生态系统服务价值空间可视化

基于网格尺度实现了生态系统服务价值空间信息的可视化、定量化和精细化, 直观映射了园区ESV空间格局。利用自然断点法, 结合研究区土地资源空间分布格局, 将园区ESV设定为低价值(ESV＜0.40×10⁶元)、较低价值(0.40×10⁶元≤ESV＜0.80×10⁶元)、中等价值(0.80×10⁶元≤ESV＜1.60×10⁶元)、较高价值(1.60×10⁶元≤ESV＜2.40×10⁶元)、高价值(ESV≥2.40×10⁶元)5个等级, 最终得到不同时期单元网格ESV空间格局分布图(图 5)。园区ESV空间异质性明显, 不同ESV等级呈现“大分散、小集聚”镶嵌交错分布格局。具体来看, 高价值区集中分布在黑河源、大通河源、托勒河源、疏勒河源等高寒河源湿地区域, 河流水系、沼泽草甸等水源涵养、水资源供给功能显著, 生态系统重要性等级高, 可提供的价值量、功能量也相应较高。3个时期均为81.00%以上的水域和93.00%以上的沼泽地以高价值等级占主导；较高价值区主要集中分布在园区东南部, 主要为青海云杉、祁连圆柏为主的森林生态系统, 属于稳定的高生产力生态系统, 生态服务功能类别齐全, 生态价值也相对较高；中等价值区在公园分布广泛且较为分散, 呈现为在高价值区和较高价值区外围分布特征, 表现出空间过渡性, 主要是因为草地ESV在中等价值区间分布比例较高所致。从土地资源结构来看, 森林和草地生态系统以中等价值及以上等级为主, 3个时期有林地较高价值及高价值面积占比约为88.90%, 灌木林中等价值及较高价值面积比例在81.00%以上；高覆盖度草地、中覆盖度草地、低覆盖度草地的中等价值及以上等级面积比例分别在93.80%、90.30%、87.40%左右；永久性冰川雪地(岗格尔肖合力雪山冰川、疏勒南山冰川群、冷龙岭冰川等)为中等价值集中分布区, 冰川作为祁连山重要的景观类型, 在淡水资源供给、气候与径流调节、美学观赏与游憩等方面具有较高的显示度, 生态价值不容忽视；裸岩石质地、裸土地、戈壁等其他土地利用类型以低价值、较低价值等级为主, 面积占比在79.00%以上, 高寒荒漠、高山流石坡、戈壁荒滩等区域植被稀疏, 生态脆弱, 生态服务价值也较低。

2.3 景观生态风险空间分异特征 2.3.1 景观生态风险数量特征

1998年、2008年、2018年园区景观生态风险指数(ERI)分别为0.2287、0.2286和0.2310, 生态安全状态整体较好(表 4)。从土地利用类型来看, 3个时期景观生态风险变化特征大体可归纳为4类：旱地、有林地、中覆盖度草地、低覆盖度草地、水域ERI呈波动上升趋势, 沼泽地、永久性冰川雪地ERI持续上升, 灌木林、高覆盖度草地ERI呈波动下降趋势, 建设用地、其他用地ERI持续下降。

可以发现, 土地利用类型ERI量级可以归结为三个梯级：旱地、建设用地ERI在1.0以上, 位于景观生态风险的第一梯级, 受人为生产、生活干扰强度较大, 生态风险等级最高；水域、永久性冰川雪地、沼泽地ERI分别为0.5421—0.5449、0.3415—0.3434和0.3004—0.3016, 属于景观生态风险第二梯级, 受自然扰动和人为干扰会的影响, 水体与湿地生态系统脆弱性较高, 稳定性较低, 局部生态功能退化；灌木林、其他用地、低覆盖度草地、中覆盖度草地属于景观生态风险第三梯级, 依次为灌木林＞其他用地＞低覆盖度草地＞中覆盖度草地, 草地灌丛化现象正在改变原有景观格局, 可能是灌木林景观生态风险变化的诱因, 以高寒荒漠生态系统为主的其他用地, 受人为因素干扰小, ERI也较低；有林地、高覆盖度草地ERI分别为0.1643—0.1645和0.1787—0.1789, 生态系统结构稳定, 自我修复、抗干扰能力较强, ERI最低, 属于第四梯级。

从景观生态风险贡献权重来看, 沼泽地、高覆盖度草地的景观生态风险贡献权重在20.00%以上, 对园区生态安全格局影响较大；其次是中覆盖度草地和低覆盖度草地, 景观生态风险贡献权重为13.66%—13.88%和13.12%—13.27%, 对于优化园区生态安全格局潜力较大；水域、其他用地、永久性冰川雪地、灌木林和有林地的景观生态风险贡献权重较小, 但也是园区景观生态格局构建不可忽视的生态资产；旱地的景观生态风险贡献权重最低, 农业用地开发利用规模不大, 不是主导生态风险源, 但其ERI呈波动上升趋势, 说明土地农业开发存在潜在隐患, 设置开发生态红线, 促进耕地集约利用, 可降低生态环境风险压力。

2.3.2 景观生态风险结构变化

采用自然断点法为将园区景观生态风险等级设定为低风险(ERI＜0.20)、较低风险(0.20≤ERI＜0.25)、中风险(0.25≤ERI＜0.30)、较高风险(0.30≤ERI＜0.35)、高风险(ERI≥0.35)5个等级。进一步统计园区不同等级景观生态风险面积发现(表 5), 1998—2018年低生态风险区面积呈现先增后降的变化趋势, 2018年低生态风险区面积最小, 占园区面积的49.43%；1998年、2008年、2018年较低生态风险区面积分别为6409.24 km²、6405.91 km²和6388.06 km², 呈现波动降低趋势, 30年降低了0.33%；中生态风险区、较高生态风险区面积呈现先降后增的变化趋势, 均为2008年的景观生态风险面积最小, 分别占园区面积的6.58%和2.32%；高生态风险区面积呈现波动上升趋势, 由1998年的143.53 km²增加至2018年的198.34 km², 增加了38.19%, 黑土滩退化草地综合生态修复工程的实施, 增加了人工牧草地种植面积，提高了天然草场的覆盖度, 提升了退化草地原料生产供给功能, 其抵御自然灾害的能力弱, 加之人工牧草地生态稳定性有待提升, 在一定程度上造成了高生态风险区面积扩张。

从生态风险等级结构来看, 3个时期园区ERI以低生态风险等级和较低生态风险等级占主导地位, 分别占园区面积的90.11%、90.21%和89.78%；中生态风险等级、较高生态风险等级面积比例较小, 分别是6.58%—6.65%和2.32%—2.33%, 是园区改善生态安全格局的突破口；3个时期高生态风险等级面积分别仅占园区面积的0.91%、0.89%和1.25%, 在人为和自然因素的共同影响下, 高生态风险等级面积增加可能对于园区生态安全格局产生不利后果。总之, 园区景观生态安全格局以低生态风险等级和较低生态风险等级占主导地位, 但其面积的波动降低和高生态风险等级面积的波动增加, 促使园区ERI存在增加的可能性, 是园区景观生态安全格局优化的潜在隐患。

2.3.3 景观生态风险空间图谱特征

依据生态风险等级阈值区间划分结果得到3个时期园区景观生态风险空间图谱信息(图 6)。园区景观生态风险呈现低生态风险区、较低生态风险区广泛连片式分布, 中生态风险区插花随机式分布, 较高生态风险区、高生态风险区零星集中式分布的图谱特征。具体来看, 3个时期高覆盖度草地、中覆盖度草地、低覆盖度草地均以低生态风险等级为主, 分别占其面积的73.00%、65.00%和54.00%左右；3个时期有林地的低生态风险等级面积占比分别为83.81%、84.07%和83.85%；1998—2018年灌木林、其他用地以较低生态风险等级为主, 面积占比分别为63.28%—63.66%和66.12%—66.74%；1998年、2008年、2018年水域的较低生态风险等级面积占比分别为48.97%、49.23%和48.46%；3个时期沼泽地以较低生态风险等级和中生态风险等级为主, 面积占比分别为88.35%、88.32%和83.34%；1998年、2008年、2018年永久性冰川雪地以中生态风险等级和较高生态风险等级为主, 分别占其面积比例的79.51%、79.48%和79.37%；3个时期旱地以高生态风险等级为主, 分别占其面积的86.14%、86.68%和86.62%；2018年人工牧草地中面积占比的98.32%为高生态风险等级。从典型生态功能区来看, 园区东南部的仙米国家森林公园以低生态风险、较低生态风险等级为主, 园中中部的黑河源国家湿地公园以较低生态风险、中生态风险等级为主, 岗格尔肖合力雪山冰川、疏勒南山冰川群(其中岗纳楼冰川最大)、冷龙岭冰川以中生态风险、较高生态风险等级为主。

2.4 生态空间分区构建 2.4.1 生态空间分区格局

通过对ESV、ERI进行Z-score标准化计算, 得到园区生态分区象限分布图(图 7)、生态分区空间图谱(图 8)和生态分区统计值(表 6), 结果表明：

(1) 园区以低生态系统服务价值—低生态风险区(III)为主, 1998年、2008年、2018分别占园区面积的52.04%、52.05%和52.65%, 单位面积ESV在89.69—89.97万元/km²之间, ERI分别是0.1928、0.1928和0.1931。该生态分区在园区广泛分布, 以中西部地区较为集中, 主要土地利用类型为高覆盖度草地、中覆盖度草地、低覆盖度草地和其他用地, 分别占其面积的44.41%、60.72%、70.86%和78.15%左右。

(2) 高生态系统服务价值—高生态风险区(I)面积呈现下降趋势, 3个时期分别为3035.45 km²、3030.82 km²、2883.58 km², 单位面积ESV分别是401.55万元/km²、400.23万元/km²和407.97万元/km², ERI分别是0.2483、0.2485和0.2556。该生态分区呈现“组团式”或“廊道式”分布格局, 与河流水系、沼泽湿地等重叠度较高, 水域面积占比的74.34—76.16%和沼泽地面积占比的91.07—91.98%位于此生态空间分区, 旱地面积占比中超过62.00%的区域分布在该生态空间分区。

(3) 高生态系统服务价值—低生态风险区(IV)面积呈现增加趋势, 从1998年的3066.89 km²增加至2018年的3224.76 km², 增加了5.15%, 3个时期单位面积ESV依次是211.01万元/km²、210.28万元/km²和211.87万元/km², ERI分别为0.1849、0.1750和0.1861。在园区东部地区分布集中, 其中有林地面积占比中80.00%以上的区域位于该生态空间分区。

(4) 低生态系统服务价值—高生态风险区(II)面积最小, 呈现下降变化趋势, 3个时期面积分别为1490.23 km²、1481.31 km²和1387.91 km², 下降了6.78%, 单位面积ESV在90.95—91.63万元/km²之间, ERI在0.2597—0.2742之间。该生态分区分布范围与园区内永久性冰川雪地的分布基本吻合, 永久性冰川雪地面积占比中98.00%以上的区域位于该生态空间分区, 主要包括岗格尔肖合力雪山冰川、疏勒南山冰川群、冷龙岭冰川等。

综上, 园区以低生态系统服务价值—低生态风险区(III)为主, 面积呈现略增加趋势, 高生态系统服务价值—高生态风险区(I)、低生态系统服务价值—高生态风险区(II)面积呈现下降趋势, 高生态系统服务价值—低生态风险区(IV)面积呈现增加趋势, 园区生态分区结构呈现稳中向好的变化趋势。

2.4.2 生态分区管控方案

根据网格单元的生态系统服务价值、景观生态风险特征, 可以将园区划分为生态保障型生境修复区、生态脆弱型特别保护区、生态改良型发展利用区和生态预防型保育涵养区(图 9)。具体如下：

(1) 生态保障型生境修复区(I)：该类区域单位面积生态系统服务价值高, 生态风险等级高, 是园区生态服务功能重要保障供给区和生境修复优先区, 主要以黑河流域、疏勒河流域的水体与湿地生态系统为主, 属于高山河源型沼泽湿地。采取封闭管理, 严格禁止放牧, 减少人类活动干扰, 逐步提升景观整合度, 降低生境破碎度, 确保以湿地为栖息地的生物生存空间、生态廊道提质优化；适时启动实施流域综合治理工程, 采取自然恢复和人工修复相结合的生态治理模式, 提高水源涵养、生物多样性保护等生态功能；全面建立河长制, 任命河长, 设立河长办公室, 运用科技巡查、牧民巡查力量, 强化河道监管巡查, 确保河流湿地生态系统健康。此外, 在园区东部零星分布的旱地(以种植油菜为主)也以高价值—高风险为主, 重视生态保护, 兼顾农牧业生产, 平衡耕地资源开发强度, 发挥国家公园社会服务功能, 延伸休闲与游憩、科研、科普教育等功能, 实现游憩展示与传统利用、严格保护与合理开发的“双赢”。

(2) 生态脆弱型特别保护区(II)：该类区域单位面积生态系统服务价值低, 生态风险等级高, 是园区典型生态脆弱区, 也是特殊的高山积雪冰川带, 主要包括岗格尔肖合力雪山冰川、疏勒南山冰川群、冷龙岭冰川等高山冰雪带和高寒荒漠、高山流石坡等植被稀疏地区, 属于气候变化敏感区, 生态系统脆弱, 景观脆弱度较高, 内部稳定性较差。此类区域应加强严格保护力度, 坚决杜绝人为干扰, 保护生态系统的原真性和完整性；建立冰川变化监测站点, 开展科学研究提升冰冻圈与气候变化理论认知水平及预警能力, 明确冰川是公共生态产品的属性和冰川保护的全球性责任, 从全球及国家尺度上进行关键领域、关键产业节能减排, 助力“双碳”目标, 减缓全球变暖。此外, 在国家公园外围周边区域谋篇布局科教游憩延伸区, 启动建立国家公园科普馆、祁连山生态文化体验中心, 通过实景监测、远程展示平台, 开展冰川、湿地、野生动物等生态气象科普及环境教育活动, 提升民众科学素养和环境保育意识。

(3) 生态改良型发展利用区(III)：该类区域单位面积生态系统服务价值低, 生态风险等级低, 覆盖范围广, 是园区全面改善、整体提质的生态改良区, 也是生态畜牧业发展利用区, 以高寒草甸和荒漠草原生态系统为主。作为生态—生产协同发展区, 持续推进祁连山地区生态修复型飞机人工增雨(雪)作业综合生态保护和治理举措, 从整体上改善园区土壤湿度, 增加植被生物量, 提升植被覆盖度, 筑牢祁连山生态安全屏障；兼顾草地生态系统的生态、生产功能, 在保护和利用中采取主次分明、划区管理的“组团式”生态治理模式, 低覆盖度草地以生态修复为主, 局部退化严重草地采取工程措施(围栏封育+人工补播)与自然恢复相结合的方式进行；高覆盖度草地、中覆盖度草地坚持发展高质量生态有机畜牧业, 建立季节性放牧制度和轮牧制度, 确保草蓄动态平衡和草地生态服务功能可持续供给。

(4) 生态预防型保育涵养区(IV)：该类区域单位面积生态系统服务价值高, 生态风险等级低, 生态资源本底好, 是园区重要的生态安全源地, 也是重要的保育涵养区。该区域以森林生态系统为主, 主要集中分布在园区中东部, 包括有林地和灌木林, 是园区调节服务、支持服务生态功能的重要支撑区。以青海云杉、祁连圆柏等为主的寒温带针叶林分布较为集中, 可通过划定生态保护“特区”采取“集中连片式”预防型生态保护模式, 确保生态涵养功能保值增值, 突出生态系统服务功能的溢出效应；加强森林火灾监测预警, 保证森林资源更新和再生, 促进林地资源合理利用和良性循环。

3 讨论

本研究中土地资源本底现状表明草地资源是祁连山国家公园青海片区的主体景观类型, 与付建新等^[36]的研究结果相互印证, 其研究表明2018年祁连山国家公园青海片区草地面积为8174.93 km²。需要指出的是, 祁连山国家公园青海片区地处青藏高原东北缘, 受自然海拔、气候因子的影响, 草地构成了高寒植被生态系统的主体, 表现出明显的高原地带性特征。本研究核算出祁连山国家公园青海片区单位面积ESV为172.94万元/km², 与祁连山南坡单位面积ESV为76.27万元/km²的研究结果^[37]相互补充, 折射出祁连山南坡ESV空间分异特征明显, 祁连山国家公园青海片区为ESV的高值区, 其蕴藏的生态价值远高于祁连山南坡其他地区, 生态系统服务价值盈余潜力巨大, 是国家生态功能格局中提供公共生态产品的关键区域, 对维系祁连山地区乃至中国北方地区生态安全有着不可获取的作用^[13—14]。与此同时, 本研究中永久性冰川雪地单位面积ESV为129.68万元/km², 与孙美平等^[38]关于祁连山区单位面积冰川服务价值的研究结果存在差异。一方面, 由于数据来源、指标选取、评价方法的限制, 冰川生态价值评估研究之间可比性存在一定的局限性, 另一方面, 本研究估算的永久性冰川雪地ESV只是一个保守值, 但依然可为掌握冰川服务价值提供重要的生态信息和决策依据。由此可见, 生态系统服务价值核算具有相当大的复杂性和不确定性^{[14, 39]}, 在谋求探究一个相对准确的价格标签的同时, 重视生态系统服务价值的边际变化同等重要^[40]。

景观生态风险定量化评价是衡量区域可持续发展的命题。本研究中祁连山国家公园青海片区ERI在0.2286—0.2310之间, 生态安全状况总体稳定, 景观生态风险等级以低生态风险和较低生态风险为主, 这与青藏高原大部分区域景观生态风险处于稳定状态的研究结果一致^[41]。祁连山国家公园青海片区景观生态安全状态整体较好, 从一个侧面也验证了祁连山生态保护与综合治理工程实施的生态效益较好。林地、草地生态系统稳定, 抗干扰能力强, 景观生态风险等级低, 是重要的生态安全源地^[42]。近年来国家公园生态保护分区的相关文献^{[2, 27, 43—45]}围绕国家公园管理目标开展了前瞻性、尝试性研究, 形成了针对性较强、易操作实施的国家公园科学管理先导性经验。目前, 从价值—风险综合角度开展生态功能空间识别还处于不断探索阶段, 研究路径尚不明晰, 可借鉴的实证案例较少, 评价方法还有很大改进和提升空间。在未来, 还需继续深入总结、归纳、完善生态分区的技术方法, 以便更好服务于国家公园建设有序推进。

本研究集成遥感技术、地理信息模型方法、景观生态学方法、GIS网格法, 结合生态系统服务和景观生态风险两个维度, 阐释土地利用、生态系统服务价值及景观生态风险的时空演变特征, 提出园区生态分区差异化管控方案, 较好地识别了生态资源空间信息, 为自然保护地生态功能分区空间识别研究提供了新思路。必须指出, 祁连山作为西北干旱半干旱地区的“生态半岛”和“气候湿岛”, 生态系统服务功能外溢效应显著, 祁连山地区与河西走廊、柴达木盆地、青海湖流域、河湟谷地等共同构成了“生态地理单元综合体”。本研究仅开展了祁连山国家公园青海片区生态系统服务价值、景观生态风险评估及生态空间分区识别, 未从“生态地理单元综合体”视角综合考虑祁连山山地系统生态价值的溢出效应、溢出边界、生态损益、盈余潜力及风险阈值等, 这应是今后进一步探索的方向。

4 结论

(1) 草地是园区最主要的土地覆被类型, 构成了高寒植被生态系统的主体, 占园区面积的55.00%以上；30年间园区土地利用之间转移总面积为102.49 km², 中覆盖度草地转出面积最大(34.83 km²), 人工牧草地转入面积最大(24.61 km²)。

(2) 1998—2018年园区ESV分别为274.36×10⁸元、273.52×10⁸元和273.44×10⁸元；湿地、草地、林地是构成了园区ESV的主体, 其中水域与湿地生态系统的ESV的贡献量最大, 占ESV比例在44.00%以上；园区单位面积ESV为172.94万元/km², ESV空间异质性明显, 高寒河源湿地区(黑河源、大通河源、托勒河源、疏勒河源等)和园区东南部寒温带针叶林区(以青海云杉、祁连圆柏为主)为ESV的高值区。

(3) 3个时期园区ERI分别为0.2287、0.2286和0.2310, 生态安全状态整体较好；3个时期园区景观生态风险以低生态风险等级和较低生态风险等级占主导地位, 分别占园区面积的90.11%、90.21%和89.78%；园区ERI呈现较低生态风险区、低生态风险区广泛连片式分布, 中生态分区区插花随机式分布, 较高生态风险区、高生态风险区零星集中式分布的图谱特征。

(4) 园区以低生态系统服务价值—低生态风险区(III)为主, 在园区广泛分布, 以中西部地区较为集中；高生态系统服务价值—高生态风险区(I)呈现“组团式”或“廊道式”分布格局, 与河流水系、沼泽湿地重叠度较高；高生态系统服务价值—低生态风险区(IV)面积呈现增加趋势, 在园区东部分布集中；低生态系统服务价值—高生态风险区(II)面积最小, 分布范围与园区内永久性冰川雪地基本吻合。