基于改进灰色关联度聚类的海水质量状况评价方法研究：以大亚湾为例

赵连玲 ，刘华雪，饶义勇，廖秀丽，戴 明，黄洪辉

1. 中国水产科学研究院南海水产研究所/广东省渔业生态环境重点实验室，广东 广州 510300

2. 西安电子科技大学 数学与统计学院，陕西 西安 710071

海洋环境污染会危害海洋生态系统，影响海洋食物产出，并对人类健康造成影响。海洋环境质量评价是开展海洋环境管理的依据，而海洋环境评价的精度取决于海洋环境的监测和评价水平。世界各国学者开展了海洋环境质量综合评价研究，其中海水水质评价方法、海洋沉积物评价方法和海洋环境综合评价方法是当前研究的重点。在有限的投入下精准评价海洋环境质量，一直是我国学者努力追求的目标。

自1998 年开始实施的海水水质标准 (GB 2097—1997)，根据海域的不同使用功能和保护目标，将我国海水水质分为4 类，并将各类水质的监测项目水质标准进行了定量划分，同时界定了各监测项目的海水水质分析方法[1]。Paul 等[2]对沿海水质的时空分布分析方法进行了总结，包括聚类分析、判别分析、因素分析和主成分分析、变差函数等方法。聚类分析被用于研究海水水质监测数据的相似性。因素分析和主成分分析主要用于沿岸和海水水质数据集的研究，通过减少冗余因子，获取主要的海水环境影响要素。

为了进一步改进评价方法，国内学者就指标体系构建和模型实验进行了大量尝试，如林小苹等[3]利用主成分分析，根据每个站点的主成分得分进行站点聚类，成功将8 个站点分为4 类。马丽[4]运用可变模糊评价模型，得到了5 个海湾的评价级别特征值和健康水平，并利用模糊综合评估模型进行了方法一致性比较。林琳等[5]从水体污染、初级生产力、生物资源等方面对大亚湾海域生态环境的质量状况进行了综合评价，并模拟和分析了环境质量的空间分布趋势。杜飞雁等[6]、方良等[7]、廖秀丽等[8]研究了大亚湾浮游动物的生物量、优势种等生态特征。党耀国等[9]通过构建面板数据的灰色关联度模型，发现基于面板数据灰色关联度模型的聚类方法能够在数据上表现出良好的效果。张燕军等[10]提出了一种基于灰色关联度的人机系统功能聚类方法，详细阐述了人机系统的分层聚类方法，并通过案例验证了该方法的可行性和有效性。岑詠霆[11]根据旅行社顾客满意度，确定旅行社样本及评价指标体系，计算灰色关联度，最后以离差平方和最小为聚类标准，进行聚类分析。郭三党等[12]定义了不同类别之间的关联度，构造了一种基于最大灰色关联度的聚类方法，结果表明新的聚类方法改进了已有的灰色关联聚类中聚类原则存在的缺陷。

大亚湾是我国亚热带海域重要的海洋生物种质资源库[13]，受周边地区工业和人为影响，大亚湾海域生态环境遭到破坏，生态系统也在快速退化[14-15]。越来越多研究者开始关注陆源性污染物的输入以及近海人工养殖造成的海湾环境污染问题[16-18]。海水养殖是大亚湾的重要产业之一，促进了该海域水体中营养盐和有机物的富集，局部海域已出现富营养化趋势[14]。由于大亚湾附近没有江河湖泊等较大的陆源径流输入，湾内海水交换条件一般，主要通过湾口与外海水进行交换[19-24]，水体交换能力总体为南部优于北部、东部优于西部，全湾水体更新时间为26 d[25]。大亚湾核电基地是我国目前在运行核电装机容量最大的核电基地，拥有6 台百万千瓦级压水堆核电机组[26]，大亚湾核电站和岭澳核电站的温排水也成为了影响研究海域生态环境的重要因素之一[27]。

本文通过改进灰色关联度聚类，结合层次聚类，构建涵盖水环境和海洋生态的指标体系，在大亚湾西部海域 (核电附近海域) 开展海水质量评价探索研究，将统计学方法引入海域质量评价，旨在建立科学有效的指标评价体系，准确分析目前大亚湾西部的海水质量状况，进而为行业主管部门提供对策和建议。

1 材料与方法

1.1 数据来源

本文所有数据均来自2020 年对大亚湾西南海域的科考调查结果，调查范围为114°30'E—114°55'E、22°25'N—22°50'N，包含2020 年四季24 个站点的数据 (图1)，主要包括海水水质环境、海洋生态(浮游植物、浮游动物) 调查。

