湿地作为地球上生产力最高的生态系统, 拥有丰富的生物多样性和宝贵的生态服务功能, 其在调节气候、涵养水源、防洪蓄洪、土壤保持等方面的作用是其他生态系统所无法取代的^[1-2]。乌裕尔河-双阳河流域(以下简称：乌双流域)作为我国东北沼泽湿地的重要分布区, 其在区域国土生态安全中的地位十分重要。然而近几十年来, 由于气候变化和人类活动的影响, 湿地景观格局和水文情势已经发生了深刻变化, 面临着沼泽湿地面积日渐减小、湿地破碎化严重、湿地功能急剧下降、盐碱化程度加剧的局面^[3-5]。

遥感技术的快速发展实现了湿地景观大范围、多时相数据的获取, 使得湿地生态系统的宏观动态监测及其驱动力分析成为当前湿地研究的热点之一^[6-7]。驱动力研究能够帮助揭示湿地景观格局演变的内部机理, 并由此可以预测未来的演变趋势及制定对应的管理政策^[8]。目前大多数研究主要通过分析气候、水文等自然因子^[9-10]和人口、GDP等社会经济因子^[11-12]的变化趋势来探究与湿地变化间的关系, 有效弥补了以往定性分析只能罗列和描述因子的缺陷。此外一些学者也尝试对湿地演变驱动因子进行综合性分析, 找出影响湿地变化的主导因素, 如Gong等^[13]利用logistics回归分析得出降水量和入水量是影响北京市水库湿地演变的主导因素；吕金霞等^[14]运用GIS空间分析和主成分分析方法探讨京津冀地区湿地景观时空变化及其驱动力, 发现城市扩张和农业发展是湿地减少的主导因素；刘吉平等^[15]运用景观指数和地理探测器模型定量分析三江平原沼泽湿地变化的区域分异及影响因素, 结果表明人为干扰和地形地貌主导了湿地的变化。

以往驱动力研究的重心在于明确影响湿地变化的主导因子, 对因子贡献程度量化及影响机制研究较少, 而如何有效遏制湿地损失以及制定对应的湿地管理与保护方案需要建立在对湿地变化驱动机制深刻理解的基础上。鉴于上述, 本文选取乌双流域作为研究区, 提取1980—2015年7期湿地分布数据, 综合运用空间分析、景观指数、转移概率矩阵方法揭示湿地景观格局的时空变化特征, 并利用增强回归树模型分析湿地变化的驱动机制, 内容包含驱动因子贡献率、主导因子影响机制和未来湿地变化预警。

1 研究区概况

乌双流域位于黑龙江省西部, 东经123°47′—127°27′, 北纬46°23′—48°25′, 范围涉及齐齐哈尔市、大庆市、黑河市和绥化市共20个区县, 流域总面积约2.49万km², 地形特征为东北高、西南低(图 1)。乌裕尔河全长587 km, 是黑龙江省最大的内陆河, 是嫩江左岸较大的无尾河流, 其上游地形主要为低山丘陵区, 有明显的河床, 是主要的产流区, 下游地势低洼, 排水不畅, 河水漫溢失去河道而形成大范围的沼泽湿地——扎龙湿地。双阳河全长89 km, 自南向北流过一段后折向西南, 在尾闾低洼沼泽地消失, 与乌裕尔河共同形成闭流区。乌双流域属温带大陆性季风区, 冬季干燥严寒, 夏季温热多雨, 年降水量在400—500 mm左右, 主要集中在汛期6—9月份, 占全年降水量的80%, 历时短, 强度大, 径流的主要补给方式是降水。

2 数据与方法 2.1 数据来源

本研究使用的数据包括1980年、1990年、1995年、2000年、2005年、2010年、2015年7期土地利用数据, 气温和降水数据, 数字地面高程(Digital Elevation Model, DEM)数据, 道路数据、自然保护区边界数据以及DMSP/OLS夜间灯光数据。

