典型网状河网区域土地利用和景观格局对地表季节水质的影响<br/>——以江苏省溧阳市为例*

典型网状河网区域土地利用和景观格局对地表季节水质的影响
——以江苏省溧阳市为例*

2022-09-05张志敏杜景龙陈德超

湖泊科学 2022年5期

张志敏，杜景龙，陈德超，张飞

(1：苏州科技大学地理科学与测绘工程学院，苏州 215009) (2：新疆大学资源与环境科学学院，乌鲁木齐 830046)

水质在生态安全、人类健康、经济发展和休闲娱乐方面发挥着关键作用[1-3]. 快速城市化和人口增长给世界各地的水生态系统带来了巨大的压力，使水污染成为一个全球性的重要问题[4-5]. 流域水质不仅受气候降雨、土壤侵蚀和水文循环等自然过程的影响，而且还受人类经济社会系统的干扰[6]. 作为自然和人类干扰综合体现的土地利用和景观格局深刻影响着区域水文过程和水质[7-9]. 河流水质与周边土地景观特征之间的关系往往受到时空变化的影响[2,10-11]. 因此，深入理解景观特征与河流水质之间复杂关系已成为当务之急.

景观特征(包括景观组成和景观格局)对水生态系统的功能具有重要影响[12-15]. 河流缓冲区是一个营养保留区，能够通过沉积/捕获、植物吸收和微生物反硝化作用极大地减少径流中的沉积物、营养物质、酸性和有毒化学物质[2,16-20]，缓冲区的生态效应因距离水体距离远近而不同. 因此，水质指标和景观特征之间的关系具有尺度敏感性[17,20]. 之前许多研究探索了不同尺度缓冲区内的土地利用/景观格局对水质的影响[2,17,19,21]，其影响方向和程度因具体区域特征而不同. 因此，在资金有限的情况下，探索区域景观因素对流域水质的最有效影响尺度对合理的景观规划设计具有重要的指导价值. 研究还表明，流域水质存在显著季节性差异，因为季节性降水和地表径流对河流的流量和污染物浓度有显著的影响[1,22-23]. Li等发现旱季河岸带土地利用对汉江上游流域水质影响最大[24]，Zhang等和Wu等也发现土地利用和景观格局对旱季水质的影响更为显著[3,25]. 然而，Shen等发现雨季景观特征与流域水质存在更强的相关性[17]. 探讨流域景观特征与河流污染物关系的季节性变化，对于优化区域景观格局和改善流域水质管理具有重要意义[5].

为了探索自然过程和人类活动与水质之间的关系，以往研究中应用了很多建模工具和多元统计技术. Zhang等通过多元线性回归模型拟合不同尺度下景观指数与水质指标的关系[26]，但这种拟合技术不能排除解释变量之间的多重共线性. Zhao等采用结构方程模型分析了不同土地利用类型对水环境质量的贡献率[27]，但该方法无法定量预测水质变化. SWAT(soil and water assessment tool)分布式流域水文模型和其他基于物理的建模工具也被用于水质预测，但这些建模工具所需参数过多，较难在缺资料地区推广使用. 与上述方法相比，偏最小二乘回归(PLSR)技术可以用较少的成分有效地解释响应变量的方差，且当存在多重共线性以及预测变量的个数大于观测样本数量时，PLSR是最有效的预测方法之一[22,28-30]. Wang等使用PLSR分析了亚热带农业流域影响河流磷浓度和负荷变化的主要因素，并证明了使用PLSR识别重要因子和构建预测模型的可行性[31].

