人口增长、城市扩张、土地利用变化等人类活动对生态系统的结构和功能造成极大干扰^[1], 加之“零风险”环境管理弱点的暴露^[2—3], 使得许多区域面临严重的生态风险。建设用地、耕地等对湿地、水系的侵占趋势加剧了生态用地的减少, 不合理的土地利用变化是生态系统服务价值丧失的主要驱动力, 而破碎化的景观格局导致生态系统内物种病变、栖息地破坏、生物多样性减少等负面影响, 严重威胁着区域生态系统安全^[4—5]。生态风险评价是研究生态系统及其组分由于暴露于多种风险压力源而导致的生态效应可能性^[6—8], 是可持续生态系统管理的有效工具^[9—11], 近年来逐渐成为热点研究话题^[12—16]。

20世纪80年代, 生态风险评价研究开始兴起, 经历了从环境风险、生态风险到区域生态风险的发展历程^[17—18]。生态风险评价研究也由单一风险源、单一风险受体逐渐发展为利用空间要素表征多种风险源和多种风险受体的影响^[19—20], 研究尺度由种群、群落等微观研究发展为区域水平的宏观研究, 研究范围由局地扩展到区域生态风险评价^[21—22]。作为生态风险评价的一个分支, 区域生态风险评价是在区域尺度上综合集成人类活动与自然环境等多种胁迫因子对风险受体造成的影响及危害程度进行描述和评价, 逐渐成为环境科学研究的重要领域与前沿热点^[23—24]。

湿地生态风险评价作为一类区域生态风险评价, 侧重于识别湿地面临的主要风险源(自然灾害、人为因素等)对湿地风险受体(土地利用、景观格局等)可能导致的危害^[25—26]。三江平原是我国最大的淡水沼泽湿地分布区, 但由于日益增长的人口压力与人类活动大面积垦殖与过度开发导致三江平原湿地面积减少、景观破碎化程度高、生态系统的结构和功能遭到严重的干扰破坏^[27], 实现三江平原湿地资源合理利用与有效保护是解决现存问题的基础。从土地利用视角, 一些学者分析了三江平原土地利用变化对区域生态环境的影响。如刘兴土等^[28]揭示了三江平原大面积开荒引发土地退化、动植物资源稀缺等环境变化；杨春霞等^[29]发现政策和人口是三江平原土地利用变化的主要驱动因素, 土地利用变化进而导致湿地生态系统服务功能下降、生物多样性丧失、水环境恶化与农田土壤退化等生态环境问题；张彪等^[30]阐述草地退化、土地盐碱化等土地利用空间格局变化现象导致生态系统基本结构和功能破坏或丧失, 降低了区域生态系统的稳定性。从景观生态视角, 一些学者从景观空间格局出发分析其对三江平原湿地生态过程的维护和控制作用。如刘吉平等^[31]基于三江平原生物地理信息对生物多样性进行规划设计, 通过计算鸟类物种运动阻力指数建立等阻力面, 判别景观生态安全格局；刘春艳等^[32]研究三江平原景观生态风险时空动态变化规律, 量化地形地貌、居民点建设等驱动因子的影响。从评价体系视角, 湿地生态风险评价尚未形成统一体系, 不同学者构建了不同的三江平原湿地生态风险评价体系。如廖玉静等^[33]从三江平原湿地生态环境、湿地服务功能、人类社会影响三方面构建湿地生态系统稳定性评价指标体系, 解释湿地稳定性的成因与机制, 协调湿地保护与农业生产之间的关系；陈红光等^[34]基于风险分析对三江平原水稻灌区建立多水源联合调度模型, 以地下水开采最小为目标, 确定联合调度优化方案；Jiang等^[35]从水资源短缺风险的角度, 运用层次分析法建立评价指标体系和评价标准, 对三江平原水资源短缺风险进行综合评价和区域差异分析；王辉等^[36—39]一系列研究分别从三江平原风险源与风险受体的尺度依存性问题、点轴体系理论、土壤侵蚀及非点源输出分析产生的干扰对湿地风险的直接与间接影响、构建压力-资本-脆弱性-响应模型等方面评价了三江平原湿地区域生态风险。Fu等^[40]基于压力-状态-响应模型分析未来气候变化对三江平原湿地植被生产力与多样性的生态风险等级评价。综合来看, 土地利用变化、景观生态格局等方面的融入丰富了湿地生态风险评价体系, 现有研究虽对空间格局已有分析, 但较少分析三江平原湿地生态风险在空间的集聚现象与辐射扩散特征。本文通过构建三江平原生态风险评价体系, 从空间相关性角度, 确定研究土地利用与生态风险空间关系的最佳阈值距离, 采用双变量空间自相关模型揭示湿地生态风险的空间集聚特征, 对加强三江平原湿地风险防控与区域管理决策具有一定的理论意义。