土地利用数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/), 空间分辨率为1 km, 该数据由Landsat影像解译判读得到, 影像时相以秋季为主, 经野外验证其分类精度达91.2%^[16], 满足本研究的数据精度要求。根据数据源的土地利用分类系统和乌双流域的实际情况, 主要统计水田、旱地、林地、草地、居民地、水域、沼泽地、盐碱地八大土地利用类型, 将沼泽地、水域纳入乌双流域的湿地范围。气象数据来源于中国气象数据网(http://data.cma.cn/), 年均温和年降水量数据由流域所处东北地区共计102个气象站点月平均气温和月降水量观测数据经过Kriging插值、裁剪及栅格计算得到。DEM数据来源于地理空间数据云平台(http://www.gscloud.cn/)。道路矢量数据来源于OpenStreetMap平台(https://www.openstreetmap.org/)。由于人口经济因子的空间分布数据难以获取, 本研究用多年平均的DMSP/OLS夜间灯光数据来反映流域内的人口分布及经济发展情况, 数据由美国国家地球物理数据中心提供(https://www.ngdc.noaa.gov/)。

2.2 研究方法 2.2.1 景观指数

景观指数能够从看似无序的景观要素镶嵌中发现潜在的有意义的规律, 即用定量化的指标来反映景观结构组成和空间配置等方面的特征^[17]。本文依据生态意义显著、对湿地格局变化敏感且易于解译、彼此相互独立的原则选取以下景观格局指数：在类型水平上, 包括斑块数量(Number of Patches, NP)、平均斑块面积(Mean Patch Size, MPS)、聚集度指数(Aggregation Index, AI)、面积加权的平均形状因子(Area-Weighted Mean Shape Index, AWMSI)4个指数；在景观水平上, 包括最大斑块指数(Largest Patch Index, LPI)、斑块密度(Patch Diversity, PD)、聚集度指数(Aggregation Index, AI)、斑块结合度指数(Patch Cohesion Index, COHESION)4个指数。上述景观格局指数均利用Fragstats 4.2软件计算得到。

2.2.2 转移概率矩阵

转移概率矩阵可以形象、直观地描述不同土地利用类型之间的相互转化情况^[18]。对两期土地利用数据做空间叠加分析, 可以得到土地利用类型间数量关系的原始转移矩阵A_ij。为了分析乌双流域湿地的损失情况, 由原始转移矩阵A_ij可以进一步计算出湿地的转移概率矩阵, 公式为：

式中，A_ij表示k时期i种土地利用类型转变为k+1时期j种土地利用类型的面积, B_ij表示k时期i种土地利用类型转变为k+1时期j种土地利用类型的比例。

2.2.3 因子选择和数据提取

(1) 因子选择：研究区1980年以后增温速率明显偏高^[19-21], 提升了湿地水体的潜在蒸发量；降水通过水资源的时空分配直接影响湿地形成和发育^[22-23]；地形起伏通过决定径流走向影响湿地空间分布格局, 同时也决定了湿地的开发难度^[24]；贯通性的工程建设会切断了湿地的天然水文联系, 导致湿地景观破碎化, 增加退化风险^[25]；耕地侵占湿地往往是从湿地外围开始的, 该部分湿地受人类活动干扰最强, 且开发难度低, 所以极易被开垦^[26]；建立自然保护区被国内外实践证明是保护典型湿地生态系统的有效措施^[27-28]；此外人口分布和产业发展对湿地的影响不容忽视。鉴于上述, 本文选取了7个驱动因子, 分别为：1980—2015年温度变化、1980—2015年平均降水量、数字地面高程DEM、距道路距离、距1980年耕地距离、自然保护区的建立、人口经济发展。

(2) 数据提取：将湿地区域划分为1 km×1 km的网格单元, 网格中心点作为样本点, 共4492个。1980—2015年温度变化为2015年气温三年滑动平均值与1980年气温三年滑动平均值的差值；距道路距离和距1980年耕地距离通过近邻分析方法得到；对自然保护区数据进行二值处理, 保护区内赋值为1, 非保护区赋值为0；人口经济发展状况用多年平均的DMSP/OLS夜间灯光数据值模拟。

2.2.4 增强回归树

增强回归树(Boosted Regression Trees, BRT)是一种结合统计和机器学习两种算法的模型, 分别为分类回归树和boosting算法, 其核心思想是通过不断地随机选择和自我学习生成多重回归树, 以此提高模型的分辨能力和稳定性能^[29]。BRT的最大优势在于不必考虑自变量之间的交互作用, 能适应复杂的非线性关系^[30], 输出的自变量贡献度及影响曲线比较直观, 易于解释。目前BRT在生物分布模拟^[31]、城市扩张^[32]、野火分布格局^[33]、流域污染^[34]、大气污染^[35]等研究中得到了一定的应用, 也取得了较好的效果。

将乌双流域1980年与2015年的湿地空间分布图叠加, 得到1980—2015年间的湿地变化值(0、1值), 1值表示由1980年的湿地区域转出为2015年的非湿地区域, 0值表示湿地未改变的区域。以1980—2015年湿地变化值作为因变量, 选取的7个因子作为自变量, 调用BRT程序包进行增强回归树分析。设置学习速率为0.005, 每次抽取50%的数据进行分析, 50%用于训练, 并进行5次交叉验证。

3 结果与分析 3.1 湿地动态变化分析

乌双流域湿地面积呈持续减少趋势(表 1), 1980—2015年间湿地面积减少了738 km², 较1980年减少了16.43%。沼泽地作为乌双流域湿地的主要组成部分, 占比始终在75%以上, 其面积也呈持续减少趋势, 1980—2015年间减少了667 km², 较1980年减少了19.02%, 且沼泽地占湿地比例持续下降, 从1980年的78.07%下降到2015年的75.76%。水域在湿地中的占比始终在20%以上, 其面积在1990—2005年间出现较大波动, 2005年后趋向稳定。

湿地与沼泽地的面积变化率在各时段均为负值, 而水域面积在不同时段变化方向不同(图 2)。1990—1995年间湿地和沼泽地面积减少速率最大, 约为-1.5%/a, 其次是2000—2005年, 约为-1.0%/a。水域面积在1990—1995年和2000—2005年间减少, 变化率分别约为-2.5%/a和-0.8%/a, 在1995—2000年间增加, 变化率约为1.5%/a。其余时段各湿地类型的面积变化率相对较小。整体上, 1980—2015年间湿地、沼泽地和水域面积都有所减少, 面积变化率分别为-0.66%/a、-0.76%/a和-0.22%/a。

3.2 湿地景观格局变化分析 3.2.1 类型水平

由图 3可见，1990—1995年间沼泽地的斑块数量由241个锐减到188个, 平均斑块面积有所增加, 结合沼泽地的聚集度指数在此期间降低, 表明1990—1995年间丧失了部分面积较小的沼泽斑块；其他时段内沼泽地斑块数量上升, 平均斑块面积下降, 表明沼泽地的破碎程度加剧。水域斑块数量在1980—2015年间呈减少趋势, 平均斑块面积除1995年偏低之外维持基本稳定。

面积加权平均形状因子反映了斑块类型的形状特征, 值越小说明斑块形状越稳定, 值越大说明斑块形状变得更复杂, 更不规则。由图 3可见, 沼泽地的聚集度指数在1990年后持续减少, 说明沼泽地的空间分布模式在1990年后趋于离散, 斑块间关联性持续降低；沼泽地的面积加权平均形状因子在1980—1990年间升高后一直维持在一个较高水平, 说明沼泽地可能受到外部人为活动干扰后, 斑块形状变得复杂且不规则, 处于不稳定状态。水域的聚集度指数和面积加权平均形状因子在1995年存在一个“低谷”, 可以说明1995年前后水域斑块有一个先离散再集聚的过程, 结合水域斑块的平均斑块面积在1995年偏低, 表明离散的原因是水域斑块萎缩导致斑块间距离增大, 同时斑块形状趋于规则和胁迫下的稳定状态。

3.2.2 景观水平

由图 4可见, 湿地的最大斑块指数在1980—2000年之前波动变化, 1990年升到最高值, 1995年降到最低值；2000—2015年间最大斑块指数持续降低, 说明2000年后最大斑块的优势度在整个湿地景观中的地位呈下降趋势。湿地的斑块密度除1990—1995年间下降外, 其余时段内均有所升高, 说明湿地斑块的破碎化程度加剧。

聚集度指数在1990年后持续降低, 说明湿地的空间分布愈发趋于离散。斑块结合度指数反映了景观类型的自然连通性。1980—1990年间斑块结合度指数增加, 湿地斑块的连通性增强；1990年后随着斑块结合度指数持续降低, 湿地斑块的连通性逐渐减弱。

3.3 湿地空间变化分析 3.3.1 湿地变化的空间格局

1980—2015年间乌双流域湿地损失集中分布在乌裕尔河两岸、双阳河中游及流域尾闾(图 5)。其中乌裕尔河中游、双阳河上游湿地损失面积相对较大且集中, 其他区域的损失面积相对较小且分布在湿地边缘。此外, 湿地损失的同时也有少量区域由非湿地转为湿地, 但面积较小。

湿地变化的空间格局存在时间尺度上的差异。1980—1990年间, 湿地损失主要分布在乌裕尔河两岸、双阳河中游及流域尾闾。1990—1995年间, 湿地损失程度最为严重, 乌裕尔河中游两岸大面积的湿地损失, 流域尾闾的湿地损失分布也较广。1995—2000年间, 湿地损失的集中分布区域同上一个五年, 但损失程度有所减缓。2000—2005年间, 湿地丧失区域集中在流域尾闾。2005—2010年与2010—2015年两个时段湿地损失范围较小。

3.3.2 湿地转移概况

1980—2015年间乌双流域湿地主要转出为水田、旱地、草地和盐碱地, 总的转移概率约为20%(图 6)。湿地转出集中发生在2000年之前的3个时段里, 总的转移概率均超过20%, 1990—1995年间甚至达到30%, 主要转出为旱地、草地、林地和盐碱地, 少量转为居民地；1990—1995年间开始有部分湿地转为水田, 这一转移概率在1995—2000年间进一步升高。2000—2005年间湿地转移概率明显减小, 但依然有部分湿地转为草地、旱地与盐碱地。2005—2010年间湿地几乎未发生转移。2010—2015年间再次出现湿地转为水田的情况。

3.4 湿地变化驱动机制分析 3.4.1 因子贡献率

BRT分析的因子贡献率结果如图 7所示, 影响流域湿地转出最大的因子是1980—2015年平均降水量, 贡献率为28.66%；其次是1980—2015年温度变化, 贡献率为22.16%；其余因子依据贡献率从大到小分别是距道路距离、DEM、距1980年耕地距离、自然保护区建立和人口经济发展。

气候因素是导致乌双流域1980—2015年湿地转出的主要因素, 平均降水量和温度变化两个因子总的贡献率达到50%, 说明流域内湿地转为非湿地50%的概率是由于气候因素造成的。距道路距离的贡献率为15.75%, 道路的修建会导致湿地斑块的破碎化, 更易造成湿地的萎缩。DEM的贡献率为15.61%, 对湿地转出存在一定影响。用距1980年耕地距离来量化耕地扩张对湿地转出的影响, 其贡献率为10.80%。自然保护区的建立对湿地转出的贡献率为3.99%, 即自然保护区范围内发生湿地转出的概率较低, 说明建立自然保护区对于湿地保护起到一定的积极作用。人口经济因子对湿地转出的贡献率较低, 为3.02%, 说明人口分布及产业发展对湿地转出的影响较小。

3.4.2 主导因子影响机制

选取贡献率排名前6位的驱动因子, 绘制其影响力变化图(图 8)。图中的曲线表示随驱动因子取值的变动, 其对湿地转出影响程度的变化, 边际效应数值越大即代表影响程度越大。

降水量与湿地转出大致呈负相关, 影响程度可大致分为3个层级(图 8)。多年平均降水量小于450 mm对湿地转出影响程度最大, 最高值出现在430 mm左右；多年平均降水量在450—500 mm之间对湿地转出的影响程度次之；多年平均降水量大于500 mm对湿地转出的影响程度较弱。

温度变化对湿地转出的影响随增温幅度的提高呈现正相关(图 8)。增温小于1.4℃对湿地转出的影响程度相对较小, 而当增温超过1.4℃, 影响程度陡然上升, 在1.5 ℃左右影响程度达到最大。

距离道路越近, 湿地转出风险越大(图 8)。由其影响曲线可以看出道路修建对湿地转出在一定距离范围内具有明显的影响趋势, 这一距离范围约为1.8 km, 超出1.8 km影响程度会逐渐减弱。

高程的提升会增大湿地转出的风险(图 8)。由DEM的影响曲线, 高程达到170 m左右, 其对湿地转出的影响程度达到最大, 超过170 m影响程度保持相对平稳, 说明相对高海拔处的湿地受地形影响更易损失。

距离耕地距离越近, 湿地越易受农业生产影响而转出(图 8)。由其影响曲线可以看出具有明显影响趋势的距离范围约为2 km, 超出2 km影响程度相对较小, 表明耕地周围2 km的湿地存在被开垦的风险。

相较于自然保护区, 非自然保护区内的湿地转出风险更大(图 8)。由其影响程度图可以看出, 非自然保护区对湿地转出的影响程度比自然保护区高出约一倍, 说明自然保护区的建立确实有效遏制了湿地损失。

3.4.3 机制验证及湿地变化预警

基于湿地变化驱动力研究可以进一步预测未来湿地的演变趋势。本文根据湿地转出主导因子的影响机制, 在空间上预测未来极易发生湿地损失的区域, 即湿地变化预警区域, 从而为湿地管理与保护政策的制定提供科学参考。

在绘制湿地变化预警图之前, 首先结合1980—2015年湿地变化情况对驱动机制结果中所得到的阈值进行验证。由于湿地变化的各驱动因子并非孤立产生影响, 而是共同发挥作用, 因此首先提取综合考虑主导因子影响机制的空间范围, 其同时满足以下六个条件：(1)多年平均降水量低于450 mm；(2)温度变化值大于1.4℃；(3)位于道路1.8 km缓冲区内；(4)高程在170 m以上；(5)位于1980年耕地2 km缓冲区内；(6)排除多年平均降水量高于500 mm的区域。结果显示, 1980—2015年间87.3%湿地转出区域落在提取的空间范围内, 由此证明驱动机制的结果真实有效。

基于2015年的湿地范围和主导因子的影响机制绘制乌双流域湿地变化预警图, 如图 9所示。湿地预警区域集中分布在乌裕尔河中游、双阳河中下游及流域尾闾, 这部分湿地在当前驱动因子条件下更易受到影响而损失, 故在未来的湿地保护工作中需重点关注。

4 讨论

影响湿地变化的因素复杂多样, 在建立驱动因子指标体系时常用自然因子和社会经济因子来概括, 两类因子相互交织作用致使湿地景观格局发生变化^[36]。本文综合考虑了影响乌双流域湿地变化的自然因子与社会经济因子, 同时加入政策因素如自然保护区的建立, 揭示了湿地变化背后的原因, 为乌双流域湿地保护提供决策依据。结果显示气候变化在乌双流域湿地景观格局演变过程中起到关键性作用, 这与佟守正等分析扎龙湿地生态变化影响因子的结论相一致^[4]。降水直接决定流域水资源的时空分配与水文循环, 而增温导致湿地水体消耗加剧, 造成湿地萎缩, 同时还引起温度为限制因子的农作物适宜种植范围的扩大^[37], 间接导致湿地被开垦。但由于数据获取及处理等限制, 驱动力因子指标体系还有待于进一步完善, 如社会经济方面, 除湿地开垦外, 占据水源优势的农民还热衷于将旱地改为水田, 由土地利用转移分析发现, 1980—2015年间有254.7 km²的旱地转出为水田, 这导致上游地区的农业用水量增加, 同时生产生活用水截留也减少了进入湿地的地表水补给量, 从而造成部分沼泽湿地因水资源短缺而变为干草地甚至盐碱地。

如何实现因子数据的空间化是基于空间尺度进行驱动力分析首先需要面临的难题。在处理大多数社会经济因子时, 按照行政单元划分进行因子赋值是常用的方法之一, 如Wang等^[38]在分析某县神经管畸形出生缺陷的环境风险因子时, 将各村的统计发病率和因子数据赋给各村的行政区域后做驱动因子分析。但该处理方式也存在一定局限性, 一是按行政单元赋值会忽视行政单元内部的分异性, 仅适用于研究尺度大于数据尺度的情况, 二是不适用研究区域边界与行政边界不重合的情况, 这一现象在生态学研究中较为普遍。另一个解决思路是用相关性较强的其他数据进行代替或模拟。如本文尝试用夜间灯光数据来反映流域的人口经济情况, 其作用与刘吉平等^[39]通过不同的土地利用类型计算人为干扰度相似。夜间灯光值强弱与人口经济因子间存在正相关关系, 如何利用这种定量关系推算出单一的社会经济因子是驱动力研究进一步的思路。此外, 在量化道路修建和耕地扩张这两个因子时所用到的邻近分析方法, 也为因子空间化处理提供一个新的思路。

政府决策对于流域短期的水文情势具有强制性的影响。以乌双流域为例, 黑龙江西部地区频频遭遇的干旱一度导致乌裕尔河、双阳河的径流量明显减少, 再加上流域上游生活用水的不断增加, 使扎龙保护区常常面临缺水状况^[40]。为解决保护区常年遭遇的“渴水”状况, 水利部和黑龙江省政府从2001年开始向扎龙湿地补水, 并实施了“引嫩”工程从松花江上游的嫩江引水补充扎龙湿地储水, 至2013年6月, 已累计向扎龙湿地补水19.52亿m³, 最多年份补水3.43亿m³, 使扎龙湿地的缺水危机得到明显缓解, 湿地损失程度减缓。因此在今后的驱动力研究中需要考虑补水量对短期湿地变化的直接影响, 针对典型补水区域的湿地变化, 对引水工程的影响进行定量评估。

5 结论

(1) 1980—2015年间乌双流域湿地面积持续减少, 减少了738 km², 较1980年减少了16.43%；沼泽地占湿地的75%以上, 面积也持续减少, 且在湿地中的占比也在下降；水域面积在1990—2005年间经历较大的波动后趋于稳定；1990—1995年是湿地面积减少速率最大的时段。

(2) 1980—2015年间湿地景观格局的变化趋势是最大斑块的优势度地位下降, 破碎化程度加剧, 空间分布离散, 连通性减弱；沼泽地斑块除1990—1995年间数量锐减外, 其他时段斑块数量上升, 破碎度增加, 集聚度和斑块稳定性下降；水域斑块数量呈减少趋势, 1995年前后有一个先离散再集聚的过程, 离散的原因是水域斑块萎缩导致斑块间距离增大。

(3) 1980—2015年间乌双流域湿地损失集中分布在乌裕尔河两岸、双阳河中游及流域尾闾, 其中乌裕尔河两岸、双阳河中游湿地损失面积较集中, 流域尾闾的损失面积较分散, 主要在湿地边缘；湿地主要转出为水田、旱地、草地和盐碱地, 总的转移概率约为20%, 且集中发生在2000年之前。

(4) 气候因素是导致乌双流域湿地损失的主要因素, 降水量和温度变化两个因子总的贡献率达到50%, 其中降水量与湿地转出呈负相关, 小于450 mm对湿地损失影响最大, 温度变化与湿地转出大致呈正相关, 温度变化值超过1.4℃对湿地损失影响最大；道路修建对湿地损失的影响范围约为1.8 km；高海拔会增加湿地损失的风险；耕地周围2 km的湿地存在被开垦的风险；自然保护区的建立有效遏制了湿地的损失。未来乌裕尔河中游、双阳河中下游及流域尾闾的湿地存在损失风险。