基于聚类分析和主成分分析的长江下游稻田氮磷监测指标筛选

胡冰涛，陈玉东

(生态环境部南京环境科学研究所，江苏 南京 210042)

长江下游是我国水稻种植最为广泛的地区，过量的化肥施用在保障水稻产量的同时，导致大量氮磷营养盐通过稻田径流和淋溶等途径进入河流、湖泊和地下水，造成严重的水体富营养化和地下水污染[1-4]。选用适宜的氮磷指标是监测稻田氮磷流失状态和评估其对水环境污染影响程度的重要环节。利用原始变量开展监测工作不仅工作量大且重叠度高[5-6]。运用多元统计技术，从众多原始变量中筛选出相对独立的敏感性指标，通过建立最小数据集(minimum data set，MDS)，可有效减少监测指标数量，降低监测工作量，且监测指标的信息量更为集中[7]。目前，该方法已经在土壤质量[8-9]以及水环境[10-12]监测和评价中被广泛应用。

金卫斌等[13]将湖北四湖流域沟渠6项水质指标进行聚类分析，筛选出代表水体氮、磷负荷的主要监测指标，并将沟渠水质划分为4个污染等级；JUAHIR等[14]运用层级聚类分析将越南冷岳河(the Langat River)上、下游按污染程度分区，采用正反向逐步判别分析从各区域23项水质指标中提取出6～7项主要监测指标，用以评估冷岳河流域水质空间分布特征；YANG等[15]、SHRESTHA等[16]、KOWALKOWSKI等[17]均运用多元统计方法筛选水质评估指标，并进一步解释水环境指标间复杂关联且有效识别污染肇因和解析水质的时空动态变化。

当前稻田氮磷面源污染评估研究多侧重于稻田氮磷流失、渗漏负荷及影响因素的分析[9,18-22]，而关于简化稻田氮磷流失监测步骤和提高污染负荷估算精度的指标筛选研究相对较少。仅有谢文明等[20]针对长江三角洲区域提出4个简易监测指标(田面水铵态氮浓度、日最高气温、土壤120 cm深度总氮浓度和降水量)，计算该区域稻田氨挥发总量、氮磷淋溶损失和氮磷径流损失，其建立的简易监测指标及构建的数学关系显著提高了面源污染发生量预测精度。

以上海市青浦区现代农业园内水稻田田面水和不同深度淋溶水为研究对象，基于1 a内全生育期田间试验，获取8个氮、磷类指标数据，通过主成分分析和聚类分析分别建立最小数据集，并采用水质综合得分方法[23-24]比较、验证最小数据集指标的合理性。研究结果可为长江下游地区稻田氮磷污染监测指标体系优化提供技术参考。

1 材料与方法

1.1 实验地概况

试验于2018年6—11月在上海市青浦区现代农业园区(31°12′10″ N，121°08′12″ E)进行，该地处于亚热带潮湿型季风气候区，年均温为15.3～16.7 ℃，年均降水量为1 150～1 550 mm。降水主要集中在3—9月，占全年总降水量的75%左右。试验年稻季平均气温为26.3 ℃，降水量为595.3 mm。试验地土壤类型为湖沼相沉积物起源的青紫泥水稻土，质地为重壤土。耕作层(0～20 cm)土壤m(土)∶V(水)=1∶5，pH值为7.33，w(有机质)为18.74 g·kg-1，w(全氮)为0.94 g·kg-1，w(全磷)为1.17 g·kg-1，w(全钾)为1.67 g·kg-1，土壤电导率(EC)为0.14 mS·cm-1。当地水田常种植单季稻或水稻-小麦两熟轮作。

1.2 实验材料

