区域生态风险评价是指从区域的尺度上, 评价自然及人类活动对区域生态系统所造成的不利影响的可能性和危害程度的评价^[1]。土地利用/覆被变化(LUCC)是人类开发利用自然环境最重要的表现形式, 其方式的变化不仅对土壤、大气、水等综合地理要素产生影响, 同时会威胁区域的生态环境健康^[2-5]。流域作为人类活动的重要场所, 是连接地表水与地下水、土地覆盖和生态系统的重要综合生态地理区域^[6]。流域生态风险评价是流域生态环境保护与管理的重要研究内容, 基于流域土地利用/覆被变化及景观生态变化的流域生态风险评价工作是当前流域生态管理工作关注的重点内容之一。

与一般的区域生态风险评价相比, 流域生态风险评价具有其独特的流域特征^[7]。流域生态风险评价方法主要分为基于单一或多源污染因子评价法和景观分析法^[8]。景观分析法是选取适当的景观指数, 以景观类型及格局作为评价受体, 对生态风险的分布和演化过程进行定量分析^[7], 目前在国内外已得到广泛应用。如：Potter等^[9]探索了美国北卡罗来纳州73个流域的景观格局与非点源污染生态风险之间的关系, Obery等^[10]评价了美国Codorus Creek流域景观格局对脆弱生境的累计风险效应；黄木易等^[11]探索了近20年来巢湖流域景观格局变化特征及生态风险时空演化规律；刘永超等^[12]基于景观干扰度指数和景观脆弱度指数评估了近30年美国Tampa Bay流域景观格局及生态风险的时空特征；刘世梁等^[13]基于景观格局和土壤侵蚀过程评估了云南省红河流域景观生态风险分布规律。这些研究成果为国内外流域生态规划、景观结构调整及格局优化、社会经济可持续发展决策等提供了理论依据。目前, 流域生态风险评价研究已涵盖了湿润区、半湿润区以及干旱区流域的湖泊、河流、河口三角洲、海岸带等地^[14], 而以开发历史悠久的东北老工业基地采矿流域作为研究对象的案例尚属空白。

细河为辽宁省太子河一级支流, 流域煤炭开采历史悠久。长期粗放式开采使矿区地下水水位大幅下降, 采矿区土地资源和森林资源遭到严重破坏, 导致流域景观生态发生变化。本研究结合细河流域的实际情况, 基于流域1985、1995、2005和2015年遥感影像解译而得的土地利用图, 采用由景观干扰度及景观类型脆弱度指数共同构建的景观生态风险指数, 对流域生态风险的时空分布特征及空间关联特征进行了评估。研究结果可为基于生态安全的流域土地利用规划与管理提供科学依据。

1 研究区概况与数据处理 1.1 研究区概况

细河流域位于辽宁省本溪市境内, 地处123°30′—124°00′E, 40°45′—41°20′N。全长120 km, 流域总面积1119 km²。细河主要流经本溪市南芬区、平山区和本溪满族自治县。流域山地较多, 地势东南高、西北低(图 1)。细河流域属于温带大陆性季风气候区, 年均降水量800 mm, 主要集中于夏季6—9月份^[15]。流域内有丰富的铁矿资源和森林资源。

1.2 数据来源与处理

