黄河流域多时空尺度土地利用与水质的关系*

郭羽羽，李思悦，刘 睿，张 静**

(1：重庆师范大学地理与旅游学院，重庆 401331) (2：中国科学院重庆绿色智能技术研究院生态水文研究中心，重庆 400714)

水资源是人类生存和发展所依赖的重要资源，随着社会的不断进步，水资源分配不均的问题日渐凸显[1]. 水质是水资源质量的决定因素，点源污染目前已基本得到控制，非点源污染成为影响水质的关键[2]. 土地利用变化对生态过程、地表径流、水文循环产生一系列深刻的影响，进而影响河流水质安全[3]. 因此，探讨土地利用与水质之间的关系，对土地利用的管理规划及水生态环境资源的有效保护具有重要意义.

近年来，土地利用对水质指标的影响机制研究已成为国内外学者研究的热点[4-8]. 目前，大部分的研究是针对河段[9-11]、河岸带[12-13]、子流域[14-16]等空间尺度，研究方法主要包括相关分析[17]、冗余分析[11,13]、SWAT模型[18]、地理加权回归模型[19]、多元线性回归[20]等. 例如，杨洁等[9]认为在1000 m大范围尺度，土地利用对水质的影响更加显著，城镇用地、旱地及水田影响显著，景观聚集程度与水质指标呈正相关，与破碎度呈负相关；Shi等[10]研究表明在1000 m缓冲区尺度下，土地利用与河流水质的相关性最高，且不同的土地利用指标对水质有不同的尺度效应；张殷俊等[12]研究了巢湖流域缓冲区尺度下土地利用与水质之间的关系及影响因素，结果表明土地利用类型对水环境影响具有尺度效应，林地对水质具有正面影响，建设用地具有负面影响，耕地对水质的影响具有复杂性；汪昱昆等[13]报道了上海市河网水质对河岸带土地利用的响应，结果发现土地利用在500 m和1000 m缓冲区尺度对水质解释最强，建设用地与大多数水质指标呈正相关；康文华等[21]发现赤水河流域景观组成在河岸带尺度对水质解释率最高，建设用地和耕地具有负面影响，林地则为正面影响；de Mello等[15]研究发现萨拉普伊河(Sarapuí River)流域的土地利用/土地覆盖在子流域尺度比河岸带尺度更好地解释了整体水质的变化，森林覆盖在保持水质清洁方面发挥着重要作用，而农业和城市地区则导致水质恶化；Zhang等[16]发现三峡库区的土地利用在子流域尺度比河岸带尺度对所有水质影响更显著；Ding等[20]认为在东江流域子流域尺度相比于河岸带或点缓冲尺度更能解释总体水质变化. 综上，土地利用对水质存在显著影响，但其对水质指标解释能力最强的空间尺度存在不确定性，部分研究结果表明土地利用在缓冲区尺度下对水质解释能力更显著，而另一部分表明土地利用在子流域尺度下对水质的影响更强. 同时，土地利用特征的差异、水文过程的不同都会导致景观与水质的相关性及相关程度存在差异[22]. 因此，尺度效应成为研究重点，且多体现为子流域尺度和不同距离的缓冲区尺度之间的差异[14,22]，研究不同空间尺度土地利用与水质的定量关系，不仅有利于优化土地利用结构，还可以控制面源污染，有效进行水质管理.

黄河是我国的母亲河，也是我国水土保持、资源保护的重点流域[23]. 由于土壤、地形、降水等自然因素的影响，黄河泥沙淤积严重，加上人类活动不断加剧对河流生态的破坏，河流水质问题日渐凸显[24]. 以往研究多针对中小流域探讨土地利用与水质之间的关系[25-26]，而黄河流域跨度大，不同河段土地利用差异明显；同时，黄河流域由于其特殊的气候特征，年际径流变化明显，导致水质指标在不同季节具有显著差异. 但目前少有针对黄河流域进行土地利用与水质关系分析的研究，刘绿怡等[27]以伊河流域为研究区，基于不同距离的点缓冲区分析景观空间异质性对水质的影响. 虽然伊河流域位于黄河中下游的中心位置，但其研究只探讨了河段缓冲区尺度，未考虑到河岸带、子流域等多种尺度土地利用对水质影响的差异性. 吕振豫等[28]以2004-2015年黄河流域上、中、下游9个监测断面和4个水质指标的监测数据为研究对象，分析了黄河流域的水质达标情况及各水质指标的时空演变特征，虽然其研究对象与本研究一致，但他们只研究了水质达标情况及各水质指标的时空演变特征，并未定量探讨形成其研究结果的原因. 因此，本研究选取黄河流域作为研究区域，综合运用GIS技术、单因素方差分析(analysis of variance, ANOVA)和冗余分析(redundancy analysis, RDA)等方法，定量探讨流域不同时空尺度下土地利用指数与水质的关系，旨在找出对黄河流域水质影响最显著的时空尺度及土地利用，为优化土地利用、加强水质保护提供科学参考.

