北京白河流域不同空间尺度下水质对土地利用的响应关系

夏明珠，薛万来，黄俊雄，李文忠，李 垒，刘可暄，李久义，田秀君*

1.北京交通大学环境学院，北京 100044

2.北京市水科学技术研究院，流域水环境与生态技术北京市重点实验室，北京 100048

随着我国经济的快速发展和人民生活水平的不断提高，水环境的供需矛盾也逐渐凸显出来[1].有研究表明，影响水环境的因素中土地利用不仅是人类可控性最高的因子之一，也是人们从长期生产生活和科学实践中认识到造成水环境问题的最主要原因，其与流域水环境之间存在着密切的联系[2].众多学者表明，土地利用可通过多种方式对水环境质量产生影响，如刘成建等[3-4]研究表明，新修水利工程、增加建设面积、降低植被覆盖率等不同人类扰动情况下，流域水环境质量降低.同时，由于不同流域的土地利用结构存在尺度效应，因此何种空间尺度下的土地利用结构对水环境产生的影响最大，目前仍未形成统一结论[5].方娜等[6-7]的研究表明，在较小尺度(<1 000 m)的空间格局下有助于探索土地利用结构与水质的关系；而易帆等[8-9]的研究结果与之相反，其认为大空间尺度缓冲区(2 000 m)的土地利用结构与水质的相关性更加显著；另有研究如段少琼等[10-11]研究认为，子流域尺度能更好地解释水质对土地利用结构的响应关系.研究确定水质对土地利用结构响应的尺度效应是开展流域水生态保护修复的重要前提.此外，由于降雨对污染物迁移的影响，水质对土地利用结构的响应存在季节性差异.如Ly 等[12]研究表明，湄公河下游土地利用类型/土地覆盖变化与总悬浮固体和硝酸盐浓度之间的时空关系在汛期更为显著；黄益平等[13]研究表明在袁河研究区，丰水期和枯水期土地利用结构对干、支流不同水质指标的解释度各异.受土地利用结构的复杂组成以及变化的影响，水质对土地利用结构响应的季节性差异仍有待进一步研究.

由于气候变化和人类活动的影响，北京市面临着水资源短缺、地下水位下降、水土流失等水环境问题[14]，同时由于北京市人口密集、污染物种类繁多、山区地域广阔地貌复杂等因素影响，其水环境问题治理难度大.密云水库作为北京市主要的饮用水源地，上游白河流域是入库水量最大、流域面积较大、人类活动较为强烈的区域，是水环境保护的重要区域.目前已有学者对密云水库上游流域水质时空变化特征开展了研究，如Guo 等[15-17]采用景观格局探寻了不同空间尺度下土地利用结构与密云水库上游汛期白河水质的响应关系.但由于白河属于典型的季节性河流，年内的水环境质量季节性差异是其显著特征，同时因土地利用结构的空间尺度效应，精准识别流域内土地利用结构的差异性对探寻其与水质的响应关系具有重要影响.该研究采用2 期水质及“国土三调”数据(精度优于1 m)，将园地和林地作为两种不同的土地利用类型，从而直接采用土地利用结构探寻其与水质的响应关系.因此，该研究基于北京市白河流域土地利用结构，结合2022 年汛期和非汛期2 次水质实测数据，研究北京市白河流域不同空间尺度下地表水水质对土地利用结构的响应关系，以期为白河流域土地利用结构优化调整，以及密云水库上游流域乃至北京山区面向流域水生态环境质量改善的土地利用结构优化提供支撑.

1 材料与方法

1.1 研究区域概况

白河流域地处北京北部山区，属燕山山脉，北京市内涉及延庆、怀柔、密云3 个区(见图1).白河流域(40°31′N～41°03′N、116°07′E～116°53′E)，北京市内涉及河段长度126.7 km，面积2 203.9 km2，属于暖温带大陆型半湿润半干旱季风气候，年均温度9～12 ℃，年均降雨量约为620 mm[18]，土壤类型主要为棕壤和淋溶褐土，主要用地类型为林地[19].流域内水系纵横交错，有水沟道分布也较为密集，多流经中低海拔山地，平均坡度15.71°，岸坡陡峻，河身狭窄，属于典型的山区谷底河流.降雨主要集中在汛期的6-9 月，非汛期降水量较少，常发生断流，季节性差异显著.