随机选择青浦区现代农业园区6块等面积水稻田，实验面积共计1 200 m2，供试水稻品种为“花优14”。种植密度为1.13×105株·hm-2。试验过程中的关键农事操作如下：6月29日施基肥，6月30日移栽水稻，记为移栽后0 d，7月12日(12 d)施分蘖肥，8月6日(37 d)追施穗肥，11月8日水稻成熟收割。水稻生长期间基肥施BB肥(又名散装掺混肥，以N计)100.0 kg·hm-2，分蘖肥施BB肥(以N计)75.0 kg·hm-2，穗肥施尿素(以N计)75.0 kg·hm-2，均以表面撒施形式施用。磷肥为过磷酸钙，施用量(以P2O5计)为120 kg·hm-2，钾肥为氯化钾，施用量(以K2O计)为130 kg·hm-2，磷钾肥均以基肥形式于移栽前与耕层土壤充分混合后一次性施入。

1.3 监测指标和分析方法

1.3.1淋溶水和田面水取样

在水稻种植前，每块实验田均在距离田埂3 m处并排垂直安装2根不同取样深度的渗漏管(直径为5 cm，长为90 cm)，底部封闭，分别在管壁30和60 cm处(上下5 cm范围内)每间隔2 cm 均匀钻 6 排直径为5 mm 的渗水孔，将收集的渗漏水视作30和60 cm处的混合水样。用塑料纱网将渗水孔段管壁包裹扎牢，以阻隔淤泥进入管中。为防止田面水顺管壁下渗，沿管壁裹缠塑料薄膜并向四周水平外延约20 cm，然后在其上回填原层次土。PVC管管口高于土壤表面20 cm，以避免田面水过多时管口被淹没。取样时，将1 m 长的软橡皮管(内径为6 mm) 插入渗漏管中，用100 mL医用注射器将管中渗漏水全部抽出作为待测水样。每次田面水取样时，不扰动土层，按照五点采样法抽取田块内5点各100 mL田面水溶液，组成混合样品，取混合样200 mL保存。水稻移栽后，分别于基肥、分蘖肥和穗肥施用后7 d内每天上午采集田面水和渗漏水样。

1.3.2检测指标

采用流动分析仪(AA3，荷兰Skalar公司)测定田面水和淋溶水样中总氮(TN)、溶解态总氮(DTN)、硝态氮(NO3--N)、铵态氮(NH4+-N)以及总磷(TP)、溶解态总磷(DTP)和磷酸盐(PO43--P)，颗粒态磷(PP)为TP与DTP的差值。

1.4 数据处理与统计分析

采用 Excel 2010进行数据处理与制图，采用PASW Statistics(SPSS)19.0对数据进行统计分析。

1.5 最小数据集的建立

最小数据集的构建有多种方法，如线性回归分析、多元回归分析、主成分分析、聚类分析[7]、判别分析[5]和评分函数等方法。建立最小数据集后，最小数据集指标与未入选指标相关性分析法[25]和质量指数比较法[26]是当前验证最小数据集合理性的主要方法。

1.5.1聚类分析

聚类分析(cluster analysis，CA)以欧氏距离来衡量评价指标之间差异性，将各指标分类并建立最小数据集，在水质评价[5,27]中应用广泛。采用SPSS 19.0运用欧氏距离最短法对评价指标进行R型聚类分析。在相应聚合水平之间，将评价指标划分为反映水质污染不同方面特征的若干分组。根据相关分析，每组中相关性显著的指标可相互代替，并与野外调查结果、文献资料和前人研究成果相结合，剔除冗余指标，选择具有代表性且相互独立的指标进入最小数据集[7]。

1.5.2主成分分析

主成分分析(principal component analysis，PCA)是建立最小数据集的主要方法[16]。通过降维将多个指标转化为少数指标，消除众多指标间的相关性，使彼此之间不相关。提取特征值≥1的主成分，将同一主成分载荷≥0.5的指标分为1组，计算评价指标的Norm值。Norm 值为该指标在由主成分组成的多维空间中矢量常模的长度，长度越长，表明该指标在所有主成分的综合载荷越大，其解释综合信息的能力就越强。Norm值计算公式为

(1)

式(1)中，Nik为第i个变量在特征值>1的前k个主成分上的综合载荷；Uik为第i个变量在第k个主成分上的载荷；λk为第k个主成分的特征值。如果所有主成分都被选择，则该Nik值为1。分别计算各组指标Norm值，选取每组中Norm值在最高总分值10%范围的指标，进一步分析每组所选指标的相关性。若同组所选指标高度相关，则Norm值最大的指标进入最小数据集(MDS)；若相关度低，则同组指标均进入MDS。相关系数>0.5，表明高度相关。