1 材料与方法

1.1 研究区概况

黄河是中国第二大河，地处31.48°～41.56°N，95.04°～117.99°E之间，流域面积达7.95×105km2，是我国生态保护、流域治理的重点对象(图1). 黄河流域幅员辽阔，地貌类型多样，由东至西海拔逐渐升高，横跨黄淮海平原、黄土高原、内蒙古高原及青藏高原4个地貌单元. 流域内草地、林地及耕地景观广泛分布，植被疏松，部分地区土壤盐碱化现象严重[29]. 黄河流域大部分地区位于干旱、半干旱地区，各地区之间气候差异显著，年平均气温在7℃左右，光照充足，湿度小，蒸发大，日夜温差明显；流域内降水量偏少，全流域多年平均降水量为451 mm[30]. 在气候和地貌的双重影响下，黄河流域水资源稀缺，年径流量仅占全国2%，且年内、年际变化大，流量呈持续减小趋势，枯水期长，水土流失严重[31]. 近年来，随着黄河流域经济的飞速发展，加速了人口城市化的进程和乡镇工业的发展，从而导致工矿企业和城镇居民生活的废污水急剧增加，且流域内耕地化肥、农药等的过量使用，致使黄河流域水质遭到严重污染. 据前人调查研究显示，1994年黄河流域干流污染物超标河长占总河长的57.4%，2003年超标河长增至78.1%，主要超标项目包括溶解氧(DO)、氨氮(NH3-N)和高锰酸盐指数(CODMn)等[28]. 从2004-2015年，黄河流域的水质达标率由48.5%增加至77.4%，整体上呈逐步改善的趋势，但黄河中游地区水质污染严重，如山西运城河津大桥监测站点和陕西渭南潼关吊桥监测站点水质达标率仅分别为0.2%和4.2%[28].

图1 黄河流域水质监测站的空间分布及子流域、河岸带、河段3个空间尺度(W1～W9代表 黄河流域内9个监测站点，依次为W1-甘肃兰州新城桥监测站点，W2-宁夏中卫新墩监测站点， W3-内蒙古乌海海勃湾监测站点，W4-内蒙古包头画匠营子监测站点，W5-山西忻州万家寨 监测站点，W6-河南济源小浪底监测站点，W7-山东济南泺口监测站点， W8-山西运城河津大桥监测站点，W9-陕西渭南潼关吊桥监测站点)Fig.1 Spatial distribution of water quality monitoring stations in the Yellow River Basin and the three spatial scales (catchment, riparian corridor and reach)

1.2 数据获取