图1 研究区域概况及采样点示意Fig.1 General situation of research area and schematic of sampling points

1.2 数据来源

小流域出水口能够反映有水沟道产汇流特点，该研究基于小流域出口选取14 个水质采样点(S1～S14)，均分布在研究区范围内，采样点信息见表1.所用水质数据分别在2022 年7 月(汛期)与10 月(非汛期)于河流中心表层(50 cm 深度)采集，用聚乙烯瓶收集1 L 水样带回实验室按照国家标准检测分析NH4+-N、TP、NO3--N、TN、CODMn浓度，并在现场使用便携式多参数水质测试仪(YSI SL1000，美国哈希公司)监测水体pH.

表1 采样点属性Table 1 Properties of sampling points

研究区域土地利用类型数据为第三次全国国土调查优于1 m 精度的数据，基本可以反映2022 年前后一段时间内的土地利用结构特征，尤其在北京山区生态保护优良的情况下，其差异可忽略不计，研究区采样点和土地利用类型分布见图2，由于研究区面积较大，选取S10 典型采样点及土地利用类型详细分布情况进行展示，可以更清晰地分辨各土地利用类型与有水沟道的关系(见图3).

图2 研究区采样点和土地利用类型分布Fig.2 Distribution of sampling sites and land use types in the study area

图3 典型采样点及土地利用类型详细分布情况Fig.3 Detailed distribution of typical sampling sites and land use types

1.3 分析方法

1.3.1 空间分析

该研究结合《第三次全国国土调查技术规范》和中国科学院资源环境数据中心的土地利用分类系统，将土地利用类型划分为耕地、林地、园地、草地、水域、建设用地、未利用地7 大类.

基于ArcGIS 10.8 软件的空间水文分析模块，将白河流域分为116 个子流域(见图2)，是白河流域研究的基础空间单元，子流域中包含14 条有水状态的河流沟道.该研究区域位于山区，坡度较大，因此选取有水沟道带状缓冲区，缓冲区半径的划分依据国内外流域水环境与土地利用结构的空间尺度效应研究结果，一般以20～2 000 m 为研究半径[20]，并且结合流域面积较小、有水沟道窄等特点，以及研究区域内土地利用类型多为林地的特征进行缓冲区的划分.因此，以小流域内有水沟道为中心，在垂直沟道方向上设置50、100、200、300、500、1 000 m 共6 个带状缓冲区尺度(见图4)，并与土地利用类型图相交得到不同空间尺度下不同土地利用类型的面积占比数据.

图4 不同空间尺度缓冲区示意Fig.4 Schematic diagram of buffer zones at different spatial scales

1.3.2 统计分析

在SPSS 26 软件中，对野外采样所获汛期和非汛期的水质指标数据进行非参数Kolmogorov-Smirnov test (K-S)检验，因水质数据不符合正态分布，故选用Spearman 相关性分析，以揭示不同尺度土地利用类型面积占比与水质指标的关系[21]，采用Mann-Witney U检验对汛期和非汛期数据的平均值进行显著差异分析.

利用Canoco 5 软件对水质指标进行趋势对应分析(detrended correspondence analysis,DCA)，结果显示水质指标数据的梯度值小于3，故选择冗余分析[22](Redundancy Analysis,RDA).将水质指标选为响应变量，土地利用类型面积占比作为解释变量，识别流域采样点环境条件和水质结果之间的相关性，揭示各土地利用类型对流域水质的贡献程度，以及土地利用类型对汛期和非汛期多个水质指标差异的解释能力.

2 结果与分析

2.1 水质时空分布特征

对14 个小流域有水沟道的主要水质指标数据进行描述性统计，结果如表2 所示.研究区域水质NH4+-N、TP、NO3--N、TN、CODMn浓度的时空分布情况如图5所示.

表2 研究区域水质指标描述性统计结果Table 2 Descriptive statistical results of water quality indicators in the study area

图5 汛期和非汛期主要水质指标的空间分布Fig.5 Spatial distribution of major water quality indicators during and after the flood season

