松花江流域哈尔滨段景观生态风险评价

张玉娟，曲建光，叶猛猛

(黑龙江工程学院测绘工程学院，黑龙江 哈尔滨 150050)

0 引言

流域生态环境变化及其引发的各种负面影响越来越受到人们关注[1-2]，改变土地的利用方式是人类开发利用自然环境最主要的表现形式，不合理的土地利用方式会威胁区域的生态环境健康[3-4]. 生态风险评价是一种能有效支持生态系统管理的工具，通过生态风险评价能够评估环境污染、 人为胁迫或者自然灾害等外界干扰对区域生态系统及其结构产生负面效应的可能性及危害程度[5-7]. 这有助于找出区域发展过程中存在的环境问题，对区域土地科学利用与合理开发、 人类经济活动和区域可持续发展具有重要意义[8-9]. 生态风险评价出现由传统生态风险评价, 到区域生态风险评价, 再到景观生态风险评价的发展趋向. 景观生态风险评价通过选择与生态风险相关的景观指数结合权重构建评价模型，如吕乐婷等[5]采用景观干扰度和景观脆弱度结合景观类型面积比例构建景观生态风险模型对1985—2015年细河流域景观生态风险时空分布特征及空间关联格局进行评价；汪翡翠等[10]采用斑块密度和蔓延度结合层次分析法构建景观生态风险模型对京津冀城市群土地利用生态风险空间分布特征进行分析；胡绵好等[11]通过对景观类型设置风险等级值结合景观类型面积比例对德兴市生态风险进行了分析. 景观生态风险评价注重风险的时空异质性和尺度的效应，致力于实现风险的空间表征及其可视化，目前得到了广泛应用.

松花江流域哈尔滨段位于松花江流域中部，地属东北平原，地域平原分布辽阔、 土壤肥沃，是重要的粮食生产区, 在行政区划上属于 “一带一路”战略节点区域—哈尔滨市. 因此，这一流域的景观格局不仅对本区域生态安全产生影响，对整个松花江流域生态格局乃至东北对外开放水平都产生重要的作用. 近20年来，该区域土地利用方式不断发生变化[12]，生态安全性也随之发生变化，目前针对该区域粮食安全和水环境评价方面有较多的研究[13-14]，但景观生态风险未见研究. 鉴于此，本文以松花江流域哈尔滨段为研究对象，利用ENVI5.3、 GS+、 ArcGIS10.3、 Geoda095i等软件提取研究区域1995年、 2005年和2015年三期土地利用类型，通过景观易损度指数及景观类型面积比重结合层次分析权重赋值法构建并计算景观生态风险指数，借助冷热点、 空间变差函数法，对松花江流域哈尔滨段景观生态风险的时空变化特征进行研究. 松花江流域哈尔滨段景观生态风险状态研究成果是景观格局优化的基础，能够明确景观格局规划的保护重点，有助于完善松花江流域哈尔滨段的生态格局，为区域生态环境的进一步治理以及土地的合理利用提供科学依据.

1 数据准备及研究方法

1.1 数据来源及处理

本文遥感数据均来源于地理空间数据云，年份分别为1995年、 2005年和2015年，行带号分别为116/28、 117/28、 118/28、 117/29、 118/29. 非遥感数据主要包括哈尔滨市行政区划图，哈尔滨市部分土地利用现状数据，外业调查数据，哈尔滨市1996、 2006和2016年鉴统计中各区县人口密度数据等. 根据2017年国家标准《土地利用现状分类》，结合本研究区域的主要土地利用方式，将研究区划分为耕地、 林地、 建设用地、 水域、 草地以及未利用地6种土地利用与覆盖类型. 利用ENVI 5.3遥感影像处理软件，对下载获得的初始影像分别进行辐射定标、 大气校正、 裁剪、 镶嵌等预处理，根据外业调查数据、 土地利用现状数据及Google earth地图等相关资料，利用最大似然法进行遥感影像监督分类，获得松花江流域哈尔滨段1995年、 2005年和2015年土地利用与覆盖类型分布图(见图1)，三期影像分类精度分别为0.871、 0.862和0.888.

1.2 研究方法

1.2.1 景观生态风险指数

为建立土地利用类型与区域生态风险之间的联系，利用景观易损度结合各景观类型面积比重及层次分析权重赋值法构建景观生态风险指数(ecological risk index, ERI)，该值的大小反映了区域景观生态风险的高低，其计算模型如下：

(1)

式中：IER为生态风险指数；n为景观类型的数量；Ai为研究区内景观类型i的面积(km2)；A为研究区域总面积(km2)；Wi为研究区内景观类型i所反映的生态风险强度权重；Ri为景观类型i的易损度，选择与景观类型易损程度相关的景观破碎度、 分维数和分离度组合构建，按下式进行计算.

Ri=aCi+bNi+cFi

(2)

式中：Ci为景观类型i的破碎度；Ni为景观类型i的分维数；Fi景观类型i的分离度[15]，为了消除不同量纲的影响，对计算所得的各景观指数进行归一化处理；a、b、c分别为Ci、Ni和Fi的权重.

