基于主成分分析与层次聚类分析的水质综合评价

邢 涛,雍 毅,侯 江, 吴 怡 ,吴 迪,刘恒博

(1.四川省生态环境科学研究院,成都 610041; 2.成都理工大学生态环境学院,国家环境保护水土污染协同控制与联合修复重点实验室,成都 610059)

前 言

泸沽属半封闭外流型湖泊，与外界河流交汇少，自净能力差[1]。同时其水体停留时间长，一旦污染短期难以恢复。近年来，泸沽湖生活污水排放量剧增，使得亮海湖滨带水质逐步降低。而草海为亮海出水泄流区，氮磷等污染物在此汇集并被水生植物所吸收。由于泸沽湖禁止打草收割，造成大量营养盐难以通过收割植物的方式清除，植物腐烂沉降后与底泥一起加剧内源污染，使得草海水质逐渐恶化[2]。因此，泸沽湖急需开展水质现状综合评价，并根据污染程度和特征分区治理。

水质评价是将定性水体状态问题通过多指标综合计算决策获得定量结果的过程[3]。因此，参与计算指标数的多寡将直接影响定量结果的科学性和准确性。目前，水质评价多采用综合营养状态指数法[4～6]，以叶绿素a(Chl.a)、总氮(TN)、总磷(TP)、高锰酸盐指数(CODMn)和透明度(SD)作为参评因子，可快速获悉水体大致营养状态。然而该方法各因子权重是根据中国26个主要湖泊调查的计算结果[7]，将其应用于某个具体水体不可避免有误差[8-9]。同时，其他水质因子也会对水体营养状态有一定的贡献率，如巢湖非离子氨和溶解氧对其水质富营养化贡献率分别为8.7%和14.5%[10]，天津滨海混盐水体清净湖的矿化度为水体富营养化藻类生长的主控因子[11]。而水质恶化后敏感指示生物数量的变化可能比Chl.a指标能更准确反映水质状态，如云贵高原湖泊特有的物种海菜花[12]。除此之外，水体植被的覆盖度、微生物、浮游动植物群落结构及生物多样性指数等生态指标也会影响水质评价结果[13-14]。如何将上述量纲不同的多种指标纳入水质评价目前尚无统一标准。

目前，水质评价常用指数法、神经网络法、灰色数学法和模糊评价法，但都无法很好地解决多指标问题[15]。主成分分析(Principal Component Analysis，PCA)优势在于通过降维将多指标提纯为少数几个综合指标，其计算过程基于相关系数矩阵，可消除各指标间量纲影响[16]，因此可将多指标纳入同一系统进行定量化研究。但PCA缺点是其结果不能反映提取出的综合指标间相互关系[17]。层次聚类分析(Hierarchical Cluster Analysis，HCA)优点在于可直观表现N维空间点与点之间的相互关系，但其缺点是聚类前需要保证各指标量纲统一，否则影响聚类结果[18]。因此将PCA和HCA结合使用能优缺点互补，多指标综合评价结果较为客观准确。

HCA是一种无监督学习过程，定义不同的类距离和相似性量度可产生不同的聚类结果，因而选择何种类距离进行最优化聚类仍无定论。目前，聚类距离常选择最短距离法和离差平方和法，但各研究对为何选择上述聚类方法并没有详细解释[19-20]。

因此，本实验以泸沽湖亮海和草海水体为研究对象，尝试将多种影响水质指标纳入同一体系进行综合评价。多指标采用主成分分析后，以主成分得分矩阵为基础进行层次聚类分析。在参考综合营养状态指数法评价结果的基础上，筛选出最优类距离下HCA结果。研究结果为多指标体系水质综合评价提供参考，并为泸沽湖水污染分区管理提供理论依据。

1 材料与方法

1.1 采样点布设

对泸沽湖亮海及草海沿线布设28个点位，其中亮海四川片区9个点位(点位1～8)、亮海云南片区5个点位(点位9～14)、草海沿线14个点位(点位15～28)。具体分布如图1所示。水样采集0.5m亚表层，部分点位水深不足0.5m则采集中间层水样。

图1 泸沽湖采样点位分布Fig.1 Distribution of the sampling sites in lake Lugu

1.2 监测项目与测定方法