流域的生态环境特征(包括土地覆被、地形地貌、地质土壤和气候水文等)和水质之间的关系是复杂的，而且具有区域特色[14,32-33]. 虽然国内外学者对区域景观特征对水质的影响进行了大量的研究，但研究区域通常是由边界清晰的树枝状河流组成的水文单元[34]，关于网状河网区域的研究较少. 由于低洼平坦的地形加上丰沛的降雨，网状河网区域形成了密集的河道和复杂的网状水系，河流流速缓慢. 在人工运河及水闸控制下，河流流向不定，因此很难界定河流上下游关系和子流域边界. 上述河流水动力学状况下的土地利用和景观格局对水质的影响可能与其他区域存在很大不同[35-36]. 溧阳市作为长江三角洲典型网状河网区域，较高的人口密度和发达的工农业生产导致河网的水质对当地土地利用模式的敏感性高于中国其他地区. 本研究旨在：1)分析典型网状河网区域多缓冲区尺度景观特征的空间变异性；2)通过主成分分析(PCA)选取主要水质因子，采用冗余分析(RDA)确定区域景观因子对水质指标的最佳影响尺度；3)应用PLSR方法确定影响季节水质的最重要的土地利用和景观特征指标.

1 材料与方法

1.1 研究区域

江苏省溧阳市位于太湖流域西部，总面积1535.87 km2(图1). 属东亚季风气候区，年平均气温为16℃，年平均降水量为1147 mm. 市域内西、北、南三面环山，中部是平原圩区，河网众多，湖泊星罗棋布. 溧阳市作为长江三角洲的典型网状河网区域，河道密集，网状结构复杂，主要大中型水库有沙河水库和大溪水库，主要湖荡有长荡湖和前马荡，主要河流有丹金溧漕河、中河、北河和南河等. 土壤类型主要为水稻土、黄棕壤和黄褐土，主要农作物有水稻、油菜、茶叶和蚕桑等. 2017年年末，全市户籍人口为79.1万人，其中城镇人口为45.91万人，城市化率为60.2%，GDP达858.04亿元，人均GDP为11.26万元. 全市共有工业企业5000多家，其中规模以上工业企业402家(2018溧阳统计年鉴). 区域主要污染源包括工业、城市污水、农业和畜禽养殖业等.

图1 研究区域监测断面分布Fig.1 Distribution of monitoring sections in the study area

图2 溧阳市月平均降雨量分布Fig.2 Distribution diagram of average monthly rainfall in Liyang City

1.2 监测断面及水质指标

1.3 土地利用和景观指数