黄河流域水质数据由国家地表水水质自动监测实时数据发布系统(http://123.127.175.45:8082/)获取，结合前人研究[28]，选取黄河上游甘肃兰州新城桥、宁夏中卫新墩、内蒙古乌海海勃湾、内蒙古包头画匠营子，中游山西忻州万家寨、山西运城河津大桥、陕西渭南潼关吊桥、河南济源小浪底和下游山东济南泺口共9个重要的具有代表性的监测站点，大部分监测点位于黄河干流附近. 于2018年2月7日-2019年2月7日每天中午12:00收集为期一年的水质数据，将此数据作为一年内每天的日均值，并分别对不同监测站点不同季节的监测数据求平均值. 前人研究结果显示，黄河流域的水体污染物未达标项目主要为DO、NH3-N、CODMn等指标[28]，故本研究的水质参数选取pH、DO、NH3-N和CODMn4个指标.

DEM数据由中国科学院资源环境科学数据中心获取(http://www.resdc.cn)，高程数据分辨率为30 m，可用于提取水系. 为了与水质数据时间对应，由地理空间数据云平台获取2018年研究区域对应的Landsat 8 OLI遥感影像数据，空间分辨率为30 m. 采取支持向量机法与人工目视解译相结合的方法，参考《土地利用现状分类》标准，得到林地、草地、耕地、建设用地、水域及未利用土地共6类土地利用类型. 经混淆矩阵精度验证，分类总精度为86.36%，Kappa系数达到0.81，符合精度要求.

1.3 不同空间尺度土地利用指数的计算

参考国内外相关研究，河段、河岸带及子流域3个空间尺度下土地利用对水质具有较好的预测效果[10,20,22]. 再结合黄河流域区域特征，黄河干流宽度多为300～500 m，且小范围内土地利用变化程度较小，故将土地利用指标分为子流域、500 m河岸带缓冲区和1000 m河段缓冲区3个空间尺度. 基于ArcGIS 10.4平台，结合DEM高程数据和中国科学院资源环境科学数据中心(http://www.resdc.cn/)发布的二级流域数据，并参考Shi等[10]划分子流域的方法，对黄河流域水系进行提取，并以监测点为出水口提取对应的子流域. 1000 m河段缓冲区以监测点上游1000 m处为圆心，设置半径为1000 m的点缓冲区；500 m河岸带缓冲区沿监测点上溯，以子流域为边界截取河流向左右两岸延伸500 m生成对应的线缓冲区. 通过对缓冲区范围和土地利用数据作相交处理，得到各空间尺度的土地利用数据，为景观指数的计算提供数据基础.

基于前人研究结果，土地利用类型百分比(PLAND)指数、斑块密度(PD)和香农多样性指数(SHDI)对水质影响较显著，且由于研究区内景观类型的差异，PD和SHDI对水质的作用结果存在不确定性[6]. 因此，本研究在景观类型水平上选用PLAND指数，用于计算各土地斑块占景观总面积的比例；在空间格局上选取了PD和SHDI指数.PD反映了景观密集程度的差异，SHDI为邻近度指标，可表征景观异质性，土地利用类型越丰富，SHDI值越高[6]. 以上景观格局指数由Fragstats v4.2.1 软件计算得到，综合体现了不同空间尺度下的个体要素和整体景观特征.

1.4 数据统计与处理

利用SPSS Statistics 19软件，首先对水质数据进行K-S检验，符合正态分布，对其进行ANOVA单因素方差分析，用于比较水质指标之间的时空差异性. 同时，对土地利用指数进行单因素方差分析，得到各指数的空间差异性.

冗余分析(redundancy analysis, RDA)作为常见的线性直接梯度排序分析方法，能够反映多个环境变量对所有物种变量的影响，同时可以得到单个环境因子对物种变化的方差贡献率[13]. 运用Canoco 5.0软件，首先对物种变量进行趋势对应分析(detrended correspondence analysis, DCA)，若分析结果得到梯度长度(gradients length)第1轴小于3.0，则选择RDA，若大于4.0，则选择典范对应分析(canonical correspondence analysis, CCA)，若值介于3.0～4.0之间，则两种分析方法均可使用[11]. 本研究中土地利用指数为环境变量，水质数据为物种变量，得到DCA结果均小于3.0，因此，选用RDA方法定量探讨黄河流域所有土地利用指数对全部水质指标的多时空尺度影响，并通过单个土地利用指数对水质变化的方差贡献率大小得到对水质影响显著的土地利用格局.

2 结果与分析

2.1 土地利用特征

以黄河流域9个监测站点为中心，提取对应的河流及子流域，得到1000 m河段缓冲区、500 m河岸带缓冲区和子流域3个空间尺度对应的土地利用类型(图1)，土地利用类型比例和景观格局指数变化如图2，表1所示. 研究区内景观分布整体以耕地和草地为主，不同空间尺度下，景观组成和空间格局存在差异. 其中，PD、SHDI、林地和水域景观空间差异性显著(P< 0.05)(图2,表1). 1000 m河段缓冲区尺度下，景观类型面积占比表现为耕地 > 草地 > 水域 > 未利用土地 > 建设用地 > 林地(图2). 500 m河岸带缓冲尺度景观类型面积占比表现为耕地 > 草地 > 水域 > 未利用土地 > 林地 > 建设用地，建设满用地面积占比大幅下降，而草地、林地面积占比上升(图2). 子流域尺度下，耕地和草地仍为主要景观，林地次之，建设用地、未利用土地和水域面积占比均较小(图2). 在多数空间尺度下，PD的值表现为1000 m河段缓冲区 > 500 m河岸带缓冲区 > 子流域，河段尺度景观斑块更加破碎，子流域尺度下曲线变化趋势相对平缓(表1).SHDI指数表现为500 m河岸带缓冲区 > 子流域 > 1000 m河段缓冲区，其中，河岸带和子流域尺度下SHDI均值较为接近，土地利用类型相对丰富，空间分布相对均衡(表1).

图2 黄河流域各空间尺度景观类型组成(平均值为在不同空间尺度下各土地利用类型总面积所占比例； 图例后*表示景观类型占比的空间性差异显著，*表示P<0.05，**表示P<0.01)Fig.2 Statistics of landscape proportions in different spatial scales in the Yellow River Basin

表1 黄河流域各空间尺度景观格局指数*

2.2 水质数据时空变化特征

流域水质参数的季节差异和空间性如图3所示. 黄河流域水质参数在时空尺度上均表现出明显的差异性. 从季节变化来看，除了CODMn，其余的水质指标均呈现出明显的季节变化(P< 0.05)(图3). 其中参数pH和NH3-N的变化趋势类似，整体表现为冬季 > 春季 > 秋季 > 夏季，且pH标准差在4个季节均相对最小，监测点数据相对稳定(图3). DO浓度变化规律为冬季 > 秋季 > 春季 > 夏季，CODMn在各季节变化不大，表现为春、夏季略高于秋、冬季(图3). 从各监测点统计值的空间变化来看，4个水质参数均呈现出明显的空间变化(P<0.05)(图3). 冬季水质指标DO浓度的空间差异最显著(P<0.05)，其次是pH在夏季空间差异显著(P< 0.05)(图3). 在多数季节尺度下，NH3-N和CODMn浓度的变化规律表现为W8(山西运城河津大桥监测站点)最大，W4(内蒙古包头画匠营子监测站点)和W6(河南济源小浪底监测站点)较小(图3).

图3 不同季节黄河流域水质参数的描述性统计(水质参数包括pH、溶解氧(DO)、 氨氮(NH3-N)和高锰酸盐指数(CODMn)；不同字母表示水质的差异显著(P<0.05)， 小写字母表示空间差异显著，大写字母表示季节差异显著)Fig.3 Descriptive statistics of water quality parameters in the Yellow River Basin in different seasons

2.3 不同时空尺度土地利用对水质的影响关系

不同时空尺度下土地利用与水质的RDA结果如表2和图4所示，结果表明，土地利用在夏季和1000 m河段缓冲区对水质变量的解释率最大，且前两轴对水质指标解释率较高，第1轴解释率至少达到61.39%，为第2轴解释率的5倍之多，说明土地利用对水质的影响关系主要取决于第1轴的解释率(表2). 对于空间尺度，在春季和冬季，土地利用对水质的空间尺度影响强度均表现为1000 m河段缓冲区 > 500 m河岸带缓冲区 > 子流域尺度，夏季和秋季则表现为1000 m河段缓冲区 > 子流域 > 500 m河岸带缓冲区尺度，但4个季节下土地利用对水质的空间尺度影响均在1000 m河段缓冲区尺度下最强(表2). 对于时间尺度，3种空间尺度下土地利用对水质影响最显著的时间尺度均为夏季(表2). 在不同的时空尺度下，对水质影响最剧烈的土地利用指数有所不同，对于1000 m河段缓冲区尺度，对水质影响最大的解释变量为耕地面积占比(30.6%～45.1%)和SHDI(23.0%～25.4%)，在500 m河岸带缓冲区尺度下解释变量主要是林地面积占比(27.4%～44.1%)和耕地面积占比(10.1%～19.4%)，耕地面积占比(11.4%～11.5%)和未利用地面积占比(10.0%～10.9%)在子流域尺度下对水质影响较大(表2). 综上，对于各个时空尺度，对水质影响最显著的土地利用指数为耕地面积占比、林地面积占比、未利用地面积占比和SHDI.

黄河流域景观指数与水质参数的RDA排序图如图4所示，RDA排序图中箭头之间的夹角表示变量之间的相关性，夹角等于零时，表示两者不相关，夹角小于90°时，表示两者为正相关，夹角大于90°则为负相关. 箭头夹角的余弦值越大，两者之间相关性越强；且箭头的长短也可以表示相关性的强弱，箭头越长，相关性越强，反之相关性越弱[11]. 不同土地利用指数与水质相关性存在差异，在多数时空尺度下，多数土地利用指数与水质指标(NH3-N和CODMn)呈负相关，如草地面积占比、水域面积占比、未利用土地面积占比和SHDI等指数；耕地面积占比和PD指数与水质指标呈正相关，分别在500 m河岸带缓冲区尺度和1000 m河段缓冲区尺度表现出较强的相关性(图4). 在多数时空尺度下，耕地面积占比、林地面积占比、未利用地面积占比与水质指标相关性最强(图4).

表2 黄河流域土地利用变量解释总体营养物质变化比例的RDA结果*

3 讨论

3.1 土地利用对水质的多时空尺度影响

本研究对于土地利用指数的选取包括土地利用类型百分比和景观空间结构两个方面，在不同的时空尺度下，不同土地利用指数对水质的影响存在明显差异，与前人研究结论相符[2,4,16]. 在多数时空尺度下，耕地面积占比与NH3-N和CODMn浓度呈现较大的正相关性(图4)，体现了“源”景观对水质的负面作用[32]. 在人为作用的干扰下，耕地景观破碎程度大，加上大量使用化肥农药的影响，污染物浓度不断增加，经地表径流等方式集中到河流中，导致水质恶化. 相较于1000 m河段缓冲区和500 m河岸带缓冲区尺度，子流域尺度上土地利用与水质的相关性较弱. 这是因为耕地面积占比在该空间尺度上相对最小，并且随着流域面积的增加，污染物迁移时间越长，降解量也就越大，对水质的恶化作用也会相应降低[12]. 在各时空尺度下，PD指数均与NH3-N和CODMn浓度呈正相关，并且在1000 m河段缓冲区和500 m河岸带缓冲区尺度下，PD与水质指标的相关性强于子流域尺度(图4). 这是因为PD值的大小可表征景观斑块的破碎程度，PD值越大，越容易造成严重的土壤侵蚀[11].

图4 黄河流域冗余分析结果排序图 (环境变量为虚线空心箭头，包括Grass-草地，Water-水域，Farmland-耕地，Urban-建设用地，Forest-林地， Unused-未利用土地，PD-斑块密度，SHDI-香农多样性指数；物种变量为实线实心箭头，水质指数包括pH， DO-溶解氧，NH3-N-氨氮和CODMn-高锰酸盐指数；(a)、(b)、(c)为春季河段、河岸带和子流域 尺度的RDA结果，(d)、(e)、(f)为夏季河段、河岸带和子流域尺度的RDA结果，(g)、(h)、(i)为秋季河段、 河岸带和子流域尺度的RDA结果，(j)、(k)、(l)为冬季河段、河岸带和子流域尺度的RDA结果)Fig.4 Ordination diagrams originated of redundancy analysis (RDA) results in the Yellow River Basin

草地在大多数时空尺度下与NH3-N和CODMn浓度呈负相关，与pH和DO浓度呈正相关，对水质有净化作用(图4)，与前人研究结论一致[2,32]. 林地在1000 m河段缓冲区尺度下，与NH3-N和CODMn的关系为负相关，与而在500 m河岸带缓冲区和子流域尺度表现为正相关(图4). 以往研究表明，作为典型的“汇”景观，林地对水化学元素具有吸收、截留的作用，从而能够保持水土、净化水质[11]. 但河段缓冲区尺度下，景观破碎程度相对较大，对污染物的吸收降解作用有所降低；且城市绿化用地距离河流较近，在种植过程中使用的化肥农药排入河流，导致水质降低[33]. 水域面积在降雨、地表径流的综合作用下，与水质呈现负相关，对水质污染有稀释作用[6]. 未利用土地由于人类活动的干扰较少，表现出与耕地相反的相关性，对河流排放污染少，对水质安全表现出有利的一面.SHDI指数在一些研究中被认为对水质存在负面作用[6,20]，但在本研究中表现为对水质的正面作用. 产生这种差异的原因可能是研究区域内主要景观类型不同，这些研究区域以林地和草地景观为主，而黄河流域内耕地占据了主导地位. 因此，SHDI值越大，景观多样性越丰富，导致林地、草地等景观的净化作用增强，而对黄河流域内的水质保护有积极作用，能够有效防止污染物扩散. 之前有研究表明随着建设用地增加，道路、房屋等城镇不透水面上的污染物会随着雨水冲刷到河流，从而造成水质恶化[10,16,20]. 而本研究结果显示建设用地在多数时空尺度下与NH3-N和CODMn浓度呈负相关(图4)，表明建筑用地面积的增加并未对流域水质产生负面的影响，原因可能是研究区内多数监测点附近的建筑用地面积较少，城镇规模较小，建筑用地面积比不足以作为影响水质的关键因子，这与王杰等的研究结果一致[4].

3.2 土地利用对水质影响的时空尺度效应

河流水质的变化与植被覆盖和地表径流量息息相关. 研究结果显示，在夏季和1000 m河段缓冲区尺度下，冗余分析得到的解释率最高，土地利用与水质的相关关系最为显著. 黄河流域以玉米、小麦等经济作物为主，主要分布在内蒙古、宁夏、陕西、甘肃、山西、青海、山东、河南等地，全面覆盖了本研究所在区域，春玉米和春小麦主要分布在黄河流域中上游，夏玉米主要分布在黄河流域中下游[34]，其生长期主要为5-9月，这一阶段农户对耕地进行大量施肥工作，污染物浓度急剧增加. 并且夏季降雨集中，降水量达到一年中的峰值，污染源通过地表径流的方式排入河流，导致夏季在各空间尺度下解释率均较大，与欧洋等[35]得到的雨季时景观对水质解释能力最强这一结论相符. 蔡宏等[22]的研究表明，相较于缓冲区尺度，子流域尺度水质对景观指数的响应更为显著. 而方娜等[14]对鄱阳湖的研究发现，缓冲区尺度土地利用对水质的解释率高于子流域尺度，与本研究1000 m河段缓冲区解释率最高一致，距离监测点越近，土地利用对水质影响越大. 产生这种差异的主要原因是在不同地区，由于土地利用空间分布的不同，水文迁移的过程存在差异. 因此，对于黄河流域，夏季和河段缓冲区尺度更能体现土地利用对水质参数的响应.

4 结论

1)在夏季和1000 m河段缓冲区尺度上，土地利用指数对水质指标的解释率最大，对黄河流域保护水资源最具有指导意义. 在夏季降雨量增多时，应有效控制农田中化肥农药的使用量，避免污染物在近河岸处的直接排放.

2)在不同的时空尺度下，不同土地利用指数对水质的影响差异明显，对水质影响最显著的土地利用指数为林地、草地和PD指数. 草地、水域、未利用土地面积占比及SHDI的增加对水质具有正面效应，能够有效保护水质安全；而耕地和建设用地、PD值的增加则会加剧水化学污染，导致流域水质恶化. 对于林地景观，需要发挥其在1000 m河段缓冲区尺度上的积极影响，通过土地利用的合理规划，增强其在河岸带和子流域尺度上的水质净化作用. 通过植树造林、退耕还林等措施提高植被覆盖率，保护现有水域水质；在进行城市土地规划时，降低各土地利用的破碎程度，减少黄河流域土壤侵蚀带来的负面影响.

3)通过对黄河流域不同时空尺度下土地利用指数与水质的定量关系展开研究，确定了对水质管理最有效的时空尺度，并反映了不同土地利用指数对水质影响的相关性. 根据研究结果可有效采取措施对黄河流域水质进行管理和优化，对保护水生态环境具有参考价值. 但本研究仍存在一些不足，由于Landsat遥感数据分辨率较低，可能对监督分类精度造成影响，在后续研究中，可采用更高分辨率影像进行解译.