本实验选取与水体富营养化相关的水质指标：Chl.a、pH、SD、溶解氧(DO)、TN、TP、氨氮(NH3-N)、CODMn。其中，NH3-N 是浮游藻类偏好吸收的氮形态。而藻类含量较多的水体，因其光合作用消耗大量CO2使得水体pH 升高。同时藻类释放较多O2，甚至能使水体处于DO过饱和状态。因此，将NH3-N、pH和DO指标纳入本次水质多指标评价体系中。Chl.a水样现场采集1L后，每升样品加入1mL碳酸镁悬浊液保存(10g/L)，随后放入0～4℃低温采样箱带回实验室分析，分析方法参照《水质叶绿素a的测定分光光度法》(HJ 897-2017)。SD采用塞氏盘法测定。pH和DO由便携式水质分析仪现场测定，其余水质指标同样取样1L，在现场加入酸保存剂后带回实验室检测，检测方法参考《水和废水监测分析方法》：总氮采用碱性过硫酸钾消解紫外分光光度法测定；总磷采用钼酸铵分光光度法测定；氨氮采用纳氏试剂比色法测定；高锰酸盐指数采用高锰酸盐滴定法测定[21]。

1.3 数据分析方法

1.3.1 总体水质判断

总体水质采用各水质参数在《地表水环境质量标准(GB3838-2002)》中所属水体类别进行判断，综合各个水质参数结果并以最差水体类别作为总体水质判断标准。

1.3.2 综合营养状态指数法

采用综合加权指数模型对水样Chl.a、SD、CODMn、TN和TP进行分析[7]，根据《湖泊富营养调查规范(第二版)》中的标准[22]，采用0～100的连续数值对各采样点水体营养状态进行分级。

1.3.3 主成分与聚类分析法

为消除多指标量纲和数量级的差异，对数据进行标准化处理。之后计算各指标相关系数矩阵，以相关系数矩阵为基础计算主成分，然后计算相关矩阵的特征值、方差贡献率、累积方差贡献率和因子在各主成分载荷矩阵。根据陡坡图以及累积方差贡献率大于85%的原则确定主成分个数。确定主成分个数后，根据各采样点位提取出的主成分得分矩阵进行层次聚类分析。采用欧式距离进行聚类分析，聚类方法综合对比最短距离法(Nearest neighbor method)、最长距离法(Furthest neighbor method)、离差平方和法(ward method)、类平均法(Group average)、中位数法(Median method)和重心法(Centroid method)，以综合营养状态指数结果为基础，筛选出直观且符合实际情况的聚类方法。以上过程采用OriginPro9.1软件进行分析。

2 结果分析

2.1 泸沽湖总体水质与综合营养状态

本次调查泸沽湖三片区总体水质如表1所示，亮海四川和云南片区大部分区域为Ⅱ类水体，但若以部分最差区域水质作为总体水质判断标准，则亮海四川片区最差为Ⅲ类水体，云南亮海片区最差为IV类水体，且两省亮海水质最差区域基本分布于亮海和草海交界处。其中，亮海四川片区S8(亮海草海交界处)点位水质为Ⅲ类，主要原因为水体CODMn含量较高。亮海云南片区点位S10(尼塞村)和S14(普洛码头)水质分别为Ⅲ类和IV类，主要原因为TN含量较高。而草海水质严重恶化，基本为Ⅲ～IV类水质，部分较严重区域甚至已达V类水质。S15(舍垮爱情码头)、S24(走婚桥阿六码头)和S27(母支码头)，三个点位TP含量均达到湖库标准V类，最严重区域S15采样点其水体DO和TN含量也达到湖库V类标准。

表1 泸沽湖总体水质情况Tab.1 Overall situation of water quality in lake Lugu

续表1

各采样点综合营养状态指数(TLI)如表2所示。亮海四川片区全为贫营养状态，而云南片区大部分区域为贫营养状态，部分区域S12和S14达到中营养状态，与总体水质判断结果相似。草海片区除S16、S17和 S18点位(接近亮海干净水体且出水水流较大)为贫营养水体外，其余点位大部分为中营养至轻度富营养状态。其中草海走婚桥附近及其下游区域污染较为严重，水体为轻度富营养状态。S15点位已达中度富营养状态，该点位位于草海南岸河口冲积扇，雨季农业面源污染较重且现场有生活污水和畜禽养殖废水汇集。

表2 泸沽湖各采样点营养状态评价结果Tab.2 Evaluation results of nutrition status at each sampling site in lake Lugu

续表2

2.2 泸沽湖水质主成分分析

2.2.1 相关矩阵

采用主成分分析法对泸沽湖28个点位水质指标进行分析，程序运行结果经标准化消除量纲后，变量间相关矩阵如表3所示。从表3可以看出，各参数间90%以上的数据绝对值大于0.1，表明全部变量至少与一个以上其他变量有较大的相关关系，适合用主成分分析法研究变量间关系。三种水质理化指标pH、SD和DO之间均有较好的相关性，其相关系数均大于0.8。营养性污染指标TN、TP和NH3-N之间也有较强的正相关关系。而有机污染指标CODMn与pH、SD和DO之间有较强的负相关关系(负相关系数绝对值大于0.8)，与营养性污染指标间有较弱的正相关关系，其中TP与CODMn正相关程度较高，相关系数为0.714。