汛期水质总体劣于非汛期.除NH4+-N、NO3--N、pH 符合GB 3838-2002《地表水环境质量标准》Ⅱ类标准外，其余指标均有不同程度超标.TN 浓度范围为0.16～7.82 mg/L，污染情况较重，在汛期和非汛期其超标率分别达100%和85.7%；NH4+-N 浓度波动范围为0.028～0.294 mg/L，汛期平均值(0.099 mg/L)低于非汛期(0.117 mg/L)；TP 在汛期和非汛期的超标率分别为7.14%和21.4%，变异系数在汛期和非汛期均属强变异类型，说明汛期和非汛期的差异较大；NO3--N 浓度的P值为0.027*，非汛期为强变异，汛期属于中等变异，汛期和非汛期NO3--N 浓度的时空分布存在显著性差异，其数值在0.02～5.54 mg/L 之间大幅波动，汛期平均值(2.01 mg/L)大于非汛期(1.01 mg/L)；CODMn浓度在汛期和非汛期存在差异，最小值和最大值分别为0.9～5.2 和0.9～4.1 mg/L；pH 属于弱变异，在汛期和非汛期其平均值差距较小，均达到GB 3838-2002 Ⅱ类标准，由此判定pH 全年数值较为稳定且受到土地利用结构的影响较小，因此pH 不作为探寻响应关系的水质指标.

2.2 土地利用特征分析

按照不同空间尺度对土地利用类型占比进行统计分析，结果如图6 所示.研究区土地利用类型以林地、园地和水域为主，且各类土地利用类型占比在不同空间尺度下有所差异.在有水沟道带状缓冲区空间类型下，林地面积占比(45.51%～78.36%)随缓冲区半径增大而增加；园地面积占比(7%～9%)较为平稳；建设用地在200 m 半径缓冲区占比(11.01%)最大，总体呈递减趋势；耕地、未利用地面积占比随缓冲区半径增大而减少，整体呈现波动趋势，耕地面积占比平均约为6.6%，未利用地面积占比平均约为1%；草地面积占比最小，由0.66%降至0.21%.在子流域空间类型下，林地面积占主导(84.87%)，其次是园地(5.45%)和建设用地(4.20%)，说明该区域植被覆盖率高，生态环境质量好.

图6 不同空间尺度缓冲区内土地利用类型面积占比Fig.6 Percentage of area of land use types in buffer zones at different spatial scales

2.3 土地利用类型对水质影响的相关性分析

汛期和非汛期土地利用类型面积占比与水质指标间的Spearman 相关性响应结果分别如图7 和图8所示.汛期未利用地与水域面积占比的负相关性在所有尺度上均较弱；建设用地面积占比与除TP 浓度外其余4 项水质指标表现出较强的正相关性，尤其是建设用地面积占比与NH4+-N 浓度在100 m、200 m、300 m 缓冲区呈显著正相关；草地面积占比与多个水质指标间相关性为负，其中NO3--N、TN、CODMn浓度与草地面积占比的负相关性更为显著，NO3--N 浓度与草地面积占比在200～1 000 m 缓冲区尺度下均呈显著负相关，虽然TP 浓度与草地面积占比在小尺度半径(50 m)下呈正相关，但随着缓冲区半径增大，二者间相关性转变为负相关，且负相关程度逐渐增强；园地面积占比与除NH4+-N 浓度外所有水质指标呈正相关，但随着缓冲区半径增大而减弱；林地面积占比与多项指标间相关性均为负；耕地面积占比与NH4+-N、NO3--N、TN、CODMn浓度在任何空间尺度下均呈正相关，与TP 浓度在有水沟道带状缓冲区时呈负相关，在子流域空间尺度呈正相关.

图7 汛期不同空间尺度下土地利用类型和水质指标的相关性分析Fig.7 Correlation analysis of landuse type sand water quality indicators at differents patial scales during the flood season

图8 非汛期不同空间尺度下土地利用类型和水质指标的相关性分析Fig.8 Correlation analysis of landuse type sand water quality indicators at differents patial scales during the non-flood season

非汛期未利用地面积占比与各水质指标间相关性较弱且正负不明确；建设用地面积占比同汛期情况相似，其与NH4+-N 浓度在200～500 m 缓冲区以及子流域尺度下正相关性均显著，但与其他水质指标间的相关性较低；草地面积占比与NH4+-N 浓度在100 m和200 m 缓冲区半径下表现为显著负相关；在各空间尺度下，园地和耕地面积占比与水质指标多呈正相关，林地面积占比与水质指标多呈负相关，但之间的相关性并不显著.

2.4 最佳缓冲区尺度识别