本文以三江平原为研究区域, 基于2000年、2005年、2010年、2015年4期遥感影像, 从土地利用、道路密度、夜间灯光等方面分析以人类活动为主导的风险源对三江平原湿地的影响。其中, 土地利用与生态风险密切相关, 不同土地利用结构和强度的生态影响可直接作用在区域生态系统变化上, 并且具有典型的区域性和累积性特征；道路密度代表道路建设对生态用地连通性的扰动；夜间灯光数据与人口分布密度呈线性相关关系可用来表征人类活动强度。三江平原湿地生态系统结构和功能受景观、生态系统服务等风险受体的影响显著, 综合考虑景观生态格局、生态系统服务价值两种受体，分析三江平原湿地生态风险强度并运用空间相关分析法剖析生态风险空间差异特征及时间变化特征, 基于土地利用与生态风险双变量选取合适的距离阈值确定空间权重进行空间自相关分析, 从而为降低局部生态风险, 提高三江平原湿地生态风险防控提供理论依据, 实现三江平原湿地生态安全与可持续发展的目标。

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

三江平原位于黑龙江东北部, 地处45°01′—48°28′N, 130°13′—135°05′E之间, 涵盖佳木斯市、鹤岗市及鸡西市等23个县市(如图 1)。作为我国东北三大平原之一, 三江平原是由黑龙江、松花江及乌苏里江冲击而成的低平原, 总面积达10.8万km², 是中国淡水沼泽湿地的集中分布区^[32]。气候类型属于温带湿润、半湿润大陆性季风气候, 夏季温暖而短促, 冬季寒冷漫长, 全年日照时数2400—2500h, 雨热同期, 总降水量的75%—85%集中在夏季, 主要土地类型为耕地、林地、草地、水域、城乡建设用地及未利用土地等。

1.2 数据来源