利用美国陆地卫星Landsat-8遥感影像(30 m分辨率)解译获得区域的土地利用数据. 一般研究区一年之内各类土地利用的分布和面积变化不大，因此从地理空间数据云(http://www.gscloud.cn/)下载了影像质量较好的2017年5月的Landsat-8 OLI影像数据(云量为0.23%). 在ENVI 5.2软件中对遥感影像进行了大气校正、几何校正和图像增强等预处理，然后利用具有最大似然算法的监督分类技术对土地利用进行分类，并根据目视判读和现场调查对分类进行了修订. 采用混淆矩阵评价土地利用分类的准确性，分类的总体准确度为85.6%，Kappa系数为0.84. 土地利用类型分为6类：耕地、园地、建设用地、林草地、水域和其他土地(图3).

网状河网区域的河流流向不定且无法确定明显的上下游关系，且监测断面水质受多方来水影响，因此宜采用以监测断面为中心的圆形缓冲区作为分析单元[37]. 缓冲区半径主要依据国内外区域水环境与土地利用/景观格局关系的尺度效应研究来确定，大多以100～500 m作为最小空间尺度，最大空间尺度则一般为3000～5000 m[2-3,14,21,26,32,34]. 根据研究区范围和河岸土地利用类型，以监测断面为中心划定了500、1000、1500、2000、2500和3000 m 6种缓冲区(图3). 然后计算各缓冲区范围内的土地利用类型面积占比和景观水平的所有景观指数.

图3 研究区域土地利用类型分布Fig.3 Distribution of land use types in the study area

为了研究景观格局对生态过程的影响，景观生态学家和其他研究人员开发了许多景观指数[38]. 然而，由于这些指标之间的多重共线性和跨尺度的不稳定性，并非所有景观指数都能提供描述景观配置模式的有用信息[16,39-40]. 通过相关分析保留相对独立且对景观格局分析比较重要的指标，最后选取的景观指数包括：斑块密度(patch density，PD)、最大斑块指数(largest patch index，LPI)、边缘密度(edge density，ED)、景观形状指数(landscape shape index，LSI)、蔓延度指数(contagion index，CONTAG)、散布与并列指数(interspersion and juxtaposition index，IJI)、斑块丰富密度(patch richness density，PRD)、香农多样性指数(Shannon’s diversity index，SHDI)和香农均匀度指数(Shannon’s evenness index，SHEI). 在以前的研究中，它们通常被用来解释水质的土地利用模式[29]. 此外，上述景观指数对于理解景观中的生态功能和人类感知非常重要[41-42]. 景观指数的计算通过Fragstats 4.2软件完成.

1.4 统计分析

主成分分析(PCA)是最重要的降维方法之一，可以有效提取重要指标和消除数据噪音. 本文采用PCA选取主要的水质指标因子[14,43-44]. 对水质指标数据进行除趋势对应分析(DCA)后发现，数据集的梯度值为0.3(小于4)，故选择冗余分析(RDA)来确定最佳缓冲区尺度[14,44]. 通过SPSS 20进行多变量统计分析，使用Canoco 5.0软件执行RDA.

采用偏最小二乘回归(partial least squares regression，PLSR)分析景观因子与水质指标的关系，并确定水质退化的关键预测因子[29]. 在SIMCA 15软件中进行PLSR分析[45-46]，并使用交叉验证(cross-validation)确定最优PLSR模型所需的最小潜在成分数.Q2(一个成分预测因变量总变异的分数)代表成分的交叉验证，当Q2大于0.5时，该模型有望具有良好的预测能力. 为了避免过度拟合问题，通过迭代移除非显著变量来优化每个PLSR模型，以最小化响应中解释方差(R2)和模型预测能力(Q2)之间的差值[47]. 每个解释变量在模型拟合中的贡献由变量投影重要性(variable importance in projection，VIP)确定. 一般认为，VIP值大于1的解释变量对响应变量有重大意义[45]. 回归系数(regression coefficient，RC)表示PLSR模型中每个解释变量对响应变量影响的方向和强度. PLSR在分析土地利用类型、景观指数与水质指标相关性和预测方面有广泛的应用[22,29-31].

2 结果与分析

2.1 多尺度缓冲区土地利用和景观格局特征分析

不同尺度缓冲区内土地利用面积组成的统计结果见图4. 研究区土地利用类型以耕地、建设用地和水域为主. 随着缓冲区尺度的增加，耕地、林草地和园地面积占比平均值呈增大趋势，建设用地和水域面积占比平均值呈减小趋势.

图4 不同尺度缓冲区内土地利用类型面积占比Fig.4 Area percentage of land use types in different scale buffer zones

监测断面S7和S10周边以耕地为主，所有尺度缓冲区内耕地面积占比均在50%以上，S9次之，耕地面积占比在40%以上，S4的耕地面积占比由1000 m缓冲区的24%增加到3000 m缓冲区的49%；监测断面S2、S3和S6周边以建设用地为主，建设用地面积占比在所有尺度缓冲区内基本都在60%以上，S1、S4、S5和S7的建设用地面积占比随缓冲区尺度的增加而迅速减小；S11的林草地面积占比在25%～30%之间，远远大于其他监测断面(基本在5%以下)；S11和S12的园地面积占比大于其他监测断面，且随缓冲区尺度的增加而增大，其中S11和S12的园地面积占比分别为15%～25%和8%～16%，其他监测断面园地面积占比基本在5%以下；大多数监测断面的水域面积占比在20%以上，尤其在500～1500 m缓冲区内S12的水域面积占比在85%以上，在500和1000 m缓冲区内S11的水域面积占比在50%以上，S1、S5、S9和S10的水域面积占比基本在30%以上. 考虑到各尺度缓冲区内其他用地的占比都很小，在后面的分析中排除了其他用地.

不同尺度缓冲区内景观指数统计结果如图5所示.LPI是景观优势度指标，反映人类活动的方向和强度.LPI的均值随着缓冲区尺度增加而逐渐减小，最小值出现在 3000 m 缓冲区内，说明随着空间尺度的增加，景观优势度逐渐减弱，人类活动的干扰作用增强.ED反映景观中异质性斑块之间物质和能量交换的潜力及相互影响的强度[14]，ED的最小值、中值、第一四分位数(Q1)均随着缓冲区尺度增加而逐渐增大，且ED值分布趋于集中.PD、CONTAG和IJI指标均反映景观破碎化程度.PD的最大值、第三四分位数(Q3)均随着缓冲区尺度增加而逐渐减小.CONTAG值随着缓冲区尺度增加表现出先增大后减小的趋势，在1500 m缓冲区达到最大值，说明景观中的优势斑块在1500 m尺度上具有良好的连通性.IJI对受到某种自然条件严重制约的生态系统的分布特征反映显著.IJI在各缓冲区尺度随机分布，其散布分布没有明显规律，说明人类活动代替了自然条件的制约作用对区域景观的影响更加复杂.LSI反映了斑块形状的复杂程度.LSI值从500 m到3000 m 缓冲区逐渐增大，表明在较大尺度缓冲区内斑块形状的复杂程度和人类活动的干扰程度最强.SHDI、SHEI和PRD等多样性指数反映景观异质性和均匀度. 在所有缓冲区内SHDI和SHEI没有显示出明显的规律，而PRD的最大值出现在500 m缓冲区内，随着缓冲区尺度增加PRD的值明显下降. 可能是随着空间尺度增加，斑块类型数目没有减少，而区域面积的增大使得PRD显著减小.

图5 不同尺度缓冲区内景观指数箱线图Fig.5 Boxplots of landscape metrics in different scale buffer zones

表1 水质指标主成分的载荷值、特征值及解释方差*

2.2 主要水质指标因子的选取

考虑到沙河水库和大溪水库的集水区为典型丘陵流域，相应的S11和S12监测断面不能归为网状河网区，并且按照地方饮用水保护等相关法规，取水口水质受到法律保护，缓冲区景观对监测断面水质的影响有限，特别是S11和S12北部的缓冲区域，实际上不在两大水库流域范围内，对这些监测断面水质基本无影响，因此在后面的分析中剔除了S11和S12监测断面.

图6 主要水质指标的空间分布Fig.6 Spatial distribution of main water quality indicators

2.3 最佳缓冲区尺度的识别

为了探索土地利用和景观指数对区域水质的最有效影响尺度，对10个监测断面(剔除了S11和S12监测断面)的主要水质指标与6个缓冲区内的土地利用和景观指数进行了RDA分析. 分析结果显示(表2)，所有缓冲区内土地利用因子的总解释变异都在70%以上，景观指数的总解释变异都在77%以上. 与土地利用相比，景观指数更加有效地解释了水质的变化. 土地利用与景观指数随缓冲区尺度变化表现出一致性趋势. 在500 m缓冲区内，土地利用与景观指数的总解释变异分别为78.7%和90.1%；当尺度增加到1500 m时，总解释变异减小至72.8%和77.7%；当尺度增加到2500 m时，总解释变异逐渐增大到79.0%和97.8%，达到了最大值；而在3000 m 缓冲区尺度，总解释变异又开始减小. 总之，随缓冲区尺度的增加，土地利用与景观指数对水质的总解释能力呈先减小后增大最后又减小的趋势，并在2500 m缓冲区尺度达到最大. 因此，可以确定2500 m是土地利用和景观指数对区域水质影响的最佳缓冲区尺度.

上述结果表明，研究区土地利用和景观格局特征在较大的空间尺度上对水质影响更大[3,21]. 城市工业废水和市政污水排放等点源污染往往在较小的河岸带尺度对水质影响较大，而农业面源污染则在较大区域尺度上影响水质[2]，说明区域点源污染得到了一定的有效控制，而面源污染成为区域的主要污染源，因此应重点加强农业管理实践(如保护性耕作及合理施肥等)和畜禽养殖污染源管理. 同时，如图4和图5所示，随着缓冲区尺度增加，耕地、林草地和园地面积占比呈增大趋势，建设用地和水域面积占比呈减小趋势，LPI和CONTAG值下降，ED和LSI值上升，进一步说明了2500 m缓冲区尺度内农业面源污染更加严重，并且人类活动加剧了景观的破碎化和复杂程度. 因此需要在2500 m等较大空间尺度上优化耕地、林草地和园地的聚集分布格局[21,40]，形成优势斑块和生态源地，改善区域河流的连通性，减少对景观结构的人为干扰，从而改善区域的水文状况和生态功能.

表2 土地利用/景观指数与水质指标RDA分析的排序轴特征值及总解释变异

2.4 2500 m缓冲区景观因子对季节水质的影响

表3 雨季和旱季水质指标的最优PLSR模型汇总

RC和VIP是反映土地利用和景观指数对特定水质指标重要性的综合表达方式. 表4列出了雨季和旱季水质指标的各最优PLSR模型的回归系数(RCs)及其关键变量(VIP>1).

建设用地、PD、LPI和ED在雨季和旱季都没有表现出对任何水质指标有重要影响. 值得注意的是，最优模型中的所有解释因子与特定的水质指标都有一定的相关性，而只有VIP值大于1的解释因子才被认为是最重要的. 总体而言，最优模型的VIP值最高的关键变量在雨季和旱季不同，景观指数在雨季的影响更大，而在旱季土地利用和景观指数重要性相当.

3 讨论

3.1 对水质影响最显著的缓冲区尺度

RDA分析结果表明，2500 m缓冲区内土地利用和景观指数对水质方差的解释能力最大，解释比例分别达到79.0%和97.8%. 从图4得知，随着缓冲区尺度的增加，耕地、林草地和园地面积占比呈增大趋势，而耕地、林草地和园地是影响区域水质最显著的土地利用类型；从图5得知，影响区域水质最显著的景观指数CONTAG随着缓冲区尺度增加表现出先增大后减小的趋势，在2000 m以后趋于稳定. 综上所述，土地利用和景观指数在2500 m等较大的缓冲区尺度对水质的影响较强. 然而不同研究确定的最佳缓冲区尺度差别很大. 王小平等研究确定4 km缓冲区是艾比湖流域水质预测的最佳尺度[14]，李艳利等分析确定浑太河上游流域300 m缓冲区景观因子对水质具有最大的解释能力[44]. 产生不同结果的原因在于各个流域独特的景观特征，可能与区域气候条件、地形地质、水文过程、数据精度和缓冲区选择的差异有关. 另外，污染源的分布、人口密度、人均GDP和产业结构差异等因素也可能影响水质空间分异[5,31].

3.2 影响水质的关键土地利用类型

不同的土地利用类型通过影响水文条件、地表径流、污染源分布及污染物迁移转化等流域特征最终影响水质状况[48]. 通过对雨季和旱季水质指标与景观因子的PLSR分析得知，耕地、林草地和园地是影响区域水质最显著的土地利用类型. 其中，耕地与雨季的NH3-N和旱季的Petrol、COD呈显著的正相关，林草地与大多数水质指标呈负相关，这与之前的研究一致[2,21]. 园地与除pH、DO和S2-外的其他雨季水质指标均表现出很强的负相关，这与某些研究不一致[10,49-50]，但吉冬青等对流溪河流域的研究也发现园地与CODCr、TP和NH3-N呈负相关[51]. 如图4所示，大多数监测断面缓冲区范围内的园地面积所占比例都很小(小于5%)，仅在S11和S12比例较高. 本研究排除S11和S12后重新进行PLSR分析得知，园地仍然与大多数雨季水质指标呈负相关，只是显著性有所降低. 据此可以推断，S11和S12监测断面缓冲区范围内较大的园地面积比例对结果产生了放大效应. 同时，园地植被对污染物有截留、吸收和吸附的作用，在严格控制园地的化肥农药使用量和推广使用有机肥、生物农药的措施下，园地在一定程度上能对水质起到有利影响[36,52]. 然而，建设用地对水质的影响不显著，可能归因于网状河网区域河道纵横，建设用地被河网分割呈破碎化分布，建设用地对水质的不利影响被极大的削弱，也可能与监测网点的分布和监测数据有关.

3.3 景观格局对网状河网区域水质的影响

PLSR分析结果表明，CONTAG、IJI、SHDI和SHEI是影响区域水质最重要的景观指数.CONTAG与大多数水质指标呈显著正相关，先前的研究同样发现CONTAG与NH3-N、TP和TN呈显著正相关[37,53-54].CONTAG可以反映不同景观类型的团聚程度或蔓延趋势，高CONTAG值表明景观中的优势景观类型(耕地和建设用地)聚集蔓延导致面源污染物集中产生与输出，从而加剧了水质恶化[54]. 蔡宏等和吉冬青等则发现CONTAG与水质指标呈显著负相关[15,51]. 不同结论的产生可能与研究区主导土地利用类型有关.

IJI、SHDI和SHEI与大多数水质指标呈负相关，与同样是典型网状河网区域的苏州市和常州市的研究结果一致[37,54-55]，但与其他一些区域的研究结果正好相反[15,29,32,51]. 究其原因，IJI、SHDI和SHEI反映了整个景观内斑块的邻接均衡和多样化程度，IJI、SHDI和SHEI的值越大，各斑块类型的比邻概率越均等，在面积和数量上的分布越均匀，建设用地和耕地等主要污染来源类型在景观中的集中程度和主导作用也会随之减弱[54]，同时考虑到网状河网区域的大面积水域和多方来水稀释作用，因此景观多样性和均衡化反而缓解了对水质的不利影响.

与以往研究结果相比[32,51,53-54,56]，网状河网区域的优势景观类型和河网水文特性决定了其与具有不同景观特征流域的分析结果相差较大. 景观格局特征与水质的关系具有不确定性，但区域景观特征对水质具有显著的影响已得到普遍认同[54]. 因此，不同景观特征的区域需要继续深入研究，有针对性地合理优化景观格局，以有效治理污染源和改善水环境状况.

3.4 土地利用/景观格局与水质关系的季节性差异

大量研究表明土地利用/景观格局与河流水质的关系随季节变化而不同[1-3,8,36]. 本研究中旱季大多数水质指标PLSR模型的R2和Q2以及提取的成分数都大于雨季，因此旱季PLSR模型的显著性和预测能力比雨季强，但景观指数在雨季的影响更大[3,25]. 这是由于网状河网区域地势平坦、水网纵横且人类活动强度大，而人类活动是旱季的主导干扰因素，土地利用和景观格局作为人类活动的综合反映对旱季河流水质的影响更为显著[48]. 因此优化非点源污染治理时尤其应加强旱季人类活动对水环境的干扰管控.

需要说明的是，本研究仅获得了2017年单月的水质监测数据，对于分析土地利用/景观格局与水质的季节性关系，虽略显数据不足，但雨季和旱季均有3个月的平均观测数据，也能进行相关的统计分析并得出有效的结果. 若能获取12个月的水质监测数据，研究结论会更加有说服力. 另外，将雨季和旱季的水质分别与5月份的土地利用/景观格局进行分析虽显得不够严谨，但也存在类似的研究，如Zhang等研究了三峡库区水系2015年土地利用对2015-2016年期间干季和湿季水质的影响[3]，Bu等研究了太子河流域2007年土地利用格局与2009年旱季和雨季水质的关系[10]. 如果能采集到多年度不同时间点的水质数据，后续将进一步开展多时间序列土地利用/景观格局变化对水质的影响机理及驱动机制的耦合研究.