景观格局是不同景观镶嵌体在空间上的排列^[1], 是不同生态系统组成的地理综合体, 能有效揭示区域生态状况和空间变异特征^[2]。景观格局变化主要是人类活动通过土地利用变化驱动的^[3]。景观的空间配置与类型组成影响水循环, 从而影响流域系统的水文过程^[4-5]。从景观生态学角度探讨流域景观变化的水沙过程响应, 能够揭示人类活动对流域水沙的影响。在水资源短缺和水土流失严重的区域, 开展此研究, 可以为区域土地利用的合理调整、水土保持治理与生态恢复提供科学依据。目前, 景观格局变化的水沙过程响应已成为研究热点。Dixon与Earls^[6]发现佛罗里达州查理溪流域径流随城市扩张而增加；王计平等发现斑块密度是影响黄土丘陵沟壑区产沙的重要指标^[1]；黄志霖等发现三峡库区景观格局对径流和泥沙影响均显著^[7]；郭军庭等发现潮河流域灌木林地和草地能够增加径流^[8]；林炳青等发现晋江流域洪水径流与景观格局相关性显著^[4]；李晶与周自翔发现延河流域景观格局变化与产沙显著相关^[9]；郝振纯等分析土地利用情景变化对海河流域典型区域径流的影响^[10]。

滦河流域距京津地区约240km, 是京津沙尘暴的重要通道和沙源地之一^[11]。北京沙尘绝大部分源自外地, 与滦河流域的丰宁和多伦的土壤属性非常相似^[12]。为优化京津地区环境质量, 2000年京津风沙源治理工程启动, 多伦县为工程发端地^[13]。滦河流域位于工程治理区内^[14], 主要通过土地利用变化影响地表水下渗、蒸散 (发) 等水文过程, 因此, 探讨景观格局变化的水文效应, 对滦河流域水土流失治理具有重要意义。目前, 滦河流域产流产沙相关研究得到广泛关注, 研究多集中于径流, 研究范围主要为潘家口水库上游^[15]与全流域^[16]；泥沙研究主要针对潘家口水库上游流域^[15], 较少涉及全流域；水文影响因素主要集中于土地利用与气候变化对径流的影响^[17], 较少从景观格局变化分析产流产沙。

本文运用景观格局指数分析1980—2010年滦河流域景观格局变化特征, 结合SWAT模型探讨景观格局变化的流域水沙过程响应；应用Spearman相关分析法分析景观格局变化对水沙过程的影响；应用SWAT模型模拟1976—2012年滦河流域径流与泥沙, 分析流域2000年治理工程实施前后多年年均产流与产沙的时空格局变化, 识别流域仍需重点治理的产沙关键源区, 以期为滦河流域的后续治理与土地利用规划提供科学依据。

1 研究区域、方法与数据来源 1.1 研究区域

滦河流域位于华北北部, 滦河全长888km, 流域总面积为44750km²(图 1)。滦河发源于河北丰宁, 流经内蒙古、河北及辽宁的7市 (盟) 的27个县 (区)^[18]。滦河流域属于典型温带大陆性季风气候, 年均气温5—12℃, 多年平均降水量为400—700mm, 主要集中分布于7—8月。滦河流域地形起伏较大, 由西北向东南倾斜, 分为坝上高原、燕山山地、南部平原3种地貌类型。流域内主要土地利用类型为林地、草地和农业用地。流域土壤类型丰富, 主要包括草甸土、森林土、栗钙土、棕壤、潮土和褐土等^[19]。

1.2 数据来源与处理

DEM：90m分辨率, 来源于美国地质勘探局 (USGS)。根据流域情况与相关研究^{[15, 17]}, 设定最小集水面积阈值为20000ha, 流域划分为116个子流域, 以子流域为基础单元研究产水产沙的空间分布特征。根据土壤、土地利用、坡度的分布特征, 面积阈值分别设定为15%、10%、15%, 生成1521个HRU, 充分体现了滦河流域下垫面的空间不均匀性。