图1 大亚湾监测站点Fig. 1 Monitoring stations in Daya Bay

海水水质环境测量用容积为5 L 的有机玻璃采水器采集表层水样 (水面下50 cm)，各调查项目按规范要求现场分装保存，部分指标于实验室进行分析测试。水深用船载渔探仪或FISH4200 型便携式渔探仪现场测量，透明度用萨氏盘法现场测定，海水水温、盐度、pH 和溶解氧 (DO) 用YSI Pro Plus型多功能水质仪现场测定，化学需氧量 (CODMn)、石油类、活性磷酸盐、无机氮、悬浮物和叶绿素a参考海洋监测规范第4 部分：海水分析 (GB 17378.4—2007) 进行测定。浮游植物用浅水III 型浮游生物网距海底约2 m 至表层垂直拖网进行采集，网口系流量计，每站垂直拖曳1 网，所采样品用5% (φ)的甲醛溶液固定，带回实验室经浓缩后，用数字流式细胞摄像系统 (FlowCam 8400) 进行镜检分类鉴定与计数，分析浮游植物数量分布和多样性。浮游动物用浅水I 型浮游生物网距海底约2 m 至表层垂直拖网进行采集，网口系流量计，每站垂直拖曳1 网，所采样品用5% (φ)的甲醛溶液固定，带回实验室挑去杂物后，以湿质量法称取浮游动物生物量，再用解剖镜进行镜检分类鉴定和计数，分析浮游动物数量分布和多样性。浮游植物密度单位为105个·m−3，Shannon-Weaver多样性指数计算公式为：

式中：Pi=ni/N；ni为单位空间内第i种的个体数量 (个·m−3)；N为某站总生物数量 (个·m−3)；S为出现生物总种数。

1.2 技术路线

首先，收集24 个站点的指标数据形成指标体系；其次，对指标体系进行主成分分析，验证指标体系的合理性；再次，运用改进的灰色关联度聚类，根据指标体系对24 个站点进行聚类，绘制聚类图；最后，根据站点聚类结果，得出结论，分析原因。

1.3 评价方法

1.3.1 主成分分析

主成分分析是采用数学上的降维思想，用少数几个综合指标来代替原来指标的一种多元统计方法。转化生成的综合指标被称为主成分，其中每个主成分都由原始变量的线性组合生成，能够尽可能多地反映原来指标的信息，且各个主成分之间互不相关。

假定有n个样本，每个样本共有p个指标 (变量) 描述，这样就构成了一个n×p阶的数据资料矩阵：

其中：

作X1,X2,···,Xp的线性组合即综合指标，也就是主成分，记新变量指标为Z1,Z2,···,Zp。则：

在上述方程组中要求：a21i+a22i+···+a2pi=1,i=1,2,···,p，且系数aij由下列原则来决定：

1)Zi与Zj(i/=j,i,j=1,2,···,p) 不相关；

2)Z1是X1,X2,···,Xp的一切线性组合中方差最大者；Z2是与Z1不相关的X1,X2,···,Xp的所有线性组合中方差最大者；Zp是与Z1,Z2,···,Zp−1都不相关的X1,X2,···,Xp的所有线性组合中方差最大者。

这样决定的新变量指标Z1,Z2,···,Zp分别称为原变量指标X1,X2,···,Xp的第一、第二、… 、第p主成分。

1.3.2 灰色关联度聚类

以灰色关联分析结果R为基础，构建基于灰色关联度的聚类方法，从而达到对研究对象进行分类的目的。该方法的实质就是生成如下二维空间到一维空间的映射f:R2→Rg，其中R为关联度集，Rg为评价对象的相似关系集。

1) 确定聚类论域的原始评估矩阵。设待评估海域有m个站点，记为聚类论域S=(S1,S2,···,Sm)；评价生态系统健康的指标集，记为C=(C1,C2,···,Cn)，得到原始评估矩阵：

式中：序列Yi=(yi1,yi2,···,yin) 代表海域内第i个站点的各项生态系统健康指标集。

2) 数据的标准化处理。由于站点各生态指标的量纲和界限值不同，保证指标间相同因素的可比性，消除不同指标间的不可公度性影响，需要对原始评估矩阵中的数据进行无量纲化处理，采用正态标准化，按下式作同一化处理：

式中：i为指标代号；j为生态系统健康指标代号；为第j个指标的均值；Sj为第j个指标的标准差；所有数据经过同一化处理后，得到评估矩阵X：

式中：序列Xi=(xi1,xi2,···,xin) 代表海域内第i个站点经过数据标准化处理后的各项生态系统健康指标集。