本研究使用数据主要有三江平原2000年、2005年、2010年和2015年4期遥感影像数据, 来源于地理景观遥感组数据处理, 得到耕地、林地、草地、水域、城乡建设用地、未利用地6个一级地类以及23个二级地类的土地利用/土地覆被数据(如图 2), 自2000年至2015年三江平原土地利用结构发生了显著变化, 2000年耕地面积占总面积的28.34%, 2015年较2010年耕地面积增加了20.11%, 草地面积由占比14.71%减少到占比12.01%, 2015年林地面积由2000年的22.45%削减了将近一倍, 2015年的农村用地面积有所减少但城镇用地面积较2000年扩大4倍, 湖泊、滩地面积分别减少了22.35%、47.52%；多等级道路数据, 来自中国科学院东北地理所数据网(http://www.igadc.cn/wetland/index.html)；DMSP/OLS(Defense Meteorological Satellite Program/Operational Linescan System)夜间灯光数据来自地理国情监测云平台(http://www.dsac.cn/DataProduct)；稻谷、小麦和玉米的价格及播种面积等统计数据出自《全国统计年鉴》、《全国农产品收益资料汇编》；基于ArcGIS 10.2软件平台处理数据, 将研究区划分为5km×5km的单元网格, 共有样区4636个, 计算每一个风险小区内综合生态风险指数。

2 研究方法 2.1 生态风险评价

基于相关研究^[41—45], 本文从风险源与风险受体两个维度分别选取人类活动强度为风险源, 景观生态格局与生态系统服务两方面为风险受体构建生态风险评价体系, 采用如下公式进行评价^[35]：

(1)

式中, R代表综合风险值；H是人类活动风险源强度；EV是景观脆弱性；ESV是各生态系统类型服务价值。

风险源是指可能对生态系统产生不利影响的风险来源, 包括自然风险源与人类活动导致的风险源^[46—49]。本文参考王辉等选取土地利用强度、道路密度、DMSP/OLS夜间灯光分布以及受人为与自然干扰引发的土地转移强度赋予权重0.2、0.3、0.3、0.2^[35]以评价风险源强度, 公式如下：

(2)

式中, w_i表示第i个因子权重, f_i为第i个风险源归一化数值, n为风险源个数。参考庄大方等的土地利用强度将区域土地利用类型人为活动划分为4个等级^[50]；依据人类活动程度, 利用各等级道路数据将其进行加权叠加计算道路密度^[51]；以DMSP/OLS夜间灯光数据表征人口的空间分布特征^[52]；综合考虑近年来耕地垦殖开发、湿地退缩及城镇化等现象导致土地利用结构失衡问题, 构建土地利用转移矩阵分析受人为与自然干扰引发的土地利用转移强度。

景观生态风险是指受人为或自然因素的作用, 景观格局和景观生态过程产生的累积性后果, 是生态风险评价在区域尺度上的重要分支^[30—31]。依据谢花林等的景观生态风险将其分为景观干扰度指数与景观脆弱度指数, 景观干扰度指数E_i是表示不同景观所代表的生态系统受外界干扰(主要是人类活动)的程度, 通过景观破碎度指数、景观分离度指数、景观分维度指数赋予权重叠加获得^[42]。

(3)

式中, C_i, N_i, F_i分别代表景观破碎度指数、景观分离度指数、景观分维度指数；a、b、c分别表示不同景观指数的权重。表 1为景观格局指数的生态学意义及计算公式, 结合三江平原实际情况及谢花林等的研究^[42], 将景观破碎度指数C_i、景观分离度指数N_i、景观分维度指数F_i的权重(a、b、c)分别赋值为0.5、0.3、0.2。

景观脆弱度指数V_i用以反映不同景观类型抵御外部干扰能力的大小。基于6种景观类型所代表的生态系统, 结合谢花林等研究以未利用土地最脆弱赋值为6, 其次水域为5、耕地为4、草地为3、林地为2, 城乡建设用地最稳定设置为1, 然后进行归一化处理, 得到景观脆弱度指数V_i^[42]。通过景观干扰度指数E_i和景观脆弱度指数V_i构建景观脆弱性评价模型, 考虑人为因素与自然因素的综合效应, 以此反映景观脆弱性问题。公式如下：

(4)

式中, EV_i代表景观脆弱性；E_i是景观干扰度指数；V_i是景观脆弱度指数。

本研究以三江平原2000年、2005年、2010年和2015年4期土地利用数据为基础, 从供给服务、调节服务、支持服务和文化服务4方面依据谢高地等的价值当量因子换算方法^[53—54], 结合三江平原土地利用类型的特征, 将土地利用类型与最接近的生态系统类型相对应(如表 2), 城乡建设用地包括城镇用地与农村居民点等, 取其ESV为0值。测算得出2000—2015年不同土地利用类型单位面积的生态系统服务价值系数, 计算生态系统服务价值公式如下：

(5)

式中, ESV为生态系统服务价值; C_ji是土地利用类型i的第j项生态系统价值系数；A_i是土地利用类型i的面积。表 2为生态系统服务价值土地利用类型对应表。

2.2 生态风险空间集聚特征分析

空间权重矩阵能够定量表达地理要素之间的空间关系, 是度量空间自相关的基础, 分为邻接权重矩阵和距离权重矩阵^[55]。邻接权重矩阵是指根据多边形的邻居关系确定, 包括以共边为邻接的车权重矩阵(Rook)和以共边或共点为邻接的皇后权重矩阵(Queen)。距离阈值定义了距离权重矩阵, 距离权重矩阵决定了空间自相关指数的大小。本文分别设置5、6、7、8、9、10、20、50km不同的距离阈值计算权重矩阵, 通过空间自相关指数的大小变化确定最优距离阈值。

空间自相关分析方法是表示事物或现象之间具有对空间位置的依赖关系, 包括全局空间自相关和局部空间自相关两方面, 常用Moran′s I指数表示, Moran′s I值介于正负1之间, 0值表示不存在空间相关性, 大于0表示存在空间正相关, 反之为空间负相关^[56]；生态风险与土地利用的空间分布、人类活动的动态发展等因素之间存在着一定的空间关联, 空间自相关模型在生态风险方面已有大量研究, 本文利用Geoda软件计算全局与局部空间自相关揭示三江平原生态风险的空间聚集现象。

全局空间自相关Moran′s I指数计算公式^[57]：

(6)

局部空间自相关计算公式^[57]：

(7)

式中, n为网格单元总数；x_i (x_j)为网格单元i(j)的风险值；(x_i-x) 为第i个网格单元上风险值与平均值的偏差；w_ij为标准化的空间权重矩阵。

可用标准化统计量z来表征空间自相关的显著性。计算公式为^[57]：

(8)

式中, E(I) 为I的期望值；VAR(I) 为I的方差, 当|z|＞1. 96, |P|＜0. 05时, 拒绝无效假设。

3 结果与讨论

图 3展示了2000年、2005年、2010年和2015年的风险源强度等级时空分布变化情况。从时间尺度看, 三江平原风险源强度呈增加趋势。其中, 2000年到2010年风险水平虽变化不大, 但2015年生态风险源强度等级愈加严重, 从网格水平均值统计得中高风险地区占总面积的70.92%, 较2000年增加了13.52%。从空间尺度看, 松花江、穆棱河、倭肯河一直处于中高风险水平, 且高等级道路网格分布与河流分布相对一致, 风险源等级与道路密度拟合度高, 随着道路密度的增加而增加, 高等级道路中心地所在网格风险等级高。DMSP/OLS夜间灯光数据较好的表征了行政区域多等级中心地的分布, 鹤岗市、佳木斯市、双鸭山市、七台河市与鸡西市5个地级市为高等级灯光值分布区, 其次县级所在地灯光等级较高, 灯光值与居住点的匹配有利于较好反映风险源的空间特征。土地利用变化与人类不合理的开发利用密不可分, 由于土地利用转移的频率和强度加大, 2015年河流两侧水田旱田面积明显增加, 2015年较2000年湿地面积减少了21.93%, 导致河流流域高风险突出, 出现以河流为中心向四周扩散趋势, 占用湿地过度围垦是风险源强度严重的重要原因。

图 4为2000年到2015年景观生态风险分布结果。由图 4可知, 景观生态风险分布随着时间变化总体上表现为增加趋势。2000年到2015年中高风险地区面积占比分别为26.94%、46.46%、61.86%、63.83%, 中高风险地区从网格水平上2015年较2000年增加了36.89%, 景观斑块破碎化程度高是导致研究区景观生态风险呈加剧趋势的重要原因, 并且较高风险和高风险等级面积持续增加。空间上, 景观生态中高风险地区重点集中在湿地与水体分布区, 景观自身稳定性被破坏, 受人为干扰后易损程度大。2000—2005年, 景观生态风险以低风险水平为主, 除低风险区和较低风险区面积有所减少外, 其他风险等级区均有不同程度的增加, 中风险区、较高风险区和高风险区面积分别增加了6.88%、3.37%、3.13%, 且集中分布在松花江、穆棱河、挠力河等区域, 此部分主要土地利用类型为水体和湿地, 受人为干扰易损程度大。2010—2015年, 景观生态风险以中高风险为主。高风险区面积有所增加, 其他生态风险等级区域面积均有减少。其中, 高风险区增加17.15%, 增加幅度大, 研究区内松花江、挠力河、黑龙江、乌苏里江及其周围湿地区高风险区突出, 受人类活动对湿地区的干扰程度加剧, 景观类型转移频率加大, 景观斑块破碎化程度高, 导致景观生态风险严重。

生态系统服务价值与生态风险呈负相关关系, 图 5为生态系统服务价值标准化后取负的生态风险等级结果。如图所示, 从时间上, 2000—2015年生态系统服务价值总体呈稳定趋势。虽波动幅度较小, 但生态风险以中高风险区域为主, 且重点集中在中部地区, 2000—2015年中高风险面积所占比例分别为67.54%、66.46%、66.52%、67.95%。土地利用变化是导致生态系统功能与结构发生变化的主要原因,, 且生态系统一旦遭到破坏, 短时间内难以快速恢复, 因此, 基于生态系统服务价值的生态风险在2000—2015年总体变化幅度较小, 且突出表现为中高风险。从空间格局上看, 生态系统服务低价值区即高生态风险区主要分布在中部水田、旱田、建设用地以及东北部与东南小范围的湿地区域, 由于城市扩张与耕地开垦消耗了大量的自然资源, 湿地、水体等生态系统服务高价值的地区被占用, 生态用地面积减少导致气候调节、水源涵养、土壤保持等生态系统服务功能受损, 生态系统服务低价值区域扩散分布, 中高风险等级区域出现集中化趋势。

图 6展示了2000—2015年三江平原综合生态风险的时空格局演变特征。时间上, 生态风险呈现明显的增加趋势, 2000年, 低风险、较低风险、中风险、较高风险及高风险等级网格数所占比例分别为39.68%、51.57%、4.81%、2.72%、1.23%。2015年, 这一比例分别为32.85%、15.53%、11.82%、14.50%、25.30%。较低风险区域下降了36.04%, 与此相反, 高风险等级网格数比例上升了24.07%。从空间看, 由于人为风险源干扰频繁, 2000—2015年湿地转移为其他生态系统类型面积共有3637.96km², 其中转化为农田面积达3235.16km², 湿地面积大幅度减少, 加之较高的景观生态脆弱性和生态服务价值的损失造成高风险上升区域主要集中在研究区中部的穆棱河、倭肯河、挠力河、松花江区域。随着水田与旱田面积增加, 高风险区由松花江河滩型湿地区与穆棱河的零星分布区逐渐向四周蔓延, 扩散区域主要为被侵占的湿地、林地区域。此外, 生态风险受城市扩张影响大, 高风险区域行政中心具有极高的拟合度, 鹤岗市、佳木斯市、鸡西市等城市显示为高和极高风险水平, 且辐射周边地区具有较高的风险水平, 以桦川县、勃利县、穆棱市、友谊县最为突出。

使用3种空间邻接矩阵, 计算土地利用与生态风险自相关的权重。如表 3所示, 在P值均小于0.01的显著性水平上, 使用二阶皇后矩阵(Queen2)所得的Moran′s I值小于车矩阵(Rook)的Moran′s I值小于一阶皇后矩阵(Queen1)的Moran′s I值, 说明在空间邻接矩阵中使用一阶皇后矩阵(Queen1)空间自相关效果好。基于距离的权重矩阵, 图 7表示分别选取5、6、7、8、9、10、20、50km为距离阈值构造距离权重矩阵计算所得的Moran′s I值和P值结果。结果发现在P值均小于0.01的显著性水平上, 当5km≤d≤7km时, Moran′s I值在2000年、2005年、2010年、2015年均为最高值, 分别是0.33、0.31、0.24、0.20。随着距离阈值的增大, 四年的Moran′s I值均呈下降趋势。由此说明, 生态风险的强弱受到空间距离的影响, 且距离阈值越大, 以某一网格为中心的相关区域越大, 整体的相关性指数降低。通过邻接权重矩阵与距离权重矩阵比较分析可知, 当5km≤d≤7km时, 所建立的矩阵能较好地表达生态风险的空间分布自相关情况。

本文对2000年、2005年、2010年、2015年三江平原的生态风险进行了空间自相关指数计算, 并绘制了空间集聚分布图。基于前文分析, 当距离阈值5km≤d≤7km时, 空间集聚性较高, 因此本研究选取5km为局部自相关分析的距离阈值, 在z检验的基础上(P=0.05)得到2000年、2005年、2010年、2015年三江平原综合生态风险空间集聚分布图(如图 8), 三江平原土地利用-综合生态风险分布的双变量空间自相关表示土地利用类型与其邻域生态风险均值间的相关性。由图 8可知, 2010—2015年两者的正相关型分布(高-高或低-低)较负相关型(低-高或高-低)分布比例较高, 高-高地区由7%增加到13%。2000—2015年空间自相关指数Moran′s I值呈现增加趋势, 且均通过了显著性检验, 其中2010年Moran′s I值最高为0.33, 这表示三江平原土地利用与综合生态风险的空间格局存在显著的空间正相关关系, 即随着土地利用程度的增强, 综合生态风险总体上呈现增强趋势, 这主要是受人类活动干扰引起土地利用的变化与转移, 导致综合生态风险空间上出现集聚现象, 综合生态风险分布范围有扩大趋势。

在空间演变上, 土地利用与综合生态风险呈高-高的地区主要集中在研究区内的松花江流域及周围滩地地区, 随着土地利用变化及转移, 空间关联逐渐增强且区域不断扩大, 到2010年, 两者成团状由松花江流域及周围滩地逐渐扩散到挠力河附近的水田及穆棱河中部地区, 到2015年, 区域扩大现象持续发生, 倭肯河区域的土地利用与综合生态风险的空间关联逐渐增强, 耕地面积的增加导致湿地水域区域被侵占, 生态服务价值降低, 景观脆弱性加剧, 造成综合生态风险呈现空间集聚与扩散分布现象。土地利用与生态风险呈低-低的地区主要分布在林地分布区, 这些地区由于地势高, 受人类干扰少, 土地利用变化小, 综合生态风险增加缓慢。不同时间尺度上, 负相关分布所占比例较少, 呈零星分布, 无明显的集聚特征。

4 结论与建议

本研究从风险源与风险受体两方面考虑了人类活动、土地利用变化、城市扩张、景观生态、生态服务价值等因素构建了生态风险评价体系, 并以三江平原湿地为研究对象, 评价了2000年、2005年、2010年、2015年的生态风险时空变化特征, 通过选取合适的距离阈值利用空间自相关模型对三江平原湿地的生态风险空间集聚效应进行了探究。基于以上研究, 本文得到的结论主要有：

(1) 从风险源角度, 高等级道路中心地所在网格风险等级高, 道路密度与风险源等级拟合度高, DMSP/OLS夜间灯光数据较好的表征了行政区域多等级中心地的分布, 土地利用变化突出表现为湿地面积减少, 耕地面积增加。三江平原湿地风险源强度呈增加趋势, 松花江、穆棱河、倭肯河地区一直处于中高风险水平；从风险受体角度, 2000年到2015年景观生态风险的中高风险地区重点集中在湿地与水体分布区；生态系统服务低价值区主要分布在中部水田、旱田、建设用地以及东北部与东南小范围的湿地区域, 中高风险等级区域出现集中化趋势。

(2) 综合风险源与风险受体的三江平原湿地生态风险时空变化可知, 从2000年到2015年湿地生态风险呈现明显的增加趋势, 高风险地区由松花江河滩型湿地区与穆棱河的零星分布区逐渐向四周蔓延, 受较高的景观生态脆弱性和生态服务价值损失的影响, 高风险上升区域集中分布在研究区中部的穆棱河、倭肯河、挠力河、松花江区域；受城市扩张影响, 鹤岗市、佳木斯市、鸡西市等城市显示高风险区域行政中心具有极高的拟合度, 且辐射周边地区具有较高的风险水平。

(3) 综合生态风险强弱受空间距离的影响, 本研究选取5km为局部自相关分析的距离阈值, 较好地表达综合生态风险的空间分布自相关情况。2000年、2005年、2010年、2015年三江平原土地利用与综合生态风险的空间格局存在显著的空间正相关关系。研究区内的松花江流域及周围滩地地区为综合生态风险突出热点区域, 随着土地利用变化及转移, 空间聚集性增强且呈现扩散分布趋势。

根据以上研究结果, 本文为三江平原湿地生态风险防控管理提出以下建议：提高重视湿地、水体等土地利用类型的侵占式利用与高强度转移是导致生态脆弱性显著的重要原因, 加强治理高风险区域集中分布的松花江流域挠力河等水体和湿地区。建议在城市路网规划与道路选址过程中加强分析道路建设与生态系统之间的交互作用机制, 通过自上而下的管理方式降低人类活动对生态环境的不良影响；严格把控城乡建设与耕地扩张的生态红线, 合理规划城市人口与道路建设, 以降低局部生态风险；优化土地利用模式, 增强区域各景观间的连通性, 恢复已破坏的景观, 保护湿地水域等景观斑块的完整性, 避免风险转移；明确三江平原湿地生态系统的生产能力与经济价值, 合理增加湿地与林地等高生态服务价值的面积, 增强系统自身稳定性, 提高三江平原湿地生态风险防控化解能力。