土地利用数据：1980s土地利用图, 分辨率1km, 来源于地球系统科学数据共享平台；2010年土地利用数据源于分辨率30m的TM数据 (轨道号为121/32, 121/33, 122/31, 122/32, 123/31, 124/31), 来源于USGS。运用ENVI 4.8和ArcGIS 9.3, 经过镶嵌、拼接等处理, 采用人机交互解译影像, 并重采样为1km分辨率。依据土地资源分类系统^[20], 采用国家标准体系, 分为林地、草地、农业用地、水域、建设用地和未利用地六类。依据数据可获得性、研究目的以及20世纪80到90年代土地利用变化相对缓慢^[21], 采用1980s土地利用数据用于1976—1999年时间段的模拟；同样地, 2010年土地利用数据用于2000—2012年时间段的模拟。

土壤数据：1：100万的世界土壤数据库 (HWSD) 来源于寒区旱区科学数据中心。土壤水文分组根据在完全湿润并不冻的条件下土壤最小渗透属性来确定^[22]；土壤容重等采用SPAW软件进行计算。

气象数据：来源于中国气象科学数据共享服务网, 包括多伦、丰宁、围场、承德、乐亭5个气象站1961—2012年的日降水、气温、风速等。利用上述气象数据建立天气发生器, 输入数据由pcpSWAT和dewSWAT软件得出^[8]。

水文数据：1976—1991年承德、下板城、李营、宽城、滦县5个水文站的月径流量、潘家口和大黑汀水库月出流数据和滦县站泥沙数据, 来源于滦河流域水文年鉴等资料。

SWAT模型空间数据统一采用Albers等积投影。借鉴相关研究^[15-17], 90m分辨率DEM满足流域水系的正确划分；1km分辨率土地利用图与土壤图适用于流域水沙过程模拟。

1.3 研究方法 1.3.1 SWAT模型

SWAT模型的开发目的是模拟与预测具有复杂变化的土壤、土地利用和管理方式的较大流域中土地利用及管理方式对流域产流、产沙及污染物负荷的影响。文中选用SCS径流曲线数方法计算地表径流, 泥沙由改进的MUSLE方法计算, 采用Penman-Monteirh方法计算潜在蒸发量, 马斯京根方法进行河道演算。采用SWAT2009的敏感性分析模块中全局参数敏感性分析方法 (LH-OAT) 进行径流和泥沙参数敏感性分析^[23]。运用SWAT-CUP软件的t检验方法进行参数的深层筛选^[24]；依次校准径流与泥沙, 通过SWAT-CUP的SUIF-2方法^[24]对径流进行初步率定后手动调参；泥沙校准采用手动调参。

Nash-Suttcliffe系数 (Ens), 决定系数 (R²) 和相对误差 (Re) 为模型适用性评价指标。当Ens>0.5, R²>0.5, Re < 20%时, 径流模拟结果是满意的；Re < 55%时, 泥沙模拟结果是满意的^[25]。

1974—1975年为预热期, 1976—1983年作为模型校准期, 1984—1991年为验证期。利用1992—2012年气象数据、1980s和2010年土地利用数据和已率定参数, 得到1992—2012年径流与泥沙模拟数据。借鉴相关研究经验^[17]：保持气象、土壤等输入条件不变, 仅改变土地利用条件时, 水文模型输出结果差异仅是由土地利用变化造成的。本文以相同的气象、地形、土壤、水库出流等条件, 分别输入1980s和2010年土地利用图, 模拟和分析景观格局变化的水沙响应, 评价工程治理前后流域水土流失变化。

针对上述假设, 采用2000—2012年气象数据, 满足假设中相同气象条件。研究区面积较大, 2000—2012年全流域尺度地形地貌变化、土壤类型变化皆可忽略, 满足假设中相同的地貌与土壤条件。滦河干流建设了潘家口和大黑汀水库, 支流的小型水库对流域整体水文过程影响较小, 因此主要考虑潘、大水库 (分别由1981年、1983年开始运行至今^{[16, 26]}), 满足相同的水利工程条件。可见假设合理, 研究方法可行。

1.3.2 景观格局分析

运用ENVI 4.8计算景观类型面积变化特征；利用FRAGSTATS 3.3软件分析流域景观格局指数变化；应用Spearman相关分析法^[4], 分析流域水沙过程变化与景观格局指数变化的相关关系。