1.6 水质评价方法的建立

(1)原始数据标准化。原始样本矩阵X为

X=(xij)n×p，i=1，2，…，n；j=1，2，…，p。

(2)

式(2)中，n为采样次数；p为基于不同数据集使用的指标数。

数据标准化处理按式(3)进行。

(3)

(2)相关系数矩阵。计算所对应的相关系数矩阵R，并计算R的特征值和特征向量。

R=(rij)p×p，

(4)

其中，

(5)

λi(i=1,2,…,p)为R的特征值，将其按大小排列λ1≥λ2≥…≥λp≥0,即为主成分的方差，其大小对应主成分对原始样本的权重。每个特征根对应的特征向量为w1,w2,…,wp,通过特征向量将标准化指标转化为主成分。根据主成分方差累计贡献率大于80%～90%或特征值大于1的原则确定主成分个数m。

(3)主因子得分和每个主因子的方程贡献率。根据各数据集前m个主成分的主因子得分和每个主因子的方程贡献率，求出指标权重。相关公式为

Fj=β1jx1+β2jx2+…+βpjxp，

(6)

ej=wj/(w1+w2+…+wm)，

(7)

(8)

式(6)～(8)中，Fj(j=1，2，…，m)为主成分；x1，x2，…，xp为各个指标；β1j，β2j，…，βnj为各指标在主成分Fj中的系数得分；ej为Fj的方程贡献率；ωi为指标xi的权重。

(4)水质综合得分(water component analysis,WQI，Iw,q)的计算公式为

(9)

式(9)中，xi为数据集中第i个指标值；ωi为第i个指标权重；p为数据集中指标个数。

1.7 水质评价精度验证

利用Nash有效系数(Ef)评价基于主成分分析最小数据集的精确程度，其计算公式为

(10)

2 结果与讨论

2.1 稻田田面水和淋溶水氮磷动态变化

如图1～3所示，在水稻生长期内，田面水以及30和60 cm深度淋溶水中各形态氮素浓度均出现3个峰，分别对应施基肥并移栽、施分蘖肥和施穗肥时间。氮肥投入后，田面水中TN、DTN和NH4+-N浓度迅速提高，1 d时即达到峰值，之后随着时间推移迅速下降，10 d内下降到较低浓度并保持稳定。而30和60 cm深度淋溶水各项氮素指标变化对于施肥的响应较慢，较田面水峰值出现时间迟3～6 d，这可能是由于在施肥初期，氮肥水解生成的NH4+-N通过水流向下淋溶，土壤胶体一般带负电，带正电的NH4+-N在稻田中易被土壤胶体通过交换吸附作用固定，迁移至下层土壤过程较为缓慢[4,20]。田面水以及30和60 cm深度淋溶水各形态氮素浓度峰值表现为依次递减。淋溶水中各形态氮素浓度均值表现为60 cm深度TN、DTN和NO3--N平均质量浓度为3.81、1.57和0.68 mg·L-1，略高于30 cm深度TN、DTN和NO3--N质量浓度(3.68、1.38和0.49 mg·L-1)，而30 cm深度淋溶水中NH4+-N质量浓度为0.83 mg·L-1，高于60 cm深度NH4+-N质量浓度(0.53 mg·L-1)。生育期内田面水中TN、DTN、NO3--N和NH4+-N平均质量浓度分别为7.04、6.41、1.54和4.54 mg·L-1。

图1 稻田田面水各形态氮素和磷素浓度动态变化Fig.1 Dynamic changes of nitrogen and phosphorus concentrations in paddy field surface water