评价指标权重赋值是否合理将直接影响景观生态风险指数的合理性，本文采用AHP[16-17]的原理和方法确定各景观类型生态风险强度权重和景观破碎度、 分维数、 分离度指数权重. 首先查阅现有的景观类型生态风险强度参数相关文献中各种景观类型权重值[18-20]，将这些权重按平均值大小进行排序，作为初始专家经验构造判断矩阵，通过对该判断矩阵进行一致性检验来判断所构造的矩阵是否具有满意的一致性. 本文通过计算获得一致性比例的计算结果为0.007 8(<0.1)，表明判断矩阵具有满意的一致性，获得6种土地利用类型生态风险的AHP权重分别为：建设用地0.366、 林地0.049、 水域0.107、 草地0.075、 耕地0.146、 未利用地0.257. 同理，获得景观破碎度、 分维数、 分离度指数权重值分别为0.539、 0.297、 0.164，一致性比例的计算结果为0.003 5(<0.1).

1.2.2 景观生态风险空间统计方法

1) 空间变异函数. 在地统计学上，通过半方差函数的理论方法分析地理变量空间上的关系及分布格局[21]，实现判断地理变量在单元网格尺度划分下是否具有空间相关性，从而获得地理变量在多尺度网格划分中，哪个尺度更为合理. 半方差函数的理论模型计算公式为：

(3)

式中：γ(h)为变异函数值；Z(xi)和Z(xi+h) 是在空间单元xi和xi+h上的景观生态风险值(i= 1, 2, 3, …,N(h))；N(h)是分割距离h的样本量.

2) 冷热点分析. 莫兰指数被用于判定地理变量空间自相关关系，包括全局自相关莫兰指数和局部自相关莫兰指数. 其中全局自相关莫兰指数用于分析地理变量空间相似性聚集度，局部自相关莫兰指数用来识别地理变量在不同空间单元的高值簇或热点区与低值簇或冷点区的空间分布[23]. 全局自相关莫兰指数和局部自相关莫兰指数计算公式为：

(4)

(5)

式中：Pe和Pf分别为地理单元e和地理单元f位置的函数值(e≠f)；Wef为权重系数矩阵；n为地理单元数. 莫兰指数的判断区间为[-1, 1]，区间(0, 1]表示正自相关，值越大，其正自相关性越高；区间[-1, 0)表示负自相关，值越小，其负自相关性越高；值等于0时表明变量值在空间上是相互独立分布的.

2 实证研究

2.1 研究区域概况

松花江流域哈尔滨段总面积5.31万km2，区域位于松嫩地块-松嫩断陷东南隆起区，东南临张广才岭支脉丘陵，北部为小兴安岭山区，中部有松花江通过，山势不高，河流纵横，平原辽阔. 区域属北温带季风气候区，大陆性气候特点明显，多年平均降水量一般在500 mm左右，汛期6-9月的降水量占全年的60%～80%，冬季12月至次年2月的降水量仅为全年的5%左右. 研究区域内的大小河流均属于松花江水系和牡丹江水系，松花江干流由西向东贯穿研究区域中部，是研究区域灌溉量最大的河道. 研究区域水资源特点是自产水偏少，过境水较丰富，时空分布不均匀，表现为东富西贫.

2.2 景观生态风险小区划分及生态风险指数计算

景观生态风险指数由景观格局指数构建，因此具有尺度效应，本文分别按照边长4和6 km正方形单元将研究区域分割成2 300和1 065个研究小单元，利用ArcGIS 10.3中的渔网工具和分割工具对三期土地利用与覆盖图进行单元分割，对获得矢量小单元进行批量转换生成为栅格小单元. 利用景观指数提取软件Fragstats 4.2批量提取三期每个网格4个类型水平的景观指数：斑块面积(TA)、 斑块数(NP)、 分维数(FRAC_MN)和分离度指数(DIVISION). 在Excel中按照式(1)编辑函数，计算三期每个小单元的景观生态风险值，并将该值作为网格中心点的属性值.

2.3 空间变异函数理论模型拟合分析

不同步长下的景观生态风险变异函数模型能准确反映景观生态风险在不同尺度上的变化特征. 确定最佳的采样尺度是进行景观生态风险空间分异及关联研究的基础与前提. 本文采用的是规则格网取样，可将格网间距4和6 km作为步长大小来分析在这两种采样尺度下的空间变异性，从而确定能在一定空间范围内反映景观生态风险的最佳尺度. 利用地统计学GS+软件进行两种尺度下景观生态风险的统计分析，构造景观生态风险空间变异函数理论模型，完成景观生态风险数据变异函数模型拟合(见表1和表2).

表1 4 km尺度变异函数拟合模型参数

注：Rc为结构方差占总方差的比值C/(C+C0)，下同.

表2 6 km尺度变异函数拟合模型参数

续表2

4 km尺度下，三期景观生态风险均为线性模型拟合效果最好，三期决定系数R2分别为0.974、 0.974和0.969，此时残差RSS均最小. 结构方差占总方差的比值Rc呈现逐渐增加趋势，表明在4 km尺度划分下，1995—2015年，研究区域景观生态风险受结构性因子的影响有所增加[24-25]. 研究区三期景观生态风险空间分异的有效变程A为18.31 km(大于4 km)，因此研究区域景观生态风险在4 km尺度划分下有较高度的空间相关性.