考虑相关文献中对水沙影响较大、意义明确且相对独立的景观指数^[4], 本文选取斑块类型面积 (CA)、斑块密度 (PD)、景观形状指数 (LSI)、Shannon多样性指数 (SHDI)、Shannon均匀度指数 (SHEI)、分离度 (SPLIT)、最大斑块指数 (LPI) 和蔓延度 (CONTAG), 分别反映景观类型面积、形状复杂度、破碎化程度、多样性与连通性。

统计1980年、2010年全流域和各子流域景观指数, 并计算全流域和各子流域1980—2010年景观格局的变化率, 计算公式为：

(1)

式中, P_i为i类景观格局指数变化率, F_{i, 2010}为2010年景观格局下i类景观格局指数, F_{i, 1980}为1980s景观格局下i类景观格局指数。

2 结果与分析 2.1 模型的校准与验证

如表 1所示, 各水文站控制的子流域中敏感性因子对径流的敏感度存在差异, 这与研究区面积较大, 地形、土地利用与土壤类型多样有关。然而在所有子流域中, 径流曲线系数CN2的敏感度均较高, 这与其他研究相一致^[15-17], 此参数直接影响径流的大小, 值越高径流量越大。

模型校准与验证结果见图 2、表 2, Ens和R²均达到0.7以上。径流模拟效果较好, 泥沙模拟与实测值的过程线也趋于一致, 但丰水年的泥沙模拟值低于实测值, 特别是洪水期。这主要由于SWAT模型在产沙量高于2000t时会产生低于实际值的模拟结果^[27], 而滦河流域年均产沙量约为112万t, 远高于2000t。此外, 强降水时段模拟值也会低于实测值^[28], 尤其是洪峰期的误差最大^[29]。

综上, SWAT模型模拟径流与泥沙的精度均达到要求。滦河流域面积较大, 研究区上游的滦河源头为冬春季的风蚀与夏季的水蚀互为动力, 但水蚀仍是全流域最主要的土壤侵蚀过程^[30], SWAT模型适用于土壤、土地利用等具有复杂变化的滦河流域水沙模拟研究。

2.2 治理工程前后产流与产沙的时空分布特征 2.2.1 时间分布特征

2000年京津风沙源治理以来, 年均径流量和产沙量比1976—1999年分别减少74.26%与90.25%, 而降水量仅减少3.84%(图 3), 表明人类活动 (流域治理) 是影响滦河流域径流与泥沙的主要因素。这是由于影响径流量变化的主要因素为气候变化和人类活动, 气候因素主要表现为降水变化, 人类活动主要为消耗用水和土地利用变化^[21]。相关研究也认为流域治理是径流量减少的重要原因^[16], 其中以水保措施改变下垫面的影响量约为消耗用水的2倍^[31]。

工程开展前后, 年径流与泥沙均主要集中于滦河流域的汛期 (6—9月), 7—8月径流与泥沙量最大, 约为年径流与泥沙量的44%和84%, 治理后7—8月的洪峰径流显著下降, 工程实施对年内洪峰流量的调蓄效果显著 (图 4)。

2.2.2 空间分布特征

工程实施后, 产水与产沙的高值区显著缩小 (图 5、图 6)。1976—1999年, 下游多数区域的径流深大于150mm；2000—2012年, 只有位于下游入海口附近的3个子流域径流深大于150mm, 而中上游出现年均产水量小于50mm的区域, 如沽源县 (34、35号子流域) 等 (图 5)。

如图 6所示, 工程实施前, 流域年均产沙模数大于10 t/hm²的区域位于下游的青龙县；工程实施后, 年均产沙模数大于10 t/hm²的区域已消失, 而2000年以后流域中上游的大部分区域年均产沙模数小于1 t/hm²。仍需重点治理的区域主要为滦河上游的沽源 (35号子流域)、丰宁 (33、38号子流域)、围场 (28号子流域) 与隆化 (42号子流域)；中游的滦平 (71号子流域)、承德市 (71、74号子流域)、承德县 (75、77号子流域)、平泉 (68、75号子流域)、凌源 (68号子流域)；下游的迁安与卢龙 (114号子流域)。

2.4 滦河流域治理与景观格局变化