如图2所示，仅在施基肥时施入磷肥，导致水稻移栽后，田面水中各项磷素指标达到峰值，之后变化幅度较大，但在生育期前期整体呈现浓度下降趋势。施加穗肥后，田面水中磷素指标出现明显峰值，这可能是由于施穗肥(仅氮肥)提高了田面水氮素浓度，田面水中微生物和藻类增殖，提高了磷素转化相关酶的活性，从而加速土壤磷的溶解和生物周转[28]。另外，土壤胶体和铁、铝氧化物表面的吸附点位被施加的氮占据，土壤对磷的吸持能力减弱，磷素向田面水中释放[29]。

图2 稻田30 cm深度淋溶水中各形态氮素和磷素浓度动态变化Fig.2 Dynamic changes of nitrogen and phosphorus concentrations in the leaching water at the depth of 30 cm in paddy field

图3 稻田60 cm深度淋溶水中各形态氮素和磷素浓度动态变化Fig.3 Dynamic changes of nitrogen and phosphorus concentrations in the leaching water at the depth of 60 cm in paddy field

由图2～3可知，30和60 cm深度淋溶水中各形态磷素浓度只在施基肥6～7 d时出现1次峰值，随后一直保持在较低浓度，低于0.5 mg·L-1。经统计，生育期田面水中TP、DTP、PO43--P和PP平均质量浓度分别为0.267、0.135、0.072和0.132 mg·L-1，30 cm深度淋溶水中TP、DTP、PO43--P和PP平均质量浓度分别为0.655、0.405、0.151和0.250 mg·L-1，60 cm深度淋溶水中TP、DTP、PO43--P和PP平均质量浓度分别为0.797、0.123、0.048和0.674 mg·L-1。

2.2 基于聚类分析建立的最小数据集

基于聚类分析法中的R型聚类，对稻田监测数据的全量数据集8个指标进行分类。由图4所示，对于稻田田面水和30 cm深度淋溶水，聚合水平在15～20之间时，指标可明显分为2类，TP、PP、DTP和PO43--P属于第1类，表征水样中磷素指标水平，DTN、NH4+-N、TN和NO3--N属于第2类，表征水样中氮素指标水平。聚合水平在10～15之间时，指标细分为4类，其中，TP和PP归为一类，磷素指标中DTP和PO43--P归为一类，DTN和NH4+-N组成氮素指标中的一类，余下的TN和NO3--N组成一类。60 cm深度淋溶水的聚类分析结果与前两者有较大差异，聚合水平在15～20之间时，TP、PP、DTP和PO43--P组成一类，TN、DTN和NH4+-N组成一类，NO3--N单独划分为一类。

图4 水稻田面水以及地下30和60 cm淋溶水氮磷指标聚类分析树状图Fig.4 Dendrogram of cluster analysis on water quality index of surface water and eluviation water at the depth of 30 and 60 cm in rice field

在太湖流域开展的相关研究[29]表明，在不同氮肥管理下，水稻生育期内田面水中PP为磷素主要形态。根据全量指标间的相关性，稻田田面水中TP和PP呈高度相关，故选取PP代替TP进入最小数据集；PO43--P为DTP的主要组成部分，两者呈明显正相关，故选取DTP进入最小数据集。氮素指标中NO3--N在10～15聚合水平上与DTN、NH4+-N、TN分类差异明显，选取NO3--N进入最小数据集；田面水中DTN、NH4+-N和TN这3者之间相关显著，叶玉适等[30]研究表明，淹水稻田的土壤和田面水硝化作用都较弱，大部分时间内水稻生育期田面水中氮素以NH4+-N为主，笔者研究中生育期内NH4+-N平均浓度占TN浓度的比例达64.49%，为氮素主要成分，故选取NH4+-N进入最小数据集。

根据30 cm深度淋溶水氮磷指标聚类分析结果，聚合水平在15～20之间时，NO3--N被单独划分为一类，直接进入指标集合。在水稻生长期内，淋溶水中溶解态总磷(DTP)为磷素主要形态[31]，选取DTP代替PO43--P进入指标集合。鉴于TP与DTP呈高度相关，选择与TP同组的PP进入指标集合；NH4+-N在土壤中迁移能力较弱，且土壤处于水分饱和强还原状态，减少了NH4+-N向NO3--N的转化，NH4+-N是氮素的主要存在形态[32]；选取PP、DTP、NH4+-N和NO3--N作为30 cm深度淋溶水的水质评价监测最小指标集合。