本文所用的1985、1995和2005年3期Landsat TM及2015年Landsat OLI遥感影像来源于地理空间数据云(http://www.gscloud.cn), 空间分辨率为30 m, 影像轨道号/行号分别为119/31和119/32, 每期2景, 成像时间均在植被生长茂盛期6—10月, 影像的云量均低于2%。借助ENVI 5.1软件对图像进行辐射定标、大气校正、影像的拼接与裁剪, 并根据实际情况选择7、5、2波段组合方案, 同时对影像进行2%的拉伸以增强解译效果。本文利用最大似然法, 参考1984年我国农业区划办组织制定的《土地利用现状调查技术规程》将流域的土地利用分为6个一级分类, 详见流域土地利用分类图(图 2), 并结合Google Earth建立混淆矩阵对结果进行精度验证。验证结果显示, 所解译的4个时段土地利用图Kappa系数均在0.8以上, 达到了中分辨率遥感影像精度使用要求^[16-17]。本文所用DEM数据来源于地理空间数据云的ASTER GDEM(V1)数据集, 栅格大小为30 m。并以DEM数据为基础, 借助ArcGIS中的水文分析模块, 生成细河流域集水区域。

2 研究方法 2.1 流域生态风险分析 2.2.1 生态风险评价单元的划分

为了能够将景观生态风险指数空间化, 在考虑景观空间异质性、斑块大小和流域面积的基础上, 本研究按照细河流域景观斑块平均面积的2—5倍进行网格采样^[18]。采用等间距的方式将细河流域划分为340个2 km×2 km的正方形网格生态风险评价单元(图 1)。

2.2.2 景观生态风险指数构建

本研究依据土地利用类型的面积比重和景观损失度指数R_i构建景观生态风险指数ERI^[19-20], 公式如下：

(1)

式中, ERI_k为景观生态风险评价单元k的景观生态风险指数, 该值越大表示该评价单元的生态风险程度越高, 反之, 生态风险程度越低。A_ki为景观生态风险评价单元k中i类景观的面积, A_k为景观生态风险评价单元k的总面积, R_i为i类景观的损失度指数, 通过景观干扰度指数E_i和各景观类型的脆弱度指数V_i构建^[21], 公式如下：

(2)

其中景观干扰度指数E_i^[21]公式如下：

(3)

式中, C、N、F分别为景观破碎度指数、分离度指数和分维数指数, 其生态学含义及计算方法见表 1；a、b、c分别为各景观指数的权重, a+b+c=1, 根据相关研究^[22]并结合细河流域实际情况, 分别赋值0.5、0.3、0.2。

景观类型的脆弱度指数V_i由专家咨询法并归一化获得, 具生态学含义见(表 1)。

2.2 空间自相关分析

空间自相关分析通过描述某一要素的属性值与其在空间上相邻的各要素属性值之间是否存在显著关联的关系, 来揭示空间参考单元和相邻单元在属性特征值方面的空间相关特征。其指标分为两种：全局Moran′s I指数用于检验某一要素的属性值在整个研究区内的空间相关性；局部Moran′s I指数用于反映某一要素的属性值和相邻空间单元的相关性^{[13, 28]}。全局空间自相关指数的计算公式如下^[29]：

(4)

(5)

(6)

式4中x_i表示第i地区的观测值, n为栅格数, W_ij是二进制的邻接空间权重矩阵, 用来表示空间对象的邻接关系。i=1, 2, …, n；j=1, 2, …, m；当区域i和区域j相邻时, W_ij=1；当区域i和区域j不相临时, W_ij=0。

局部空间自相关是将Moran′s I指数分解到各个独立的空间要素公式为：

(7)

(8)

Moran′s I值介于-1—1, 当Moran′s I＞0时, 表明研究区存在正相关, 研究单元的属性值呈趋同集聚；当Moran′s I＜0时, 表示负相关, 呈离散分布；当Moran′s I＝0时, 表示不存在空间相关性。

3 结果与分析 3.1 土地利用变化分析 3.1.1 土地利用结构变化分析

细河流域1985—2015年土地利用变化状况如图 2和表 2所示。在1985—2015年期间, 细河流域各用地类型中林地、耕地和水域的面积减少, 工矿及城镇建设用地、未利用地和草地的面积增加。其中林地的面积减少最多, 减少量为65.9 km², 面积比例从78.47%减少至72.58%；工矿及城镇建设用地的面积增加最多, 增加量为103.87 km², 面积比例从1.96%增加至11.24%；未利用地和草地的面积比例分别由0.64%和0.25%增加至1.58%和1.87%；耕地和水域的面积比例分别由17.38%和1.30%减少至12.13%和0.59%。

3.1.2 土地利用转化过程分析

细河流域1985—2015年土地利用转移矩阵如表 3所示。可见, 1985—2015年, 林地是细河流域土地利用的主要转出类型, 其面积主要转出为工矿及城镇建设用地、耕地和未利用地；工矿及城镇建设用地是细河流域土地利用的主要转入类型, 其面积的增加主要来源于耕地、林地和水域；耕地部分转出为建设用地, 部分由林地转入, 转出面积大于转入面积, 造成其总体面积减少。其他土地利用方式也产生了不同程度的转化。

3.2 流域景观生态风险时空变化分析

为了探索细河流域景观生态风险的空间分布特征, 将每个生态风险评价单元的ERI作为其中心点的属性值^[30], 利用ArcGIS地统计模块中的普通克里金插值得到细河流域景观生态风险的空间分布。并结合细河流域的实际情况和4个时相的ERI指数, 利用自然断点法将细河流域的ERI划分为5个生态风险等级：高生态风险区Ⅴ(ERI＞0.145)、较高生态风险区Ⅳ(0.115＜ERI≤0.145)、中生态风险区Ⅲ(0.093＜ERI≤0.115)、较低生态风险区Ⅱ(0.071＜ERI≤0.093)、低生态风险区Ⅰ(ERI≤0.071)。结果如图 3和表 4所示。

细河流域景观生态风险整体呈上升趋势, 其中较低、中、较高景观生态风险区所占面积较大, 而低、高生态风险区所占面积较小。从各级景观生态风险区面积变化来看, 1985年细河流域低、较低景观生态风险区面积比例分别为18.3%和33.1%, 占细河流域总面积的51.4%, 流域景观生态风险整体呈现较低的趋势。其中低景观生态风险区主要分布于流域北部海拔较低的本溪市平山区, 此区域的基质景观类型以城镇建设用地为主, 受到人类干扰的后的损失度低；高景观生态风险区主要分布在流域中东部地区, 此区域经济欠发达以农耕为主, 耕地和未利用地面积较广且有零散破碎的工矿用地, 开凿矿山, 破坏了山地景观自身的稳定性, 景观受人为干扰后易损程度大, 耕地景观的斑块较破碎。

1985—2005年, 细河流域低、较低、中景观生态风险区面积减少而较高、高景观生态风险区面积增加。其中较高、高景观生态风险区的空间分布格局也发生了变化, 自流域中东部地区向南部本溪满族自治县扩散。在研究时段内本溪满自治县为了发展农业, 将大面积林地开垦为耕地, 使得原有较完整的景观变得破碎, 景观生态风险增加；低景观生态风险区面积逐渐减少并在流域北部地区集聚, 在研究时段内流域北部的本溪市平山区随着经济的不断发展, 城镇建设用地面积增加且从无序状态向有序方向演变, 系统稳定性逐步提高。

2005—2015年, 细河流域中等、较高、高景观生态风险区的面积比例分别为33.9%、20.9%和7.7%, 占流域总面积的62.5%, 较1985年相比呈明显的上升趋势。较高、高景观生态风险区较2005年有所下降, 流域南部本溪满族自治县的景观生态风险有所改善。在这10年间, 本溪满族自治县大力加强城市建设, 县内的城镇建设用地面积增加且向有序方向转变, 受人为活动干扰后易损程度低, 且流域中东部地区开始施行重建矿区生态环境等措施, 降低了矿区生态风险程度。

3.3 景观生态风险空间相关性分析 3.3.1 全局自相关分析

依据1985—2015年细河流域景观生态风险空间分布数据得到Moran′s I散点图(图 4)。细河流域的景观生态风险在1985年、1995年、2005年和2015年的全局Moran′s I值分别为0.3254、0.3656、0.3529和0.3665, 呈波动上升的趋势。这说明细河流域的景观生态风险在空间上呈现集聚效应；景观生态风险值高的区域, 周边的景观生态风险值亦高；景观生态风险值低的区域, 周边的景观生态风险值亦低。

3.3.2 局部自相关分析

进一步对细河流域景观生态风险局部空间相关性进行分析, 得到局部空间自相关集聚图(图 5)。可见, 1985—1995年, “高-高”值区皆集中分布在流域中东部地区, 此区域经济欠发达, 用地模式以“开垦-废弃-再开垦”的农耕经济为主, 耕地和未利用地面积较广。另外, 零散破碎的工矿用地及矿山开凿, 破坏了山地景观自身的稳定性, 使景观损失度增加；“低-低”值区呈零散分布, 主要分布于本溪市平山区内, 这与同期的景观生态风险克里金插值空间分布格局较为一致。1995—2015年细河流域景观生态风险的“高-高”值集群结构发生了变化, 向流域南部扩散, 而“低-低”值区变化不大。位于“高-高”值区内的本溪满族自治县, 地处低山丘陵区, 基质景观类型以林地为主, 但由于追求发展, 开荒毁林较为严重。城镇建设用地和耕地都呈现较小且分散的分布模式, 景观连通性、系统内部的稳定性差, 景观损失度高。

4 结论与讨论

本文从土地利用变化及转移、景观生态风险的时空分布及空间相关3个方面对细河流域的景观生态风险进行了综合评价, 结果表明：

1) 1985—2015年, 细河流域6种用地类型的面积皆发生了变化, 其中草地、工矿及城镇建设用地和未利用地面积增加, 林地、耕地和水域面积减少。林地是流域的优势景观类型, 对流域景观格局的变化起着重要作用。