在相关性分析的基础上，为进一步确定影响研究区域内水质的关键环境因素，计算不同土地利用类型面积占比对水质变化的解释能力，结果如表3 所示，据此进一步筛选对水质影响最大的最佳缓冲区空间尺度.冗余分析结果表明，对于轴1 的解释率，汛期和非汛期都随缓冲区尺度的增大呈先减小后增大的波动上升趋势，汛期在1 000 m 缓冲区尺度下时达到最大值(29.43%)，非汛期在子流域尺度下时达到最大值(32.93%).

表3 土地利用类型与水质指标RDA 分析的排序轴特征值及总解释变异Table 3 Ranking axis eigenvalues and total explained variance for RDA analysis of land use type and water quality indicators

汛期土地利用类型对水质的总解释率大于非汛期.汛期1 000 m 缓冲区尺度下土地利用类型面积占比对水质指标的解释率(50.8%)最大；非汛期子流域尺度下土地利用类型面积占比对水质指标的解释率(43%)最大；从汛期和非汛期的平均解释率来看，1 000 m 缓冲区尺度下，两个时期的总解释率平均值为第一，且轴1 和轴2 对水质指标变化的解释率均较高，能够较好地反映分析结果，因此划定1 000 m 缓冲区尺度是土地利用结构对地表水水质影响的最佳缓冲区尺度.

最佳缓冲区汛期和非汛期土地利用类型与水质RDA 分析结果如图9 所示.所有的土地利用类型面积占比解释了50.8%的汛期水质差异和38.6%的非汛期水质差异.在汛期和非汛期，第一排序轴的解释率分别为29.43%和29.27%，其相关性与Spearman相关性分析结果一致.

图9 最佳缓冲区尺度下汛期和非汛期土地利用类型与水质RDA 分析结果Fig.9 Results of RDA analysis of flood and non-flood land use types and water quality at the optimal buffer scale

汛期第一排序轴所代表的水质指标NO3--N 浓度与草地、水域、林地、未利用地面积占比均呈负相关，与园地、耕地面积占比均呈正相关，与建设用地箭头夹角接近90°，因此二者间相关性极弱.非汛期第一排序轴所代表的水质指标TP 浓度与耕地、园地、草地、建设用地面积占比均呈正相关，与水域、林地、未利用地面积占比均呈现不同程度的负相关.因此，冗余分析不仅与Spearman 相关性分析相互佐证，同时也揭示了不同时期和空间尺度下白河流域土地利用结构与水质响应关系的最佳缓冲区尺度.

3 讨论

3.1 水质对土地利用类型的响应关系

土地利用主要从两方面影响河流水质，一方面是不同方式和强度的人类活动扰动，向水体中输入各种类型的污染物；另一方面是土地利用类型的变化使得地表形态发生改变，进而对地表径流和污染物迁移转化造成影响[23].该研究显示，建设用地面积占比是影响水质关键的类型，草地面积占比均与多个水质指标间均呈显著负相关，耕地面积占比均不是影响水质的主要因素，这与刘金宝等[24-26]的研究结论一致.这可能归因于山区谷底居民较少，并非以耕种为生，结合实地调研情况，耕种类型多以家庭为单位，就近种植日常生活所需的粮食、蔬菜等，耕种面积较小且连通性低，因此对水质影响小.

园地面积占比在两个时期均与大多数水质指标呈正相关，在近岸处正相关性更为显著，这与张志敏等[27]的研究结果相反.江苏省溧阳市网状河网研究区内园地面积占比小于5%，而白河流域园地面积占比较大，且大多沿有水沟道种植果树，虽然其根系对污染物能够起到积极的截留作用，但过近的距离和不合理施加化肥等行为，会带来大量氮、磷污染，因此受到人类活动的干扰明显，也存在污染物的输出效应.

3.2 水质对土地利用响应关系的空间尺度效应

冗余分析结果显示，从汛期和非汛期两个时期的平均总解释率来看，1 000 m 有水沟道带状缓冲区尺度下，土地利用结构对水质的解释率最高，这与海河流域(天津段)[28]、巢河流域[29]的研究结果不一致.有水沟道两侧的河岸带区域是陆地生态系统和水体生态系统的交界，是环境要素、生态过程、土地利用等呈现显著梯度变化的过渡带，其作为河流水环境的最后一道屏障，对污染物起着截留、净化等作用，同时不合理的土地利用结构规划，将会打破河岸带保护阈值，对水环境质量造成影响[30].海河流域天津段地处平原区，河网分布密集复杂，潮河流域地势也为低山丘陵和平原，而白河流域位于山区，地势起伏大，其建设用地多分布在谷底有水沟道两侧范围内，多为村民的住宅用地(见图3)，加之下游是北京市饮用水源地的密云水库，因此城市化程度较低，且工业产业较少，多以种植、旅游业为主.