3) 聚类基础的构成。根据分析问题的实际情况，需要设定基准序列X0=(x01,x02,···,x0n)，为减小误差，对于极大型指标：x0j取该指标列的90%分位点处数据，对于极小型指标，x0j取该指标列的10%分位点处数据，对于区间型指标，x0j取该指标列的50%分位点处数据或均值。

对于比较序列X1,X2,···,Xm，令Rij=|x0j−xij|，则有：

式中：xij与x0j的相对差值ξ0ij定义为关联系数；ρ是分辨系数，在 (0,1) 之间取值，可以人为给定，ρ越小分辨力越大，通常取ρ=0.5。关联系数很多，信息较为分散，它的每一个值表示了比较序列Xi对基准序列X0的关联程度强弱。为更加直观地表现2 个数列的关联度，对所有关联系数取平均值：

式中：ri代表了序列Xi对基准序列X0的关联度，称为绝对关联度，ri值越大，说明序列Xi与基准序列X0的相似性越高，即对应站点的生态系统越健康。因此可以得到聚类论域S的关联度集R：R=(r1,r2,···,rm)，关联度集R是聚类分析的基础。

4) 建立灰色相似关系矩阵。根据关联度集R，构建聚类论域S中各元素之间的关联度差异矩阵ES：

式中：eij(i,j=1,2,···,m) 为Xi相对于Xj的差异系数：

由关联度差异矩阵ES得差异距离矩阵DS：

式中：dij为差异距离，dij=eij+eji。

由于矩阵ES的主对角线为零，所以矩阵DS是主对角线为零的对称矩阵。

根据差异距离矩阵DS计算可得灰色关联度矩阵Rg：

且有：

式中：max(DS) 表示取矩阵DS中的最大值元素。对矩阵Rg，显然满足：自反性gii=1，对称性gij=gji。因此Rg是灰色相似关系矩阵。经过上述步骤，便完成了映射f:R2→Rg。

5) 聚类分析。灰色相似关系矩阵Rg反映了分析论域S中各元素相互间的关系亲疏程度，可以按照矩阵Rg利用最大树方法对论域S进行聚类。首先以所有分类的对象为顶点，从大到小依次连接灰色相似关系矩阵Rg中的元素gij，在不产生回路的基础上将所有顶点连通，生成最大树并绘制谱系图。

6) 改进的灰色关联度聚类。基于灰色相似关系矩阵Rg的灰色关联度聚类形式单一，缺点明显，类似于单链接的层次聚类，为了改进这一缺点，将灰色相似关系矩阵Rg转化为距离矩阵W：

式中：

然后利用层次聚类的最小平方差算法，计算每个簇中每个点到合并后的簇中心的距离差的平方和，得到最终的聚类结果。

2 结果与分析

2.1 主成分分析结果

将表层海水水质环境数据和浮游生物数据、表层叶绿素数据进行主成分分析 (表1)，其中KMO 值为0.727＞0.5，说明所用数据适合做主成分分析。Bartlett 球形度检验，近似卡方=812.282(P＜0.01)，说明因子变量之间高度相关，能够为主成分分析提供合理的基础。按照特征值大于1 的标准，确定主成分数，前6 个主成分方差累计贡献率为66.78%。由主成分分析结果，按照因子载荷绝对值＞0.4 可以提取4 个主成分 (表2)，第一个主成分PC1 为海洋水文和生物生态相关指标，第二主成分PC2 为石油类、叶绿素a、无机氮指标，第三主成分PC3 为透明度和活性磷酸盐指标，第四主成分PC4 为浮游动物生物量和CODMn指标。

表1 KMO 和Bartlett 的检验Table 1 KMO and Bartlett tests

表2 海水质量评价的因子载荷矩阵Table 2 Factor load matrix for seawater quality assessment

2.2 聚类分析结果

聚类选取的指标体系与主成分分析保持一致。首先需要确定每个指标的类型，指标数值越大则水质越好的确定为极大型指标；指标数值越小则水质越好的确定为极小型指标；指标数值波动于某个区间的确定为区间型指标。pH、DO、CODMn、无机氮、活性磷酸盐、悬浮物、石油类、水温参考海水水质标准 (GB 3097—1997)，浮游植物多样性、浮游动物多样性、浮游植物密度、浮游动物生物量参考近岸海域海洋生物多样性评价技术指南 (HY/T 215—2017)，盐度参考近岸海洋生态健康评价指南(HY/T 087—2005)，叶绿素a、透明度参考廖秀丽等[28]。此外，根据海水水质标准 (GB 3097—1997)，一类、二类水质是人为造成的海水温升夏季不超过当地1 ℃，其他季节不超过2 ℃，所以，温度相对越低越好，故温度取极小型指标。评价指标及指标类型总结如表3 所示。