2) 1985—2015年, 细河流域高、较高和中景观生态风险区的面积比例分别增加了2.4%、8.7%和2.8%, 且在空间上由流域东中部地区向南部转移；低景观生态风险区和较低景观生态风险区面积比例分别减少了8%和5.9%, 且在空间上向流域北部地区集聚。流域景观生态风险整体呈现增高趋势。

3) 在研究时段内, 流域的全局自相关Moran′s I值分别为0.3254、0.3656、0.3529和0.3665, 景观生态风险呈正相关现象, 在空间分布上趋于集群。且景观生态风险局部自相关格局较为一致, 高-高值分布在流域东部地区, 而低-低值分散于流域周围。

本研究以景观格局指数为指标, 通过土地利用类型的面积比重、景观干扰度和景观损失度指数构建景观生态风险指数, 仅从景观空间结构的一个角度出发来评价细河流域的生态安全格局, 没有考虑流域气候、地形、地貌和社会经济等因素, 所以并不具有绝对性。但流域景观格局是流域生态环境管理的基础单元、是人类活动的主要作用面, 其变化必然会引起生态系统功能的改变。因此, 利用景观格局指数的方法研究流域的生态安全格局是可行有效的。针对本研究, 建议细河流域结合自身实际情况, 对现有土地利用类型的组合情况进行思考, 在结合流域自然因素的前提下, 改善不合理的土地利用类型。在资源开发的过程中应坚持生态优先的原则, 在保护河流水域和林业发展的同时, 限制对地表扰动较大的工矿建设用的大规模发展, 保留一定范围的生态缓冲区, 并实施有效的生态补偿制度。受研究者知识、经验和主观因素的限制, 对于细河流域景观生态风险时空关联格局的驱动力解释还不够完善, 对此进行深入探索并进而提出可行有效的流域生态管理建议是本研究今后努力的方向。