与20世纪80年代相比, 2010年林地、建设用地分别增加41.40%、577.82%；农业用地、草地、未利用地和水域分别减少26.99%、41.76%、0.50%与36.88%(表 3)。这与流域治理工程关系密切, 1980—1989年, 流域上游开展人工造林4700km², 建设草场200km^{2 [32]}；2000年国家水土保持重点项目京津风沙源治理工程启动后, 采取人工造林、改良草场、提高植被覆盖度等措施建立高质量的流域生态防护, 治理成绩显著, 如多伦县造林565km^{2 [13]}；正蓝旗造林262km²、草地治理445km^{2 [33]}；承德县造林148 km², 治理252km^{2 [34]}；丰宁造林138km², 治理248km^{2 [35]}；围场退耕还林665km^{2 [36]}；滦平县森林覆盖度提高15.4%, 治理面积950km^{2 [37]}等。

与80年代相比, 2010年流域分离度和形状复杂度显著降低, 景观多样性和均匀度略有降低；斑块密度变化不大；最大斑块指数、蔓延度显著上升 (表 4), 说明林草建设与修复工程促进流域景观破碎化降低, 景观分布更集中, 斑块形状更规则, 景观连通性提高、景观类型向非均衡方向发展, 林地景观优势更显著。

2.5 景观格局变化的水沙过程响应

如图 7所示, 与20世纪80年代相比, 2010年流域年均地表径流减少9mm；109个子流域地表径流减小, 面积约39971km², 占全流域的89.32%。流域年均产水量增加5.44mm；85个子流域产水量增加, 面积约34326.18km², 占全流域的76.71%。流域年均产沙减小1.59t/hm²；104个子流域产沙减小, 面积约42015.37km², 占全流域的93.89%。可见, 滦河流域水土保持治理能够很好地调节地表径流, 减少泥沙。

2.5.1 景观格局变化的地表径流响应

与20世纪80年代相比, 2010年地表径流减少最显著的区域主要位于上游16、35号子流域, 中游63、74号子流域, 下游97、98、101号等子流域 (图 7、表 5)。这些区域地表径流减少与林草植被建设密切相关。例如63、74号子流域位于承德县境内, 实施京津风沙源治理工程以来, 承德县调整产业结构退耕还林还果、造林育林, 全县森林覆盖率达到61%^[34], 有效拦蓄降水, 减少地表径流。

2010年仅12个子流域地表径流增加 (图 7)。如113号子流域, 位于迁安境内, 地表径流增加最显著。迁安市经济依托于铁矿资源, 矿山排土场和尾矿库在堆放过程中导致绿化率较低。虽然进行综合治理, 但是由于植物种类与苗龄选择不当、栽植方式与时间不合理, 因而植被覆盖率仍然较低, 而且河道两侧山体上的排土场无法进行清除, 难以采取水保措施^[38], 水土流失愈加严重, 因而地表径流增加。

2.5.2 景观格局变化的产水响应

与1980年相比, 2010年产水量增加最显著区域主要位于研究区上游2、3、6、10、13号子流域, 中游66、107、108号子流域和下游106、109、111、113—115号子流域 (图 7、表 5)。产水减少区域主要位于滦河上游的16、17、19号等子流域, 中游的48、51、58号等子流域, 下游的101、104、105、116号子流域。在产水量增加的子流域中, 如1—9等子流域地表水是减少的 (图 7), 这些区域多数由于水土保持治理程度提高, 地表径流减少时大量转化补充土壤含水量与地下径流。在产水量减少的子流域中, 如流域上游的围场、沽源与丰宁境内16、17、29号等子流域地表水也是减少的 (图 7), 这是由于沽源县植被覆盖以草地为主、丰宁与围场的林地建设仍需继续开展, 地表径流仍为产水的主要组成部分, 壤中流与地下径流的增加程度尚显不足。

2.5.3 景观格局变化的产沙响应

与20世纪80年代相比, 2010年景观格局下, 产沙减少相对较显著的区域主要位于滦河流域上游的21、23、30、35号子流域, 中游63号子流域和下游的105、111号子流域 (图 7、表 5)。结果表明, 这些区域水土保持治理效果最显著, 土壤侵蚀程度大大降低。例如63号子流域, 位于承德县东北部, 产沙量与地表径流明显减少, 这与承德全县森林覆盖率增加有关^[34], 森林覆盖率增加能够改善地表径流、抑制产沙。