表3 泸沽湖水质参数相关矩阵Tab.3 Correlation matrix of water quality parameters in lake Lugu

2.2.2 主成分提取

将泸沽湖28个采样点的9个(增加TN/TP比值)水质指标进行主成分分析后得到9个主成分，如表4所示。根据主成分提取一般要求，即方差累积贡献率超过85%的原则，提取前4个主成分作为泸沽湖水质参数污染状况信息概括，其方差累积贡献率高达93.59%。由陡坡图(图2)中本征值在第4个点出现凸点也可判断主成分数量选为4较为理想。而前两个主成分方差累积贡献率达到72.45%，表明前两个主成分代表了原始水质参数大部分信息。

表4 泸沽湖水质参数主成分提取分析Tab.4 Principal component extraction analysis of water quality parameters in lake Lugu

图2 泸沽湖水质参数陡坡图Fig.2 Scree Plot of water quality parameters in lake Lugu

根据主成分载荷矩阵(表5)可知，主成分1中最高正相载荷为CODMn，其次较高正相载荷为TN、TP和NH3-N，表明第一主成分基本反映了泸沽湖水质中的有机污染指标和营养性污染指标。主成分2中最高正相载荷为pH，其次依次为TN和SD，表明第二主成分基本反映了水质理化指标信息。而水体Chl.a含量则在第三主成分上有最高的正相载荷。根据载荷矩阵和主成分得分图生成的双标图(图3)可以看出营养性污染指标TN、TP和NH3-N三者和有机污染指标CODMn对主成分1有相似贡献率，水质理化指标pH、SD和DO对主成分2 有相似的贡献率。S15号采样点在主成分1和2上具有最高得分，表明该采样点各方面水质参数均受到严重污染，这与综合营养状态指数法计算出该点为最严重的中度富营养化结果一致。而S24和S26点位在主成分1上得分仅次于S15，与综合营养状态指数法计算结果判断该区为轻度富营养状态一致。

表5 泸沽湖水质参数主成分载荷矩阵Tab.5 Principal component load matrix of water quality parameters in lake Lugu

图3 主成分分析双标图Fig.3 Biplot of principal component analysis

2.2.3 层次聚类分析

根据主成分分析结果，计算出28个采样点四个主成分得分，代替原始水质指标进行层次聚类分析。采用欧式距离表示因子得分值间的相似性，并依次选择离差平方和法、最短距离法、类平均法、重心法、中位数法和最长距离法，以综合营养状态指数判定结果为基础，筛选最优聚类方法，聚类结果树状图如图4所示。

2.2.3.1 聚类法筛选采用重心法和中位数法进行聚类时，由于数据结构不适合，导致origin软件系统无法继续进行运算，因此运算结果只剩下离差平方和法(图4A)、最短距离法(图4B)、类平均法(图4C)和最长距离法(图4D)四种聚类方法。

图4 不同层次聚类分析方法树状图Fig.4 Dendrogram of different Hierarchical cluster analysis methods

结合综合营养指数结果，确定以类间距分别不大于7、2、3和6作为上述四种聚类法划分标准。以S13、S15、S25和S26点位作为参考，离差平方和聚类结果中污染最严重的S15号点位于中间，不能直观反映出各点位污染程度。若无综合营养指数结果参考，很难判断各点位大致污染状况。最短距离法将污染程度相似的S25和S26划分为两类，且对绿色聚类结果各点位不能很好的区分(该类别点位实际上包含了水体贫营养、中营养以及轻度富营养三个状态)。类平均法结果中处于贫营养状态的S13点位处于横轴右侧，与污染最严重S15号点位距离较近，同样不能直观反映各点位污染程度。最长距离法聚类结果不仅与综合营养状态指数法相符合，且能直观反映出各点位污染程度，因此作为本次聚类法最优筛选结果。

2.2.3.2 最优聚类结果

根据上述筛选结果，本实验采用最长距离法进行聚类分析，聚类树状图横轴从左至右反映了各点位从亮海贫营养向草海中度富营养过渡的趋势。28个点位总共可聚为5类：

S1～S12、S14和S16总共14个点位聚成1类(用类1表示)。类1代表了亮海大部分水质较好的贫营养区域。