由于长江下游地区稻田土壤肥沃，有机质和黏粒含量较高，对PP有较强的吸附作用，难以随水分运动进入淋溶水[4]，DTP成为该地区淋溶水中磷素主要形态。在60 cm深度淋溶水的各项指标中，聚合水平为15～20时DTP与PO43--P组合为一类，选取DTP进入监测最小数据集；TP与DTP高度相关，以PP代替TP进入最小数据集；聚合水平在15～20之间时，NO3--N被单独划分为一类，故最小数据集应包含NO3--N；而NO3--N与DTN相关性较高，选取TN进入最小数据集。

基于聚类分析确定水质监测指标，田面水和30 cm深度淋溶水监测指标具有一致性，均为DTP、PP、NH4+-N和NO3--N，60 cm深度淋溶水的聚类分析筛选指标则为DTP、PP、TN和NO3--N。

2.3 基于主成分分析建立的最小数据集

由表1～3可知，田面水和60 cm深度淋溶水水质评价指标中特征值大于1的有3个主成分，累计贡献值分别为87.913%和85.126%；30 cm深度淋溶水指标中特征值大于1的有4个主成分，累计贡献值为85.780%，均满足信息提取要求。田面水中TP、DTP、PP和PO43--P进入第1组，TP的Norm值最高，且与DTP、PP和PO43--P相关性高(决定系数R2分别为0.671**、0.561**和0.515**)，则TP进入最小数据集。第2组中DTN与TN高度相关(0.822**)，结合Norm值考虑剔除TN和NH4+-N。第3组中仅有NO3--N直接进入最小数据集。

表1 田面水水质指标载荷矩阵和Norm值Table 1 Load matrix and Norm value of the surface water quality evaluation

30 cm深度淋溶水中全量指标依据Norm值计算，可划分为4组。第1组由TN、DTN和NH4+-N组成，DTN的Norm值为组内最高，且分别与TN和NH4+-N高度相关(0.639**和0.911**)，故代表第1组入选最小数据集。第2组指标中PP和TP相关性高(0.896**)，且TP的Norm值大于PP，故PP被剔除，DTP与PO43--P组成第3组，其中，PO43--P的Norm值高于DTP，且与DTP高度相关(0.853**)，并且DTP位于PO43--P的Norm值10%范围内，选择PO43--P入选最小数据集。与田面水相似，NO3--N被单独划分为一组，进入最小数据集。

60 cm深度淋溶水中8项氮磷指标按照载荷矩阵和Norm值计算，可分为3组，第1组包含TN、DTN和NH4+-N，其中，TN与NH4+-N高度相关(0.921**)，且Norm值均位于DTN的10%变化范围内，故选取DTN代表第1组入选最小数据集，NO3--N、TP和PP构成第2组，根据Norm值剔除PP，NO3--N与TP相关性较低(0.257**)，故NO3--N与TP共同进入最小数据集。结合Norm值和指标间相关性，第3组中DTP被确定为最小数据集组成指标。

表2 30 cm深度淋溶水水质指标载荷矩阵和Norm值Table 2 Load matrix and Norm value of the water quality evaluation at the depth of 30 cm

表3 60 cm深度淋溶水水质指标载荷矩阵和Norm值Table 3 Load matrix and Norm value of the water quality evaluation at the depth of 60 cm

综上分析，基于主成分分析结合Norm值计算筛选所得田面水监测指标为TP、DTN和NO3--N，30 cm深度淋溶水监测指标为TP、DTN、NO3--N和PO43--P，60 cm深度淋溶水监测指标为TP、DTN、NO3--N和DTP。

2.4 最小数据集合理性的验证

比较基于不同数据集的水质综合得分变化范围、均值，水质综合得分与全量数据集水质指数的相关性，以及Nash系数，从而验证最小数据集评价指标体系的合理性，是检验指标筛选结果的重要环节。田面水全量数据集水质指数(WQI-TDS)变化幅度为-0.838 6～1.722 3，平均值±标准差为0.002 62±0.003；基于主成分分析的田面水水质指数(WQI-PCA)在-0.661 2～3.035 8之间变化，平均值±标准差为0.002 32±0.002；基于聚类分析的田面水水质指数(WQI-CA)变化区间为-0.529 1～2.828 7，平均值±标准差为-0.004 49±0.05(图5)。

