生态系统服务作为人类赖以生存的自然环境与效用, 是实现生态安全的前提保障, 其价值的损益直接影响到区域可持续发展水平^[1-2]。随着我国城镇化进程提速, 生态环境问题日益突出, 重要生态系统服务不断下降, 人类活动对生态系统的强烈干扰逐渐成为保障区域生态安全和实现区域可持续发展的最大威胁^[3]。在此背景下, 基于“反规划”思想和“格局-过程”互馈理论的生态安全格局(Ecological Security Patterns, ESP), 以提升生态系统服务与人类福祉为目标, 提出改善区域生态问题的对策与措施并落实于空间地域, 被视为缓解当前生态保护与土地开发冲突、保障区域生态安全的有效途径^[4]。因此, 以维持生态系统服务为基础的生态安全格局构建, 已成为当今地理、生态以及城市规划等众多学科共同关注的热点^[5]。

生态安全格局研究源于20世纪60年代Warntz^[6]提出的物种基于潜在表面流动过程建立的点线面模型。自20世纪90年代开始, 以俞孔坚等为代表的国内学者对生态安全格局构建展开研究, 并将成果应用于城市边界增长预测与城市空间发展预警等方面^[7-8]。综合来看, 经过二十余年的发展, 生态安全格局构建的理论基础与方法日趋成熟, “源地识别-廊道提取”的组合方式已成为构建生态安全格局的基本范式^[9]。生态源地提取作为生态安全格局构建的基础, 国内研究多趋于选取具有既定边界与面积的自然保护区^[10]、绿地或水体^[11], 或建立综合评价指标体系识别斑块重要性^[12]。然而生态系统服务是提高人类福祉和实现区域可持续性的基础, 各生态系统服务之间相互作用、关系复杂, 已有研究鲜有从多种生态系统服务价值(Ecosystem Services Value, ESV)综合评估角度进行源地与廊道识别以及生态系统空间格局优化。目前, 国内学者多引用谢高地等^[13]制定的单位面积生态系统服务价值当量因子表核算生态系统服务价值^[14-15], 这种方法易于操作、结果便于比较^[16], 但其研究视角是基于全国平均状态, 如果应用于小尺度区域研究, 估算结果难以反映研究区实际情况。基于植被净初级生产力(Net Primary Productivity, NPP)的遥感估算模型, 克服了传统生态系统服务价值估算方法难以反映真实空间分布状况的缺点, 可以直观地判断不同地带在生态系统服务上的贡献大小^[17], 从而为土地利用优化配置决策的制定提供依据。

瓦房店市地处东北亚经济圈节点位置, 是大连市联系辽中城市群以及东北经济走廊的必经之地, 同时也是大连市北部重要的生态屏障。面对发展经济和保护生态的双重选择, 科学评估瓦房店市ESV损益并进一步构建合理的生态安全格局, 对于缓解城市扩张与生态保护的矛盾具有重要意义。基于此, 本文尝试以遥感估算模型结合生态经济法评估2000—2014年瓦房店市九项ESV损益情况, 并以此为基础构建两期生态安全格局, 通过对比分析各组分的演变来明确瓦房店市生态屏障的退化轨迹, 以期为研究区土地利用开发和生态系统服务提升提供科学参考。

1 研究区概况与数据来源 1.1 研究区概况

瓦房店市位于辽宁省南部, 辽东半岛中西部(图 1, 121°13′E—122°16′E, 39°20′N—40°07′N), 西濒渤海, 总面积为3881 km²。瓦房店市地处温带季风气候区, 兼有海洋性特征。全市地形由东北山地向西南沿海倾斜, 形成了山地、丘陵、平原与滨海湿地相结合的地貌类型, 蕴含多种生态系统服务。瓦房店市凭借其优越的地理位置, 雄厚的工业基础, 成为辽宁省县域经济的排头兵, 享有“中国轴承之都”的美誉。然而, 海陆兼备的地缘优势带来巨大发展机遇的同时, 也让瓦房店面临双重生态问题。一方面, 随着辽宁沿海经济带建设步伐的加快, 瓦房店市大力推进全域城市化战略, 生态压力逐年增大, 滨海滩涂等自然生态系统被侵占现象严重；厂矿向海排污造成近岸海域水质恶化；海水养殖的无序发展逐渐打破海洋原有的生态平衡。另一方面, 陆域生态问题同样严峻, 大量承接大连市重化企业转移已埋下空气质量安全隐患；东部山区水土流失问题不断加剧；主要河流季节性断流更是加剧了本已严峻的用水紧缺。瓦房店市迫切需要统筹陆海, 综合评估全域生态系统服务水平, 构建海陆联结互通的生态安全格局。

