北运河北京段水生态承载力研究及关键控制要素识别

杜 芳, 彭剑峰, 王亚俊, 张 旺

1.中国矿业大学， 江苏 徐州 221116

2.清华大学环境学院水质与水生态研究中心， 北京 100084

3.中国矿业大学(北京)， 北京 100083

北京市副中心承担着疏解非首都功能、推进京津冀协同发展的功能，区域内主要水系属于北运河水系[1]. 北运河北京段承载了北京市主要人口，多年来水质状况较差，水生态现状不容乐观[2-3]. 随着北京市副中心的快速建设和发展，区域内水资源-水环境-水生态面临严峻的承载压力. 生态系统作为城市生态系统的承载基底，水生态承载力制约着城市发展的速度和规模. 如何解析北运河北京段水生态承载力演变过程，并识别关键控制要素，是提升副中心地区生态承载力的前提要求[4-5].

水生态承载力(WECC)是近年新兴的，基于承载力理论，综合水资源、水环境等要素对水生态系统进行研究的一个综合性的概念，目前学术界对其定义及研究方法仍在不断探索中，尚未达成统一共识[6]. 水生态承载力指标体系应综合水生态系统的水资源、水环境等要素. 国内外学者多采用供需平衡法、生态足迹法、综合指标法等方法评估流域水生态承载状况[7-8]，其中，综合指标评价法可综合考虑各类因素，结构体系清晰，便于厘清各因素间关系及确定关键控制要素[9]. 水生态系统是城市生态系统的重要组成部分，同时水生态承载力制约着城市发展的速度和规模. 北京市副中心所属北运河北京段，承担着北京市排水、灌溉等多项功能，当前关于水生态承载力的相关研究主要集中于水资源承载力、水环境容量以及水生生物结构功能完整性等内容[10-11]. 现有研究中仍缺乏综合水资源、水环境、水生态指标长时间序列下的水生态承载力评估和关键要素识别[12-13].

该研究通过综合指标法构建北运河北京段的水生态承载力评估体系，分别计算了1995年、2000年、2005年、2010年、2015年和2018年北运河北京段区域整体和子流域的水生态承载力，分析了近30年水生态承载力演变特征，并采用多层感知器神经网络模型识别水生态承载力的关键控制要素，以期为副中心生态承载力提升和压力减负提供科学参考.

1 材料与方法

1.1 研究区概况

北运河发源于北京市军都山南麓，自西北向东南流入通州区，后流经河北省香河县、天津市武清区，并在天津市大红桥汇入海河. 北运河上游以山区丘陵为主，中下游为冲积平原，北运河北京段区域面积约 4 348 km2，主河道长89.4 km. 北运河流域属于暖温带大陆性季风气候，降水集中于6—9月，多年平均降水量为581.7 mm，多年平均天然径流量为4.8×108m3，年内来水量变化明显，非汛期以再生水排水为主. 北运河北京段包括沙河、蔺沟河、清河、通惠河、温榆河、凉水河、坝河、小中河、北运河、凤河10个子流域(见图1)，承载了北京市70%以上的人口，该区域经济总量占北京市的80%以上，容纳了北京市90%以上的污水排放量. 根据《北京市水资源公报(2019年)》，2019年北运河水质以GB 3838—2002《地表水环境质量标准》Ⅳ类和Ⅴ类为主，水环境质量较差. 此外，就流域水生态状况而言，京津冀对北运河流域水资源利用率较高，流域季节性断流频发，河流自净能力差，生态基流难以保障.

图1 北运河北京段子流域示意

1.2 研究方法

1.2.1水生态承载力评估体系构建及权重确定