该研究发现草地面积占比与一些污染物浓度在小尺度半径下呈正相关，在大尺度下转变为负相关，这表明在近岸带草地对径流污染物的迁移截留作用较为微弱，可能是因为近岸带草地被居民生活生产破坏，而远岸带草地生态效益较好，且污染物经过较长距离的草地截留，有效降低了污染物含量，因此在1 000 m 缓冲带尺度下负相关性更强.另外，有的污染物浓度在有水沟道带状缓冲区与子流域尺度下，与土地利用类型面积占比的相关性相反，可能是因为带状缓冲区尺度下各土地利用类型面积占比较子流域尺度更为均匀，子流域尺度下林地面积占比大于所有半径带状缓冲区尺度(见图6)，较大的林地面积占比改变了其相关关系.

3.3 水质对土地利用类型响应关系的季节性差异

多个流域的研究表明水质对土地利用结构的响应关系存在季节性差异，如鄱阳湖流域[31]、黄河流域[2].该研究中汛期土地利用结构对水质解释率均高于同空间类型尺度下的非汛期(见表4)，这与Johor River流域[32]、赣江[33]等地区的研究结果相同.降雨在影响流域水质季节性差异中具有两面性，高强度的降雨使得径流量增加，水体中的污染物得以稀释，水质变好；强降雨也会增强地表冲刷力，采样点上游的污染物随水流冲刷，汇流输入水体，污染水环境.汛期污染物浓度较高，表明白河流域水质污染情况主要受到地表冲刷影响，白河流域汛期降雨量普遍高于非汛期，高降雨强度与多土壤含水量有利于地表径流的形成和养分的输出.白河流域汛期降水量占全年的66%，而非汛期降水较少，汇流过程不显著.

该研究表明水质对土地利用结构的响应关系具有空间尺度效应和季节性差异，由于研究区域和空间数据差异所得结果不同.因此，在后续研究中将扩大研究区范围从整体考量，并且使用多年的水质数据，开展更适宜北京山区的研究，进一步突出其在时间跨度和空间尺度的差异.根据密云水库上游流域水质指标控制要求需达到GB 3838-2002 Ⅱ类标准，建议白河流域在乡村建设过程中，尽可能远离河道，严格高效处理生活垃圾，加强工业废水和生活污水的排放整治，同时增强草地、林地等土地利用类型的面积占比；同时在种植果园过程中做到合理施肥和水土保持，从而减少化肥带来的氮、磷污染对水环境造成破坏.

4 结论

a) 北京山区白河流域土地利用类型在有水沟道带状缓冲区和子流域空间尺度下，均以林地、园地、建设用地为主要用地类型，林地面积占比在半径为1 000 m 的带状缓冲区达最大，为78.36%，在子流域为84.87%，为最主要的土地利用类型.

b) 水质指标在汛期和非汛期两个不同时期存在差异.汛期水质整体劣于非汛期，NH4+-N、TP、TN 浓度汛期较低，而NO3--N、CODMn浓度非汛期较低.与GB 3838-2002 Ⅱ类标准限值相比，TP、TN 和CODMn均存在不同程度的超标，TN 超标尤为严重.TP、NO3--N属强变异类型，空间差异显著.

c) 土地利用类型对河流水质有显著影响，汛期影响更明显.建设用地和园地面积占比与各水质指标均呈正相关，与NH4+-N 浓度的相关性尤为显著，草地、林地、水域面积占比与水质指标整体呈负相关，其中草地面积占比与多个水质指标的负相关性显著.

d) 冗余分析结果表明，汛期土地利用结构对水质解释率高于非汛期.汛期各空间尺度的总解释率相似，1 000 m 有水沟道带状缓冲区尺度对汛期水质变化有更好的解释能力；子流域尺度下能更好地反映土地利用结构与水质在非汛期的响应关系；两个时期综合来看，1 000 m 缓冲区尺度是白河流域水质对土地利用结构响应的最佳缓冲区尺度.