1.2 数据来源与处理

本研究两期遥感影像来源于美国地质调查局网站(https://www.usgs.gov/), 分别为Landsat7 ETM(2000年9月16日)和Landsat8 OLI(2014年9月8日), 空间分辨率30m。借助ENVI 5.1, 结合Google地图高分辨率影像进行人机交互解译, 结果经混淆矩阵精度检验, kappa值均达到0.85以上, 符合研究要求。

NDVI数据来源于美国航空航天局网站(https://www.nasa.gov/)MODIS13Q1数据集的16天合成产品(2000年1—12月, 2014年1—12月), 空间分辨率250 m。经双线性内插后将NDVI分辨率精度提升至30 m, 以满足研究需求。通过最大值合成法将数据两两合成, 生成逐月NDVI数据, 使其保持与气象数据一致的时间分辨率。

气象数据来自中国气象数据网(http://data.cma.cn/), 获取瓦房店市以及周边13个站点的逐月降水量、蒸发量、月平均温度、太阳总辐射和净辐射等数据。利用ArcGIS 10.2对各气象数据进行反距离权重空间插值, 获得空间分辨率为30 m的气象空间分布栅格数据集。

其他数据主要包括数字高程模型(DEM), 土壤数据和社会经济统计数据。研究区DEM数据来源于地理空间数据云(http://www.gscloud.cn/), 空间分辨率为30 m。土壤数据来自1:100万土壤数据库。本研究使用的社会经济统计数据主要来自《大连市统计年鉴2015》和《瓦房店市统计公报2015》。

2 研究方法 2.1 NPP估算方法

植被净初级生产力是衡量植被生产能力和覆盖度的重要指标^[18], 用来表示单位时间内植物在单位面积上所累积的有机物能量, 即植被在光合作用过程中扣除自氧呼吸后创造出的有机物质剩余量^[19]。本研究采用朱文泉等^[20]改进的CASA(Carnegie-Ames-Stanford Approach)模型对NPP进行估算。具体计算公式为：

(1)

式中, NPP(x, t)为植物t月份在像元x处的有机物质累积总量(g C m^-2 a^-1)；APAR(x, t)为t月份像元x所吸收的有效光合辐射(MJ m^-2 a^-1)；ε(x, t)为植物t月在像元x处的光能实际利用率。

(2)

式中, SOL(x, t)表示太阳t月份在像元x处的总辐射量(MJ m^-2 a^-1)；FPAR(x, t)为植被t月份在像元x吸收有效光合辐射的百分比；固定值0.5表示植被可利用有效辐射与总辐射之比。

(3)

式中, ε_max为最大光能利用率(%)；T_ε1(x, t)和T_ε2(x, t)分别表示低温和高温胁迫影响系数；W_ε(x, t)为水分胁迫影响系数。

2.2 ESV估算方法

本研究参考Costanza等^[21]的生态系统服务分类体系与千年生态系统评估(MA)^[22], 并结合研究区生态环境特点, 将瓦房店市生态系统服务划分为4大类和9小类(表 1)。目前应用较为广泛的当量因子法通过专家评分方式建立生态服务当量的价值关系, 这种关系由于受专家经验影响, 具有一定局限性^[23]。基于NPP的遥感估算模型技术先进, 实用性强^[24], 不仅实现了ESV在时空尺度上的连续性, 而且结合典型样点的实测数据进行“本土化”修改, 能更好地体现研究区本土生态系统特征。因此本研究选取遥感估算模型、当量因子法和生态经济法相结合的估算方法, 对研究区2000年与2014年ESV进行测算。各项ESV的估算模型详见表 1, 各年份的ESV总价值计算公式如下：

(4)

式中, i表示各项生态系统服务类型, 其中V₁为有机物质生产价值, V₂为气体调节价值, V₃为气候调节价值, V₄为水源涵养价值, V₅为土壤保持价值, V₆为废物处理价值, V₇为养分循环价值, V₈为生物多样性保护价值, V₉为文化服务价值, 单位：元m^-2 a^-1。

2.3 生态安全格局构建方法

本研究依据景观安全格局“廊道-斑块-基质”理论^[33], 借鉴前人研究成果确定了由生态源地、廊道、缓冲区等组分构成的生态安全格局。

2.3.1 生态源地识别

生态源地作为生态系统服务流动的起点和基地, 是整个区域ESV的高值集聚之地, 其完整性和准确性对区域生态过程至关重要^[8]。本研究运用ArcGIS 10.2中Getis-Ord Gi^*模块对2000年和2014年两期生态系统服务总价值进行冷热点分析。Getis-Ord Gi^*模块用以判别区域高值(热点)和低值(冷点)的空间聚类。通过冷热点分析, 提取两期ESV总价值90%置信度以上的热点集聚区作为瓦房店生态安全格局的源地。同时, 为了保证生态源地的完整性以及维持生态源地正常的辐射功能^[34], 本研究剔除面积较小且分布零散的碎小斑块, 选取热点聚集区内面积在1 km²以上的斑块作为最终生态源地。

2.3.2 阻力面设置

本研究采用Knaapen^{[7, 35]}等提出的最小累积阻力模型(Minimum Cumulative Resistance, MCR)建立生态系统服务流动的阻力面。MCR模型通过计算物种从源地到目的地过程中所需代价来模拟物种的移动路径, 其计算公式为：

(5)

式中, MCR为最小累积阻力值；f为从i处到j处的最小累积阻力与生态系统服务流动过程的正相关函数；D_ij表示物种从源地j移动到目的地i的距离, R_i为物种移动到目的地i的阻力系数。

本研究依据最小阻力模型, 分别将2000年与2014年两期单位面积ESV取倒数, 以取倒数后的单位面积ESV作为生态系统服务流动的“成本”, 从而确立阻力面, 即阻力值随ESV升高而减小。

2.3.3 生态廊道及其他组分的提取

生态廊道可以将散布的生态斑块连接起来^[20], 是物种在源间移动的通道, 也是生态系统服务流动的通道。识别关键生态廊道并对其加以保护对保障生态要素正常流动具有重要意义。本研究依据最小累积阻力模型的阻力面设置, 借助CIRCUITSCAPE软件识别源地间的最小耗费路径, 进而获得生态廊道的空间分布。此外, 基于MCR模型以及最小累计值与面积关系曲线的突变情况, 以阻力阈值确定生态安全分区, 将研究区除生态源地以外的所有区域分为生态缓冲区、过渡区和其他区域。

3 结果分析 3.1 NPP估算结果分析

基于CASA模型获得研究区2000年和2014年NPP变化情况(表 2)与空间分布(图 2)。2000年与2014年瓦房店市的全年NPP总量分别为1.61万亿g C和1.65万亿g C。如表 2, 2000年林地、耕地、海水养殖盐田、其他用地等土地利用类型所累积的NPP总量均高于2014年。但瓦房店市建设用地的面积的骤增, 导致2014年建设用地NPP总量高出2000年约1142亿g C。总体来看, 得益于2014年良好的水热条件, 2014年植被生长状况好于2000年, 使得2014年NPP无论是从总量还是从平均量都稍高于2000年。空间分布上, 研究区两期NPP均呈现出东北高、西南低的分布趋势, 高值区多集中于东部老帽山、龙潭山和中南部大北山等山地林区。低值区多分布在西南沿海地区以及复州镇和瓦房店市区, 这些区域人类活动频繁, 植被覆盖度较低, 且多为人工栽培作物。经过近15年的高强度城市扩张与围填海, 瓦房店市建设用地面积大幅增加, NPP也由高低值均匀分布演化为低值集聚。

3.2 ESV损益分析 3.2.1 ESV时间变化

以表 1及式(4)所示的模型为基础, 利用ArcGIS 10.2栅格计算器计算得到瓦房店市2000年与2014年的ESV总价值表(表 3)以及各项ESV构成变化图(图 3)。2000年与2014年价值总量分别是117.38亿元和95.52亿元, 价值总量下降趋势明显。两期生态系统服务构成中, 水源涵养总价值分别为40.22亿元和22.48亿元, 贡献率分别达34.26%和23.53%, 为各服务类型之首。土壤保持总价值最低, 贡献率仅4.21%和3.33%。就单项ESV变化来看, 2000—2014年瓦房店市除养分循环变化不明显外, 其余各项服务价值均持续降低。其中水源涵养服务损失量最大, 达到17.74亿元, 降幅高达44.07%。其次是生物多样性保护, 降幅为24.75%。其他各项服务价值损失量均小于1.2亿元, 降幅稳定在1%—10%之间。

2000—2014年, 瓦房店市各土地利用类型面积变化显著。如表 4所示, 林地、其他用地和耕地面积呈减少趋势, 其中林地面积减少最为严重, 达到262.56 km²。从空间分布来看(图 4), 14年间研究区东部林地分布范围急剧收缩, 格局趋于破碎化, 西部沿海未利用土地(其他用地)也逐渐被建设用地侵占。受此影响, 瓦房店市林地ESV损失约8.74亿元(表 3), 占到损失总值(21.86亿元)的40%, 为各地类损失价值中最高。而建设用地、海水养殖与盐田面积增加趋势明显, 尤其是建设用地面积在近14年中增加了347.93 km², 土地变化率高达74.71%。如图 4所示, 相比于2000年, 2014年瓦房店市区和西南沿海地区建设用地面积成倍扩增, 挤占林地和耕地现象明显。受此影响, 研究区建设用地ESV增加了63.48%, 但远低于其土地利用面积的增长速度。此外, 草地ESV虽然随其面积增加而有所增长, 但1.57%的价值变化率远不足弥补林地生态价值的损失。可见, 城建区面积不断扩大, 林地遭到持续破坏, 农业用地逐步收缩, 是瓦房店市ESV损益变化的主要原因。

3.2.2 ESV空间分布变化

研究区两期单位面积ESV空间分布均表现为东北和西南高、中部低, 水域(包含海域)高、陆地低的总体趋势(图 5)。高值区多集中于东部和中南部林地面积广大、湖泊水库密布、植被覆盖度高的山地。而低值区域则与建设用地分布格局相一致, 且2000—2014年间, 随着城市化与围填海活动的加剧, 低值分布区呈现急剧扩增趋势, 特别是沿海方向生态系统服务退化现象最为明显。2005年后, 西南沿海长兴岛和北部太平湾两个经济技术开发区相继设立, 瓦房店市兴起一场持久性、高强度的围填海运动, 城市边界快速增长和沿海高强度开发影响了生态系统服务的正常运转, 构建生态安全格局刻不容缓。

3.3 生态安全格局构建及演变分析

利用ArcGIS 10.2软件中Getis-Ord Gi^*模块对瓦房店2000年与2014年两期ESV进行冷热点识别, 得到研究区两期ESV高低值集聚空间分布(图 6)。相比于2000年, 2014年瓦房店市ESV高值集聚区大幅缩减, 除近海海水养殖区域变幅较小之外, 陆地高值集聚区大面积消失, 仅剩河流、水库、湖泊以及东北部一些植被覆盖度较好的山地地区。

在ESV冷热点空间分布基础上, 提取置信度90%以上的热点区并剔除1 km²以下的碎小斑块作为瓦房店市两期生态源地(图 7)。研究区2000年与2014年生态源地面积约为1111.40 km²和356.96 km²(表 5), 分别占研究区总面积的28.63%和9.2%。两期源地均以水域(含近岸海域)和林地为主, 相比于2000年集中连片的分布特征, 2014年生态源地分布离散破碎, 整合度低, 对空间胁迫缓冲能力下降。生态源地作为支撑全市生态系统服务的核心区域, 应加强对原有生态资源保护, 严禁与生态保护无关的建设活动, 守好生态底线。

基于最小阻力模型低阻力值谷线, 提取出瓦房店市两期生态廊道(图 8)。研究区2000年与2014年廊道长度分别为136.63 km和231.43 km(表 5), 长度增加了近1倍, 其主要原因是源地面积收缩加大了生态服务在源间流动的难度。生态廊道整体呈横纵网状分布, 提高了各生态源地之间的空间连通性。生态廊道虽为研究区生态服务流动的低阻力值谷线, 但并不等同于其拥有高水平的生态系统服务就能保证生态“安全”。随着城市扩张导致绿地系统破碎化现象不断加剧, 生态廊道受到人类活动的干扰威胁也将增大。应对其加大保护力度, 按照廊道长度和生态服务流动阻力进行廊道风险评估, 以确定廊道建设的优先次序。对人类活动频繁区域附近和已遭破坏的廊道要列为修复重点, 以恢复其“绿色通道”功能。科学、有针对性地评估复州河、浮渡河流域实施水利工程的可行性, 切实解决二河季节性断流问题, 以恢复水体廊道。

根据最小累积阻力模型中阻力值与面积曲线突变情况, 通过设定阈值确定瓦房店市两期生态缓冲区与过渡区(图 8, 表 5)。生态缓冲区面积分别为816.51 km²和882.15 km², 主要包括湖滨湿地, 山地边缘人工疏林地以及部分海岸带, 该地区土地利用应以生态保护为主, 原则上禁止任何形式的开发活动。生态过渡区面积分别为711.59 km²和1115.78 km², 主要为生态用地与建设用地的交界区。该地区是研究区的限制开发区, 作为潜在的生态后备用地, 应该以保护性开发为主, 严格控制城市开发对生态用地的干扰, 同时应该加强该区生态基础设施建设, 以遏制生态源地持续缩小的趋势。

将上述生态源地、生态廊道、生态缓冲区以及生态过渡区组合叠加, 共同构成了瓦房店市生态安全格局。总体来看, 14年间城市化快速发展以及沿海地区高强度围填海活动, 导致瓦房店市生态用地退化轨迹明显。本文研究结果可为瓦房店市未来人类开发建设活动提供科学警示, 为推进瓦房店市高质量全域城市化的“生态红线”划定提供重要科学参考。

4 结论与讨论 4.1 结论

本文通过建立遥感估算模型对瓦房店市ESV进行核算, 并以此为基础构建生态安全格局。通过对比分析2000年与2014年研究结果来揭示沿海地区在土地利用变化影响下ESV的损益以及生态安全格局演变。结果表明：

(1) 瓦房店市2000年与2014年NPP总量分别为1.61万亿g C和1.65万亿g C, 均呈现出东高西低的分布趋势, 且这种分布不均的趋势2014年更明显。NPP高值区多集中于东部和中南部山地, 低值区则多集中于拥有大面积建设用地的沿海地区和市区。

(2) 2000—2014年, 瓦房店ESV总值减少了约21.86亿元, 降幅达18.12%。ESV空间分布表现为水域高于陆地, 山地高于平原的趋势。14年间, 瓦房店市的高值集聚区面积不断缩小, 低值面积逐渐扩大, 以西南沿海地区和中东部山地最为突出。

(3) 土地利用变化是导致2000—2014年ESV损失的主要原因。快速城市化以及高强度的围填海导致瓦房店市林地、耕地和未利用地面积大幅度减少, 特别是林地减少了约262.56 km²。而建设用地、海水养殖与盐田面积成倍增长, 增加面积达到近500 km²。城市边界快速增长和沿海高强度开发造成了瓦房店市ESV大量损失, 鉴于此, 政府应加强生态保护意识, 慎重决定围填海及土地改造开发实施地点及方向。

(4) 2000—2014年, 生态源地面积减少了近800 km², 受此影响生态廊道长度增加了近一倍。研究区ESV损失造成生态安全格局各组分退化轨迹明显。基于2014年所构建的生态安全格局, 后续城市建设规划工作中应切实转变“向海索地”的工作思路, 禁止审批新增围填海项目, 并加强内陆林地的保护。在东北部老帽山、龙潭山和中南部大北山等林地连片山地以及西部和西南部沿海毗邻海水养殖与盐田的近岸海域划定自然保护区, 以开发与保护并重的视角, 统筹陆海国土空间开发保护, 以期综合提升瓦房店市生态系统服务, 保障区域生态安全。

4.2 讨论

分析生态安全格局时空动态演变, 对明确区域生态退化方向和轨迹、指导生态恢复和建设具有重要意义。通过构建生态安全格局达到对生态过程的有效调控, 已成为目前缓解生态保护与经济发展的重要空间途径。目前基于静态生态安全评价的安全格局构建渐趋成熟, 而基于生态系统服务供需的安全格局构建以及区域生态安全格局的动态变化逐渐成为研究的趋势和学界的共识^[36]。此外, 不同的生态系统服务间的具体相互作用关系复杂^[37-38], 不同时间和空间尺度的权衡与协同关系也不尽相同, 识别生态系统服务间此消彼长的权衡作用与同增同减的协同作用对实现生态系统可持续管理至关重要^[39]。然而在现有研究在区域生态安全格局构建中, 假定不同生态过程之间相互兼容, 彼此之间不存在协同或权衡关系。但事实上, 不同生态过程之间存在密切的相互作用、相互影响, 其与整体生态功能的耦合关联可能大相迥异^[40]。因此, 有必要深入解析不同生态过程的耦合关联, 构建面向山水林田湖生命共同体一体化管理的区域生态安全格局。

本文在研究区域生态安全格局的动态变化过程中, 对安全格局的各组分变化的准确量化仍有不足。此外, 采用基于遥感模型的生态系统服务价值评估方法虽然拥有强大的计算能力和科学性, 但是简单的数字累积和公式模拟并不能完全体现生态安全的内涵, 部分公式当中经验系数也需要结合研究区当地实测数据不断模拟和调试才能获得接近研究区实际水平的生态系统服务。