产沙增加的子流域仅12个, 最显著的位于流域上游的20、32、40号子流域, 中游75、85号子流域 (图 7、表 5)。2010年, 20、32、40号子流域林地减少被农业用地占用；75、85号子流域草地减少被农业用地占用。林地、草地减沙, 农业用地增沙, 因而水土流失严重、产沙增加。

2.6 景观格局对水沙过程的影响

景观类型与水沙过程相关性见表 6, 林地与地表径流、产水量和产沙量呈负相关关系, 说明滦河流域林地面积增加对径流和泥沙都具有积极抑制作用, 这是由于森林植被的下渗与蒸发量较大, 可以有效的调节地表径流与洪峰流量, 且土壤含水量的增加也将减少泥沙的产生与输移^[39]。草地面积与产水量呈正相关, 与产沙呈负相关, 表明草地面积的增加会促进产流、减少产沙, 主要由于草地能将坡面径流流速减小40%、将起流时间推迟3倍, 增加径流入渗率、土壤含水量和地下水, 从而增加产流^[40]；草地的植被根系存在于土壤表层, 提高土壤对径流侵蚀产沙的抵抗力。农业用地与地表径流、产沙量均呈正相关, 由于农业用地的土壤结构被破坏, 降低了径流下渗速率、土壤含水量和地下水补给^[41], 促使地表径流增加；而耕地植被覆盖度较低, 水土流失也相对容易。建设用地与未利用地与产沙相关性较弱, 可能与两者面积小有关。

地表径流、产水、产沙与景观形状复杂程度、Shannon′s均匀度、景观分离度呈正相关, 与蔓延度、最大斑块指数呈负相关；景观斑块密度和Shannon′s多样性指数与产水产沙呈正相关 (表 6)。表明大力营造林草地后, 景观形状规则化、连通度增加、流域优势景观类型林地增加, 导致流域降水截留能力与水源涵养能力增强, 地表径流、产水和产沙降低。

3 结论

(1) 与80年代相比, 2010年滦河流域林地、建设用地增加, 其他用地类型减小；流域斑块密度变化不明显, 景观集中连片度提高、连通性变优、优势斑块类型明显、破碎化降低、景观形状越来越规则、多样性程度减少, 景观类型向非均衡方向发展。

(2) 不同景观格局下, 滦河流域水沙过程差异显著。与1980年相比, 2010年流域年均地表径流减少9mm, 产水量增加5.44mm, 产沙减小1.59t/hm²。年均地表径流减少区域占全流域的89.32%, 主要位于流域上游的沽源和围场, 中游的承德县和承德市双桥区, 下游的青龙、建昌、迁安与入海口的三角洲地区；流域年均产水量增加区域占全流域76.71%, 主要位于流域上游的正蓝旗、克什克腾旗、多伦, 中游的承德市双桥区和宽城, 下游的迁西、青龙、迁安、卢龙与滦县；年均产沙减少区域占全流域的93.89%, 位于上游沽源和围场、中游承德县与下游的青龙和迁西。

(3) 地表径流、产水、产沙与林地面积呈负相关；产水、产沙与草地面积分别呈正、负相关；地表径流、产沙与农业用地呈正相关；地表径流、产水、产沙与景观形状复杂程度、Shannon′s均匀度、分离度呈正相关, 与蔓延度、最大斑块指数呈负相关；产水产沙与景观斑块密度、Shannon′s多样性指数呈正相关。

(4) 京津风沙源治理工程开展以来, 年均径流量和产沙量显著减少, 汛期径流与泥沙明显降低；产水与产沙高值区显著缩小。中上游大部分区域年均产沙小于1 t/hm², 大于10 t/hm²的区域已消失。产沙关键源区仍需重点治理, 应继续加强林草建设, 提高植被质量与覆盖率；对陡坡耕地实施退耕还林, 对严重沙化耕地实施退耕还草；集约化用地等；更应重视有利于防治泥沙输出的景观类型 (如林地、草地) 的空间聚集程度、延展性、连接度与规则形状。