该研究从水资源、水环境、水生态3个层面收集了30多个指标，并进行了因子有效性分析[14]. 在SPSS 19.0软件中对指标进行主成分分析并进行最大方差法因子旋转，指标体系通过了KMO度量、Bartlett球形度检验(KMO>0.6，Sig.=0.000). 在此基础上划分出水资源禀赋、水资源利用、水环境纳污、水环境净化、水生生境、水生生物6个分指标，人均水资源量、年降雨量、水资源开发利用率、农业污染强度指数、工业和生活污染强度指数、水环境质量指数、河岸带林草覆盖率、生境连通性指数、生境整体性指数、生境聚集度指数和叶绿素含量11个子指标(见表1)，构建了北运河北京段的水生态承载力评估体系.

为避免人为赋值带来的主观性影响，对各项指标进行标准化处理后采用熵权法确定指标权重[15-16]. 熵权法计算权重过程中，如果指标信息熵越大，则该指标提供的信息越充分，权重赋值越大. 正向指标和负向指标标准化的计算分别如式(1)(2)所示，第i个样点第j项指标权重的计算公式如式(3)～(5)所示.

(1)

(2)

(3)

(4)

(5)

通过指标标准化和权重的确定，北运河北京段水生态承载力评估体系及指标权重如表1所示.

表1 北运河北京段水生态承载力评估体系

1.2.2水生态承载力分值计算

水生态承载力指标体系中水资源、水环境和水生态等3个总指标的计算如式(6)所示. 水生态承载状态分值T值由总指标加权获得，计算公式如式(7)所示.

(6)

T=A×WA+B×WB+C×WC

(7)

式中，Z为水生态承载力承载总指标，Ij为分指标第j个子指标指标值，Wj为第j个子指标权重，Wk为该总指标第k个分指标权重，WA、WB、WC分别为水资源、水环境、水生态总指标，A、B、C分别为总指标WA、WB、WC的权重.

根据水生态承载力状态分值进行承载力状态等级划分，结果如表2所示.

表2 水生态承载力等级划分

1.2.3水生态承载力差异性识别

单因素方差分析仅考虑一个因素的影响而固定其他因素，也可用于检验样本量为两个以上时均数差别是否具有显著性[17-18]. 为了对子流域间水生态承载力的差异性进行识别，首先应对各子流域近30年的水生态承载力综合分值进行方差齐性检验和正态分布检验(Shapiro-Wilk检验和K-S检验). 在方差齐性检验和正态分布检验均通过后进行单因素方差分析，这一过程在R 4.0软件中运用multcomp包完成.

1.2.4区域承载力关键控制要素识别方法

为识别承载力关键控制要素，在IBM SPSS 19.0软件中对各项指标进行正态分布检验，检验方法为K-S(Kolmogorov-Smirnov)法，对满足正态分布和不满足正态分布的指标分别用Pearson相关系数和Spearman相关系数进行分析[19]. 然后在IBM SPSS 19.0软件的神经网络模块中按照训练集和测试集7∶3 比例建模，对各子流域1995年、2000年、2005年、2010年、2015年和2018年共6期的水生态承载力综合指标和子指标进行多层感知器神经网络模型分析，并通过标准化后的重要性分析得出区域承载力关键控制要素[20]. 综合分析相关系数和神经网络模型重要性，最终得出区域承载力关键控制要素.

1.3 数据来源

1.3.1指标数据来源