表3 海水质量评价指标及指标类型Table 3 Seawater quality assessment index and index type

根据数据可以计算出4 个季节及年度水环境评价的绝对关联度 (表4)，其中绝对关联度越大，综合评价越优。

表4 海水质量评价的绝对关联度Table 4 Absolute correlation degree of seawater quality assessment

海水水质标准 (GB 3097—1997) 中，将水质分为4 个等级；近岸海域海洋生物多样性评价技术指南 (HY/T 215—2017) 中，将海洋生物划分为5 个等级；近岸海洋生态健康评价指南 (H Y/T 087—2005) 中，将生态健康划分为3 个等级；廖秀丽等[28]将海湾养殖环境划分为3 个等级。本文将海水质量状况划分为3 个等级 (良好、中等、较差)。根据距离矩阵W，进行最小平方差算法的层次聚类，4 个季节及年度的聚类结果如图2 所示。

图2 海水质量评价聚类图及站点分区Fig. 2 Cluster map and site partition of seawater quality assessment

根据春夏秋冬四季和年度水环境评价的聚类图，可将站点划分为3 类 (表5) ，一类 (良好) 站点占比20.83%，二类 (中等) 站点占比66.67%，三类 (较差) 站点占比12.5%，总体结果表现中等。4 个季节相比，平均绝对关联度为夏季 (0.747 9)＞冬季 (0.734 5) ＞春季 (0.729 0) ＞秋季 (0.709 7)，夏季海水质量状况相对较好，秋季表现相对最差。

表5 海水质量评价的分类Table 5 Classification of seawater quality assessment

2.3 整体评价结果

大亚湾海洋环境质量受近岸人类开发活动和南海外海水入侵共同影响，其中大亚湾海域潮流基本呈往复流流态。大亚湾湾口的大辣甲和黄毛山将大亚湾湾口分隔成3 条通道与外海相连。涨潮时，湾外潮流沿3 条通道进入大亚湾内，然后朝北往大亚湾顶上溯。大辣甲西侧上溯的水流在大坑附近海域辐散，分成两支潮流分别向西进入大鹏澳和向东北随主流向北进入澳头湾。落潮时，大亚湾内的水体基本上沿涨潮流相反方向流出湾口。根据评价结果，结合大亚湾水流方向及人类活动影响方式的不同，可将大亚湾24 个站点分成3 个区域 (A 区4 站：S9、S10、S19、S20，B 区12 站：S4、S6、S7、S11—S13、S15—S17、S21—S23，C 区8 站：S1—S3、S5、S8、S14、S18、S24，图2)，而这3 个区域的划分与大亚湾周边实际生产生活状况和笔者的实地踏勘情况相符。本文评价指标并未包含重金属、有机污染物和渔业资源等成本较高的监测指标，说明本方法在基于有限指标的情况下，展现出了一定的操作性。

A 区临近大亚湾北部近岸海域，北部沿岸城镇密集、人口众多，并建有大亚湾石油化工基地，该基地是我国为数不多的世界级石化基地，工业产能较大，对近岸海域环境有较大影响，因此A 区可定义为大亚湾北部工业和城镇影响区。大亚湾核电基地附近海域呈现出明显的辐散、辐合流特征，所以B 区可定义为辐散、辐合流及温排水影响区。C 区S2 站点附近是大亚湾重要的水产养殖区，主要为网箱养殖与贝类养殖，所以S2 站点主要受到养殖场的影响。其他站点均离岸较近，且S8、S14 站点位于旅游景点杨梅坑附近，在旅游旺季陆源输入的影响较大。S24 站点位于大亚湾湾口的西侧海域，代表了外海水进入大亚湾的区域，受外海水影响较大，所以C 区可定义为人类活动影响区。

总结四季各站点质量评价的分类结果见表5。

3 结论

根据构建的指标体系，通过改进灰色关联度聚类，结合层次聚类，对大亚湾海域24 个站点进行聚类评价，成功将海域分成3 个区域，而这3 个区域的划分与大亚湾周边实际生产生活状况以及笔者的实地踏勘相符，实现了在有限指标情况下的评价可操作性。

在聚类结果的展示图方面，层次聚类只能展示出同类站点，不能展示出站点的相对优劣顺序，而灰色关联度聚类可以展现站点的相对优劣顺序，所以在将灰色关联度与层次聚类结合时，可考虑在改进算法中增加站点相对优劣顺序，增强聚类结果的展示性。