图5 全量数据集水质综合得分(WQI-TDS)分别与基于主成分分析(WQI-PCA)和聚类分析的最小数据集水质综合得分(WQI-CA)的相关性Fig.5 The correlation between the comprehensive score of water quality of the total data set (WQI-TDS) and the comprehensive score of the minimum data set based on principal component analysis (WQI-PCA) and cluster analysis (WQI-CA) respectively

30 cm深度淋溶水的WQI-TDS变化范围为-0.578 4～1.969 4，平均值±标准差为-0.004 37±0.001 2，WQI-PCA范围为-0.589 9～2.034，平均值±标准差为-0.003 39±0.000 9；WQI-CA范围为-0.653 4～3.138 1，平均值±标准差为-0.003 7±0.007。

60 cm深度淋溶水WQI-TDS、WQI-PCA和WQI-CA区间分别为-0.639 5～2.849 7、-0.676 1～2.159 6和-2.179 18～7.508 6，平均值±标准差各为-0.001 25±0.004、-0.001 47±0.001和-0.007 7±0.09。

上述结果表明田面水以及30和60 cm深度淋溶水WQI-PCA变化幅度及均值较WQI-CA更接近于WQI-TDS。从拟合结果(图5)来看，WQI-TDS分别与WQI-PCA、WQI-CA呈现显著正相关，其中，田面水R2分别为0.861和0.717，30 cm深度淋溶水R2分别为0.877和0.840，60 cm深度淋溶水R2分别为0.861和0.865，总体上WQI-PCA拟合结果较WQI-CA更优。田面水以及30和60 cm深度淋溶水WQI-TDS与WQI-PCA的Nash系数分别为0.861、0.876和0.856，而WQI-TDS与WQI-CA的Nash系数分别为0.712、0.796和0.787，结果表明采用基于最小数据集-主成分分析(MDS-PCA)对稻田田面水和淋溶水水质评价的准确性较基于最小数据集-聚类分析(MDS-CA)更高，故选择MDS-PCA代替全量数据集对稻田氮磷污染和水质情况进行评价。

3 结论

由于受肥料施用、降雨、地形、土壤性质和水稻生产方式等自然 或人文因素的影响，稻田面源污染的产生和排放具有分散面广、随机性强、不易监测和难以量化等特点，使得稻田面源污染水质监测变得复杂和繁琐。目前，在科学研究过程中稻田面源污染监测实践仍然存在着监测方法繁琐、监测指标众多和监测工作耗时费力等缺点，使得监测工作局限于小范围、短时间内进行，无法满足实际农业环境的监测和管理需要。

对上海市青浦区现代农业园内水稻田田面水以及地下30和60 cm深度淋溶水氮磷指标的全量数据集和最小数据集监测指标进行分析和评估比较。分别应用主成分分析和聚类分析从稻田田面水以及30和60 cm深度淋溶水的8项氮磷污染监测指标中筛选出3～4项作为监测指标最小数据集。基于不同数据集的稻田田面水和淋溶水水质评价结果及相关性差异明显，应用基于主成分分析最小数据集的水质评价结果在变化范围、均值方面比基于聚类分析的最小数据集表现得更接近于全量数据集结果，同时，前者与全量数据集的Nash系数及相关性系数均高于后者。研究结果表明采用基于主成分分析的最小数据集(MDS-PCA)对稻田田面水和淋溶水水质评价的准确性较基于聚类分析(MDS-CA)的最小数据集更高，可选择MDS-PCA代替全量数据集对稻田氮磷污染和水质情况进行评价。基于最小数据集筛选出的监测指标数量明显少于科学研究中氮磷监测指标，能简化监测环节，减少人力、物力和时间的投入，为环保部门实施农业面源污染水质监测提供便捷的方法。