6 km尺度下，三期景观生态风险均为指数模型拟合效果最好，三期决定系数R2分别为0.985、 0.981和0.990，此时残差RSS均最小. 结构方差占总方差的比值Rc呈现逐渐增加趋势，表明在6 km尺度划分下，1995—2015年，研究区域景观生态风险受结构性因子的影响也是有所增加的. 研究区三期景观生态风险空间分异的有效变程A为37.08 km(大于6 km)，因此研究区域景观生态风险在6 km尺度划分下具有非常高度的空间相关性.

4 km和6 km两种尺度划分都适合于研究区域景观生态风险研究，二者相比较，6 km尺度下决定系数R2值更大、 残差RSS值更小、 结构方差占总方差的比值Rc更大些，景观生态风险空间相关性更高些，因此，选择6 km网格划分尺度进行研究区景观生态风险评价. 6 km尺度下，区域景观生态风险网格分布图(见图2).

2.4 生态风险时空变化分析

根据自然断点法将景观生态风险指数划分为5个景观生态风险等级，即低风险级、 较低风险级、 中等风险级、 较高风险级、 高风险级. 采用ArcGIS 10.3地统计模块的克里格插值法对划分等级后的景观生态风险离散点进行插值，形成三期景观生态风险等级连续分布图(见图3)，分别对三期景观生态风险等级进行统计(见表3).

表3 三期景观生态风险分级统计

近20年间，各风险等级分布面积不断变化，但分布的主要区域位置无明显变化. 高风险区主要分布在主城区，较高风险区主要分布在高风险区周围以及各县域城镇中心，原因为主城区分布较大面积的建设用地景观和极小面积的林地景观，县域城镇中心主要以建设用地景观和耕地景观为主. 低风险区主要分布在通河县、 木兰县、 方正县和宾县境内，原因为这些县域分布着很大面积的林地景观. 中等风险区面积最大，占研究区域41%～47.6%，主要分布于大面积的耕地景观范围内. 较低风险区主要分布在中等风险区与低风险区的过渡区域范围. 1995—2015年，高风险区和较高风险区面积呈现先增加后减少的变化，高风险区分布面积比例由6.60%到11.00%再到7.90%，较高风险区，分布面积比例由13.80%到19.90%再到19.30%. 低风险区、 较低风险区和中等风险区面积呈现先减少后增加的变化，低风险区分布面积比例由14.50%到14.00%再到14.50%，较低风险区分布面积比例由17.50%到14.10%再到16.10%，中等风险区分布面积比例由47.60%到41.00%再到42.20%. 从整体水平上看，研究区景观生态风险呈现先上升后下降的变化.

2.5 区域生态风险等级的空间关联格局分析

表4 全局莫兰指数及验证值

利用ArcGIS 10.3空间统计工具分析研究区域景观生态风险空间关联性. 利用全局莫兰指数计算工具计算景观生态风险全局莫兰指数值和Z值(见表4). 三期全局莫兰指数分别为0.512 5、 0.555 2和0.605 5，均大于0；Z值分别为24.25、 25.69和27.44，说明研究区域景观生态风险等级存在很强的正相关性，呈聚类模式的空间分布.

通过对景观生态风险等级分布进行热点分析，获得相邻区域之间景观生态风险的空间关联模式. 利用ArcGIS 10.3中聚类和异常值分析工具进行相似景观生态风险的空间聚类分析，结果见图4. 从图4分析可知，1995—2015年，研究区域景观生态风险高-高聚集、 低-低聚集分布区域并未发生较大变化，但逐渐表现为较分散的分布. 高-低、 低-高聚集出现的非常少，只在1995年和2015年聚集分布于研究区域内.

3 结语

研究区域景观生态风险具有高度的空间相关性，结构性因素是三期景观生态风险空间相关性的主要影响因素，且结构性因素对景观生态风险的影响程度逐渐加深.

从整体上看，1995—2015年松花江流域哈尔滨段景观生态风险先升高后下降. 源于1995—2005年间，区域忽视生态环境保护，转林为耕、 开荒建城，使得生态用地不断转换为非生态用地，导致区域景观生态风险增加. 2005—2015年间，随着全球生态环境保护意识增强，区域制定了多项自然保护区发展战略，逐步恢复松花江沿岸生态环境，退耕还林、 退耕还草，使得部分非生态用地逐渐转换为生态用地，区域景观生态风险降低.

1995—2015年，研究区域生态风险高-高聚集、 低-低聚集分布区域并未发生较大变化，但逐渐表现为较分散的分布. 高-低、 低-高聚集出现的非常少，只在1995年和2015年聚集分布于研究区域内.

本文通过景观易损度指数及景观类型面积比重结合层次分析权重赋值法构建景观生态风险指数，具有一定的创新性. 生态系统服务是链接人类福祉和生态系统的桥梁，将景观生态服务价值与景观格局相结合进行生态风险研究将是课题进一步研究内容.