S13号点单独聚为类2，该点位虽属贫营养水体，但由于其在主成分2所代表的水质理化指标上有远高于类1各点位得分，因此单独聚为一类，类1和类2可代表水体贫营养状态。S17～S22和S28总共7个点位聚为类3。类3除S17和S18为贫营养外，其余5个点位为中营养。从采样点分布图(图1)可以看出，S17～S22点位处于亮海干净水体逐渐深入草海沿线上，该区水质逐渐恶化，因此类3可以代表水体由贫营养逐渐向中营养恶化类别。

S23～S27总共5个点位聚为类4.类4中除S23和S27为中营养外，其余S24、S25和S26为轻度富营养状态。因此类4可以代表水体由中营养逐渐向轻度富营养状态恶化类别。

S15号点单独聚成类5。该点位为泸沽湖草海污染最严重区域，属于中度富营养状态。因此类5代表水体中度富营养状态。

3 讨 论

3.1 主成分分析优势

水质的综合评价由于涉及多个监测点位、多个相关水质指标和多个评价对象，同时每个指标在不同区域反映水质的角度和对总体水质变化的贡献率多有差异，因此它是一个多元复杂的体系[23]。主成分分析法由于采用相关系数矩阵而非协方差矩阵进行计算，消除了多指标间量纲不同对结果的影响，因而有望将富营养化水体水质理化指标、指示生物指标以及生态指标或其他量纲不同的相关因子纳入评价体系，获得更加综合全面的水质评价结果。同时，主成分正相载荷大小能反映出何种指标对水质影响较大，各污染区域是哪方面指标需重点关注治理，并以此制定有针对性的解决方案。目前已有较多的研究将主成分分析法应用于受多因素影响的水库[24]和江河[25]水质评价中。

3.2 层次聚类分析优势及筛选

虽然主成分分析法在多指标综合评价中有其优势，但其结果没有明确的范围，只能反映提取出的综合指标强弱大小关系，不能反映其所处位置[17]。部分研究试图通过提取出的各主成分方差贡献率作为权重，并乘以相应的主成分得分矩阵，得到各采样点综合得分，然后依照得分大小对所有采样点进行综合排名[26-27]。但该方法目前尚有争议[28]，需采用其他统计方法解决实际问题[19]。采用层次聚类分析法以主成分得分矩阵为基础，虽能很好的将主成分计算结果加以分类，明确各采样点所处位置。但由于层次聚类方法较多，如最短距离法、中位数法、最长距离法、重心法、类平均法和离差平方和法等，各聚类方法得出的结果一般情况下不会完全一致。目前，水质评价层次聚类分析法中应用离差平方和法较多，但对于为何采用该聚类方法没有详细说明[29-30]。因此如何衡量聚类结果是否最优目前尚无明确的标准，该聚类衡量标准的制定仍有待研究。本次实验中，我们以水质评价应用广泛的综合营养状态指数法为基准，筛选出最长距离法作为最优系统聚类法，获得的结果与实际水质营养状态相符，并能直观表现出各点位水体污染程度及相对位置，类1到类5依次代表水体从贫营养状态过渡到中度富营养状态。因此，泸沽湖分区治理应首先关注类5和类4水域及其周边污染情况。

3.3 适用性评价

泸沽湖水体由于呈现出明显的亮海水质较好，草海水质逐步恶化的趋势，因而28个采样点水质离散程度较好，采用较少的水质指标就能以基于主成分分析的层次聚类法加以区分。而对于总体水质相似，各采样点离散程度较差的区域(如某些水质均优良的西藏高原湖泊或均较差的云南滇池等湖泊)，能否以本实验采用的多指标综合水质评价方法，通过纳入更多的相关水质理化指标、指示生物指标和生态指标等加以区分，达到较完整和准确水质综合评价的目的，有待进一步检验。

4 结 论

4.1 泸沽湖亮海虽大部分处于贫营养水体，但总体已处于Ⅱ类水质。而草海沿线水质逐步恶化。部分区域如S15点位已达到中度富营养化状态，主成分分析显示其有机污染指标CODMn、营养性污染指标TN、TP、NH3-N和等水质理化指标DO均污染严重。

4.2 采用基于主成分得分矩阵的层次聚类分析法，以欧式距离下最长距离法能获得最优的聚类结果，并与综合营养状态指数法获得的结果相符，可作为水质层次聚类分析法筛选参考标准。同时，该聚类结果能直观反映出各点位水体污染程度及相对位置，聚类树状图横轴从左至右依次代表类1贫营养水体到类5中度富营养水体。分区治理应首先关注类5和类4代表的中度和轻度富营养化水域。