北运河北京段区域1995年、2000年、2005年、2010年、2015年和2018年6期各子指标量化方法和数据来源如表3所示. 其中，水资源相关数据来自1995—2018年《北京市统计年鉴》《北京市通州区统计年鉴》和《北京市水资源公报》. Landsat TM/ETM+/OLI遥感影像(分辨率为30 m)、土地利用类型数据、相应年份的年降雨量、人口、GDP等空间栅格数据(分辨率为1 km)来自中国科学院资源环境科学与数据中心(http://www.resdc.cn/Default.aspx)[21-22]. 水生生境指数中的河岸带林草覆盖率由ArcGIS(500 m)缓冲区分析得到，生境连通性指数、生境整体性指数、生境聚集度指数分别由Frgstats 4.2软件中的CONNECT指标、COHESION指标和AI指标计算得到. 生境连通性指数用来表示景观间的连接状况，景观间连接状况越好则生境连通性越强；生境整体性指数表示景观斑块间的整体性，景观间整体性越强则生境整体性指数值越高；生境聚集度指数表示景观聚集程度，即景观中不同斑块类型相邻的概率；生境聚集度指数随着景观聚集程度的增加而增大[23-24].

表3 北运河北京段水生态承载力数据来源及量化方法

生境连通性指数(CONNECT)、生境整体性指数(COHESION，%)和生境聚集度指数(AI，%)的计算公式分别如式(8)(9)(10)所示：

(8)

(9)

(10)

式中：dijk为一定临界距离内，与斑块类型i相关的斑块类型j、k的连接状况；ni为斑块类型i的斑块个数；pij为斑块ij的周长，m；aij为斑块ij的面积，hm2；N为栅格总个数；eii为基于单倍法的斑块类型i像元之间的结点数；maxeii为基于单倍法的斑块类型i像元之间的最大结点数；CONNECT、COHESION、AI分别为生境连通性指数、生境整体性指数和生境聚集度指数，取值范围均为[0，100],除CONNECT无单位外，其余两个单位均为%.

遥感水质反演技术较常规实地采样水质监测方式在监测面积和监测时长上均具有显著优势[25]. 传统的实地采样水质监测方式的不足之处主要有三方面：一是所需人力成本和物质财力成本较高；二是采样所得数据在时间范围和空间范围上不连续，无法满足长时间序列下，区域大空间范围内持续监测的实际需求；三是因特殊原因在部分年份或部分采样点数据缺失时，单靠传统采样方式难以弥补部分数据缺失问题. 因此，传统监测方法难以满足流域管理与决策的实际需求. 近年来，为解决传统采样获取数据时空连续性不足、部分区域水质历史数据缺乏等实际问题，遥感水质反演技术被广泛运用于内陆水质监测[26-27]. 由于北运河水质数据来源限制，部分年份存在数据缺失问题，现有数据多为监测断面数据，缺乏空间和时间连续性，无法满足水生态承载力评估需要. 因此，该研究通过遥感反演的方法获取了6期水质及水生生物数据，反演指标包括COD、TN、TP、氨氮浓度以及叶绿素含量.

1.3.2水质指标遥感反演

已有较多研究运用实测水质数据和遥感数据进行水质遥感反演，主要反演的水质指标包括COD、TP、TN、BOD、氨氮浓度以及叶绿素含量等. 该研究通过2019年9月和2020年10月北运河的实测数据和部分历史数据，运用遥感反演模型获得1995年、2000年、2005年、2010年、2015年和2018年6期的COD、TP、TN、氨氮浓度以及叶绿素含量数据，从而完善水生态承载力评估体系. 2019年设置水样采样点34个，2020年加密到47个采样点，并补充采集了样点水体光谱信息. 光谱信息运用美国ASD公司(Analytical Spectral Devices,Inc.)生产的便携式地物光谱仪(ASD Field Spec 4)收集[28]. 监测样品的采集、保存、测定均按照相关技术规范处理[29-31].

在MATLAB 2017和Envi 5.3软件中结合实测光谱数据建立遥感影像和水质数据的反演模型. 在MATLAB 2017软件中对实测光谱数据进行平滑降噪、散射校正等预处理，并结合实验室实测水质数据进行特征波段的选取，然后对实测光谱和实测指标值运用最小二乘法进行建模，得到基于实测光谱数据的反演模型，并运用反演模型得到水质指标的预测值[32]. Landsat遥感影像数据经过辐射定标、Flash大气校正、裁剪等预处理，将所选择的反演模型特征波段范围与Landsat遥感影像数据波段范围进行对比，选择遥感影像中对应波段与实测光谱预测值运用线性回归方法进行建模，从而完成实测光谱数据对遥感影像数据的校正；将反演模型运用于各年份研究区遥感影像，从而得到北运河北京段范围内各年份水质数据[33-34]. 得到各年份水质数据后，参考常用的精度检验方法，选用决定系数(R2)和均方根误差(RMSE)并结合已有数据对反演精度进行评估[35]. 将反演出的水质指标值与从北京市北运河管理处收集到的历史水质资料进行对照，并且与《北京市水资源公报》中该区域的水质标准等级指示的水质指标范围进行对比[36].

2 结果与分析

2.1 水质指标遥感反演精度评价

实测指标数据与光谱数据的相关系数在0.85以上，遥感影像与实测光谱预测值建模得到的相关系数在0.8以上，COD、TP、TN、氨氮浓度以及叶绿素含量的相关系数(R2)和均方根误差(RMSE)如表4所示. 由表4可见，水质指标反演整体相关性较好，反演模型精度较高，因此可用于评价指标体系的计算.

表4 水质指标反演精度评价

《北京市水资源公报》显示，1986—2018年，北运河流域水质基本为Ⅴ类或劣Ⅴ类，将流域水质范围与各样点历年水质指标等级以及北京市北运河管理处部分水质等级历史数据进行对比，判定反演的水质指标等级基本与公布水质等级符合，因此反演精度与流域实际相符. 水利部和原环境保护部对北运河水质考核要求为Ⅴ类水质标准，由于承载力评估体系中水质指标反映的是与北京市市控断面等水质控制标准相比较的达标情况，通过水质超标(Ⅴ类水质标准)倍数来表征，因此反演精度与《北京市水资源公报》中区域水质等级基本一致即可满足评价精度要求.

2.2 北运河北京段水生态承载力时空变化分析

通过构建的北运河北京段水生态承载力评估体系，计算得到1995年、2000年、2005年、2010年、2015年、2018年6期区域整体和子流域的水生态承载力，结果如图2所示. 区域水生态承载力评估结果表明，北运河北京段水生态承载力在近30年间呈周期性变化. 北运河北京段水生态承载状态变化大致可分为3个阶段：1995—2000年、2000—2015年、2015—2018年. 1995年区域综合承载力分值为51.60，2000年区域综合承载力处于超载状态，分值为30.53，因此1995—2000年区域综合承载状态下降，2000年左右区域综合承载状态为近25年间最差. 2005年起区域承载状态有所恢复，区域承载力综合分值为52.43，为临界超载状态；2010年和2015年区域承载力综合分值分别为42.19和51.70，处于临界超载状态，2015年区域综合承载状态较2005年有所下降. 2015—2018年区域综合承载状态有所恢复，至2018年，区域综合承载分值达到63.91，已基本处于安全承载状态. 2018年区域整体水生态承载力较1995年和2015年高，说明近年来区域水生态承载力有上升趋势.

图2 北运河北京段总体水生态承载力时序图

从北运河北京段子流域承载状态(见图3)来看，北运河上游承载状态好于中下游. 近30年上游的沙河子流域、蔺沟河子流域以及中游的坝河子流域承载状态较好，均处于临界超载状态；而承载状态较差的子流域主要有凉水河、清河、温榆河，近30年间基本处于超载状态. 就具体阶段而言，蔺沟河自2010年以来承载状态下降趋势明显. 至2018年，北运河承载状态最好的子流域为坝河、沙河，水生态承载力分值分别为53.72、48.12；承载力最差的子流域分别为凉水河、北运河、清河，水生态承载力分值分别为28.59、30.33、30.97.

图3 北运河北京段子流域水生态承载力时序图

对各子流域综合承载力指标值作单因素方差分析以检验子流域间综合承载力的差异显著性，结果如图4所示. 进行单因素方差分析前首先进行了正态分布检验，Shapiro-Wilk检验和K-S检验结果显示，P值分别为 0.403 9 和0.956；方差齐性检验P值为0.894，在α=0.05的检验水平下，P>0.05，因此数据满足方差齐性及正态分布，可进行单因素方差分析. 单因素方差分析结果显示，P值为2.48×10-11，小于 0.000 1，表明各子流域水生态承载力差异性显著. 沙河、坝河的水生态承载力显著高于其他子流域，蔺沟河和小中河水生态承载力显著高于北运河、凉水河、温榆河和清河，清河的水生态承载力分值最低，长期处于超载状态. 沙河、蔺沟河子流域处于北运河水系上游，早期受人为扰动较小，因此整体承载状态好于其他子流域；坝河近年来实施了一系列的综合整治措施，生态恢复效果明显.

注： 不同的小写字母表示差异性显著.

2.3 区域水生态承载力关键要素识别

运用多层感知器神经网络进行建模、结合Pearson相关系数和Spearman相关系数方法得到区域水生态承载力的影响因子重要性，识别出区域水生态承载力关键要素. 结果(见图5)显示，相关系数与神经网络建模识别出的关键因子基本一致，重要性最大的影响因子为河岸带林草覆盖率. 多层感知器神经网络得出河岸带林草覆盖率标准化的重要性为100%，其次为河流连通性、人均水资源量、COD超标倍数、生境整体性指数，标准化的重要性分别为62.4%、52.4%、38.1%、34.3%. 从总指标层和分指标层来看，区域水生态承载力主要控制因子为水生态总指标中的水生生境指标，其次为水环境总指标中的子指标，水资源总指标的重要性较水生态和水环境指标低.

图5 水生态承载力影响因子重要性

3 讨论

3.1 社会经济因素对北运河水生态承载力的影响

该研究运用综合指标法从水资源、水生态、水环境3个层面构建了北运河北京段区域水生态承载力评估体系，结果显示，北运河北京段区域水生态承载力在近30年间呈周期性变化，即1995—2000年呈下降趋势，2000—2005年呈上升趋势，2005—2010年呈下降趋势，2018年较2015年呈上升趋势，并且2000—2010年水生态承载力的变化趋势与京津冀地区水生态承载力的变化趋势基本一致[37]. 已有研究[38-39]表明，水生态承载力反映的是生态系统对区域社会经济系统的承载能力，区域所能承载的人口、GDP必须与当地的资源、环境相协调，因而历年人口、社会经济指标的变化趋势也可以反映区域水生态承载力变化. 从区域人口变化来看，根据《北京市统计年鉴(1995—2018)》，2018年末北京市常住人口为 2 154.2×104人，较2015年的 2 170.5×104人有所下降，常住人口下降对缓解环境压力起到了一定作用，人口带来的环境压力变化与北运河整体水生态承载力变化趋势一致. 从区域经济发展速度来看，GDP环比增速与区域整体水生态承载力变化趋势基本一致，均在1995—2000年有所下降、2000—2005年有所上升、2005—2015年持续下降，但是2015—2018年GDP环比增速下降，区域整体水生态承载力有所提升，说明经济增速水平在一定程度上反映了区域生态承载能力，但是经济增长不一定就会造成区域水生态承载力的下降，这与关键控制要素识别结果相符. 因此，从环境、人口、社会经济等各角度综合分析，北运河北京段区域水生态承载力变化趋势与区域实际状况基本相符. 区域整体和部分子流域水生态综合状况的改善可能与北京市近年来采取的一系列水环境综合整治措施有关. 根据北京市北运河管理处相关资料，发现北京市从2005年开始进行黑臭水体修复、污水设施升级改造以及施行水环境区域补偿机制，因此近几年北运河部分区域水质得到了明显改善[36].

3.2 生态环境因素对北运河北京段水生态承载力的影响

沙河子流域和蔺沟河子流域均处于北运河上游，有研究[40-41]表明，这两个子流域地势海拔相对较高，土地利用类型以林地、草地为主，受人类干扰较小，生态环境基础状况较好，生态系统的自我维持与净化能力较强，水生态承载力较高，其中沙河子流域历年植被覆盖度均为各子流域中最高，除2000年以外均在60%以上. 据《北京市统计年鉴》数据显示，沙河子流域和蔺沟河子流域历年地均GDP和地均人口规模在所有子流域中属中等偏下，城市开发强度相对较低，说明沙河子流域和蔺沟河子流域受人为扰动较小，生态环境状况相对较好[40,42]. 但2010年后蔺沟河子流域水生态承载力呈明显下降趋势，北京市北运河管理处已有资料显示，2018年蔺沟河日排污量在 4×104m3以上，明显高于其他子流域，说明部分子流域可能存在治理不彻底的状况，从而导致子流域水生态的恶化. 此外，2010—2015年北运河水系中游坝河子流域水生态承载力呈上升趋势，2018年左右水系上游沙河子流域和中游坝河子流域水生态承载力分值明显较其他子流域高，这与已有研究结果[40,43]基本一致. 2010年前北运河水系中下游子流域承载状态较差，最差的子流域主要有清河子流域和温榆河子流域，较差的有凉水河子流域和小中河子流域，这些子流域均处于超载状态. 已有研究[44-45]表明，温榆河存在大量未达标污水排水口，清河和凉水河是北京市区主要排水河道，清河承载着城市化程度较高的北京市海淀区的人口、社会经济增长压力，而小中河沿岸工业污染排放较大，这些因素给区域水生态综合承载能力带来了较大压力. 2010年后除凉水河外，其他子区域水生态承载力均有不同程度改善，这一结论基本与北京市北运河管理处持续监测结果一致.

北运河北京段水生态承载力最主要的影响因素为水生态，如河岸带林草覆盖率. 《北京市通州区统计年鉴》显示，2005—2018年，通州区绿化覆盖率从34.87%升至45.48%，提升趋势显著. 已有研究[46-48]表明，生态环境的改善，尤其是区域植被覆盖率的提升对于区域生态输水能力具有重要影响. 也有研究[49-51]表明，区域植被覆盖度能调节气候，有效控制水土流失，有助于水质净化，并能进行雨洪滞蓄、涵养水源. 北运河水生态承载力受水环境影响也较为显著，如COD和TN浓度. 根据《通州区水环境质量限期达标规划(2018—2020年)》，通州区北运河流域污染来源主要为城镇生活污水、畜禽养殖污染、农村生活污水等，主要污染物为COD、氮和磷，各类污染物中COD排放量较为突出，这与关键控制要素中污染物排放规模等要素相符合. 此外，北运河流域水生态承载力也在一定程度上受到水资源变化的影响，主要影响因素为人均水资源量，其反映了流域水资源禀赋. 1995—2018年北京市水资源总量与北京市人均水资源量变化趋势基本一致，2000年北京市水资源总量为16.86×108m3，较1995年的47.62×108m3显著下降，2005—2015年呈周期性变化，2018年为35.46×108m3，较2015年的26.76×108m3有所增加. 2000年北运河北京段区域整体水生态承载力较1995年下降明显，2005—2015年呈周期性变化，2018年水生态承载力较2015年有所提升(见图3)，变化趋势基本与北京市水资源总量和人均水资源量一致. 北运河北京段区域水资源禀赋不佳，水资源补给以降雨为主，同时水资源总量年内变化明显，水资源十分宝贵；所以区域水资源综合利用效率的改善也可能是区域水生态承载力提升的重要途径[43,52]. 总体而言，北运河北京段区域仍需加强对河流沿岸湿地及植被的保护，并通过人工恢复与自然恢复相结合的方法进一步改善流域生态环境[52-53].

4 结论

a) 近30年来，北运河北京段水生态承载力呈周期性变化. 2000年区域整体水生态承载力分值为30.53，为历年最差，到2018年区域整体水生态承载力分值为63.91，呈上升趋势，但仍需加强水环境污染防治及水资源节约利用.

b) 北运河北京段各子流域间历年水生态承载力差异显著. 1995—2010年，水生态承载力最好的子流域为沙河子流域；2015年、2018年水生态承载力最好的均是坝河子流域. 除2015年水生态承载力最差的子流域为凉水河子流域外，其余评价年份水生态承载力最差的子流域均为清河子流域.

c) 通过对北运河北京段区域历年承载状态影响因素进行分析可知，水生态要素对区域承载状态影响最为显著，其次为水环境要素.