近年来，生态环境监测领域取得了显著进展，特别是地表水环境质量监测网得到了全面升级，已建成1 837座水质自动监测国控站[1]，同时各省市及地方配备了相应数量的省控站和地方站，形成了覆盖广泛、功能完善的监测网络。水质监测内容涵盖水温、pH、溶解氧（DO）、电导率（EC）、浊度（Tur）5项基本参数，以及氨氮（NH3-N）、高锰酸盐指数（CODMn）、总氮（TN）、总磷（TP）4项重要指标，为实时了解水质变化趋势提供了重要数据[2]。然而，面对如此庞大的监测数据集，如何高效地进行数据分析，并深入挖掘其中的潜在价值，揭示水污染的特征及变化规律，成为亟待解决的重要课题之一[3]。
随着研究的深入，水环境评价方法日益丰富，而面对海量的水质连续监测数据，常用方法如单因子评价法、综合污染指数法、主成分分析法及聚类分析法等存在评价指标受限、评价结果不全面、评价指标选择及权重系数分配主观，以及解释性不足等局限[4-5]。因子分析法作为一种有效的数据分析工具，能够对原始变量信息进行重新组合，并利用旋转技术增强因子变量的可解释性，用于识别河流污染因子及其贡献率[6-7]。该分析方法的计算过程相对复杂，需要专业的数学知识和特定的软件支持，常规的数据分析软件难以满足其评价需求。为此，本研究利用Python语言，结合因子分析法，对江苏某省控站2021—2023年的水质监测数据进行深入分析，提取水质污染关键因子，为水生态环境保护及区域可持续发展提供参考。
1 材料与方法
1.1 数据来源
选取2021—2023年江苏某省控站的水质监测数据作为研究对象，排除异常值和缺失值（Nan值）后，共获得6 237条有效数据。
1.2 研究方法
1.2.1 因子分析数学模型
因子分析是一种多元统计分析方法，旨在以最少的信息丢失，将大量原有变量精简为少数几个相互独立的因子，并使因子具有一定命名解释性。该方法能够实现高维数据向低维数据转化。假设原始数据中有p个观测变量x 1，x 2，x 3，…，xp，每个观测变量用k（k<p）个因子f 1，f 2，f 3，…，fk 的线性组合来表示。
$\begin{array}{l}x_{1=} a_{11} f_{1}+a_{12} f_{2}+a_{13} f_{3}+\cdots+a_{1 k} f_{k}+\varepsilon_{1} \\x_{2=} a_{21} f_{1}+a_{22} f_{2}+a_{23} f_{3}+\cdots+a_{2 k} f_{k}+\varepsilon_{2} \\x_{3=} a_{31} f_{1}+a_{32} f_{2}+a_{33} f_{3}+\cdots+a_{3 k} f_{k}+\varepsilon_{3} \\\cdots \\x_{p=} a_{p 1} f_{1}+a_{p 2} f_{2}+a_{p 3} f_{3}+\cdots+a_{p k} f_{k}+\varepsilon_{p}\end{array}$
式中，xi （i=1，2，…，p）是观测变量，fj （j=1，2，…，k）是彼此不相关的公共因子；aij （i=1，2，…，p；j=1，2，…，k）是载荷矩阵，是第i个原有变量在第j个因子上的载荷；εj （j=1，2，…，k）为特殊因子，表示原有变量不能被因子解释的部分，其均值为0，独立于fj （j=1，2，…，k）。
1.2.2 因子分析评价方法
本研究以Python3.12作为开发工具，并应用NumPy、Pandas、Matplotlib、Scikit_learn、Factor_analyzer等第三方库进行开发[8-9]。NumPy是Python进行科学计算的基础包，包含强大的N维数组对象，支持高效的数组和矩阵运算。Pandas是开源的数据分析和处理库，提供了易于使用的数据结构和数据分析工具。Matplotlib库提供了丰富的绘图功能，是Python的数据可视化库。Scikit-learn是开源机器学习库，提供了大量的机器学习算法和数据处理工具。Factor_analyzer库用于执行因子分析，提供计算因子载荷矩阵、解释方差、获取方差贡献率等功能。
（3）特征根及方差贡献率计算。特征根又称方差，是表示因子影响力度的指标，通常用于评估因子的有效性或重要性，一般按照特征根大于1的原则提取主因子[14]。方差贡献率表示主因子保留并解释评价指标信息的百分比，采用Factor_analyzer库的FactorAnalyzer方法计算特征根、方差贡献率及方差累计贡献率。
2 结果与分析
2.1 KMO及Bartllett's检验
表1 各指标数据信息统计 |
| 统计项 | 水温/℃ | pH | DO/（mg/L） | EC/（μS/cm） | Tur/NTU | CODMn/（mg/L） | NH3-N/（mg/L） | TP/（mg/L） | TN/（mg/L） |
|---|---|---|---|---|---|---|---|---|---|
| 数据量（count）/条 | 6 237 | 6 237 | 6 237 | 6 237 | 6 237 | 6 237 | 6 237 | 6 237 | 6 237 |
| 均值（mean） | 20.03 | 8.07 | 8.24 | 474.15 | 14.43 | 5.48 | 0.17 | 0.05 | 1.25 |
| 标准差（std） | 8.85 | 0.40 | 2.80 | 90.45 | 8.99 | 0.98 | 0.15 | 0.02 | 0.74 |
| 最小值（min） | 3.90 | 7.22 | 0.07 | 207.60 | 2.40 | 0.42 | 0.01 | 0.01 | 0.02 |
| 25% | 12.20 | 7.74 | 6.03 | 418.50 | 8.30 | 4.89 | 0.08 | 0.03 | 0.76 |
| 50% | 20.70 | 8.04 | 8.40 | 467.70 | 12.60 | 5.38 | 0.13 | 0.04 | 1.06 |
| 75% | 28.10 | 8.37 | 10.47 | 530.60 | 17.70 | 5.98 | 0.20 | 0.06 | 1.65 |
| 最大值（max） | 36.90 | 10.32 | 18.76 | 669.60 | 102.80 | 10.51 | 2.09 | 0.23 | 14.05 |
2.2 特征根及方差贡献率
按照特征根大于1的原则，本研究提取的3个主因子以64.29%的累计方差贡献率保留并解释原始评价指标（表2）。因子旋转后，3个主因子对全部指标信息的方差贡献率分别为25.39%、20.36%、18.54%，其中，第一个主因子占主导地位，保留并解释的信息最多，对全部指标信息的影响最大。
表2 旋转前后特征根、方差贡献率及累计贡献率 |
| 主因子 | 因子旋转前 | 因子旋转后 | ||||
|---|---|---|---|---|---|---|
| 特征根 | 方差贡献率/% | 累计贡献率/% | 特征根 | 方差贡献率/% | 累计贡献率/% | |
| F1 | 2.233 | 27.91 | 27.91 | 2.032 | 25.39 | 25.39 |
| F2 | 1.509 | 18.86 | 46.77 | 1.629 | 20.36 | 45.75 |
| F3 | 1.401 | 17.52 | 64.29 | 1.483 | 18.54 | 64.29 |
2.3 主因子载荷
根据主因子载荷矩阵绘制散点图（图1），F1主因子为水质综合污染指标，反映水体中有机物及营养元素的污染状况，是研究断面水质变化的重要影响因子，其高载荷值表明F1主因子在水质变化中起着核心作用。F1相关联的评价指标TP（loading为0.874）、CODMn（loading为0.648）、NH3-N（loading为0.552）是评价水体富营养化和有机物污染的重要指标，其浓度升高通常与有机物分解有关；Tur（loading为0.58）间接反映了有机物和营养元素的污染[22]，Tur增加通常与水体中悬浮物增多有关。TP、CODMn、Tur、NH3-N与F1主因子高度正相关，表明TP、CODMn、Tur、NH3-N是研究断面水质变化的主要影响因子，TP、CODMn、Tur、NH3-N值越高，表明研究断面的水质越差。F2主因子为水质基本污染指标，反映水体的基本理化性质，相关联的评价指标pH（loading为0.861）是衡量水体酸碱性的重要指标，也是影响水生动植物成活率及生长率的重要因素[23]，与F2主因子高度正相关。研究断面pH在7.22~10.32，呈现弱碱性至碱性，pH越大，对水体中微生物活动的抑制作用越强[24-25]，从而导致水质劣变。DO（loading为-0.522）反映水生态系统结构及水体污染程度，DO值低于2 mg/L会导致鱼类及大多数藻类死亡[26]，研究断面DO与F2主因子存在较强的负相关关系，说明溶解氧减少通常反映出水质状况的恶化。F3主因子为溶解性盐类污染指标，相关联的评价指标TN（loading为0.827）包括氨氮、硝酸盐氮等多种形式；EC值（loading为0.746）较高可能反映出水体中溶解性盐类或有机物含量增加，TN和EC值越高，说明水质越差。
2.4 综合评分
根据因子得分系数矩阵（表3）及观测变量的监测值，建立主因子的计算方程。
$\begin{array}{l}\quad F_{1}=0.040 x_{1}+0.004 x_{2}-0.181 x_{3}+0.260 x_{4}+0.386 x_{5}+ \\0.247 x_{6}+0.430 x_{7}+0.106 x_{8} \\\quad F_{2}=0.540 x_{1}-0.319 x_{2}-0.050 x_{3}-0.106 x_{4}+ \\0.395 x_{5}-0.267 x_{6}+0.008 x_{7}+0.006 x_{8} \\\quad F_{3}=-0.042 x_{1}-0.013 x_{2}+0.493 x_{3}-0.161 x_{4}+ \\0.125 x_{5}+0.272 x_{6}-0.020 x_{7}+0.566 x_{8} \\\quad F=\left(\lambda_{1} F_{1}+\lambda_{2} F_{2}+\lambda_{3} F_{3}\right) /\left(\lambda_{1}+\lambda_{2}+\lambda_{3}\right)\end{array}$
式中，F 1~F 3为3个主因子得分，x 1~x 8分别表示pH、DO、EC、Tur、CODMn、NH3-N、TP、TN标准化后的监测数据，F为综合评分，λ 1~λ 3为3个主因子的特征根。相应数据代入主因子计算方程后得到6 237个综合评分。
表3 因子得分系数矩阵 |
| 水质指标 | F1 | F2 | F3 |
|---|---|---|---|
| pH | 0.040 | 0.540 | -0.042 |
| DO | 0.004 | -0.319 | -0.013 |
| EC | -0.181 | -0.050 | 0.493 |
| Tur | 0.260 | -0.106 | -0.161 |
| CODMn | 0.386 | 0.395 | 0.125 |
| NH3-N | 0.247 | -0.267 | 0.272 |
| TP | 0.430 | 0.008 | -0.020 |
| TN | 0.106 | 0.006 | 0.566 |
综合评分变化趋势显示（图2），研究断面水质波动明显（Score=-1.60~2.38，std=0.37），显示出水质的不稳定性和易受外界因素影响的特点。研究断面2021—2023年逐年的水质变化趋势总体相似，每年存在2次明显的水质波动，主要分布在2—3月（冬春季）与6—10月（汛期），特别是在6—10月，波幅明显增大（Score=-1.60~1.97，std=0.46），是全年水质较差的时段，表明研究断面水质受到相似因素的影响，且季节性超标特征明显[27]。研究断面位于闸前，水体流动性较差，藻类生长茂盛，周边分别有大片农田及养殖塘，其水质波动可能与季节变化、藻类暴发、养殖尾水排放、农田退水及汛期降水等因素有关[28-30]。从逐年看，研究断面2021年水质波动最大（std=0.38）；2022年水质波动最小（std=0.34），水质最好，跨度时间短（Score=-0.11，Score<0的时段达65.87%），可能与2022年降水量较常年偏18%有关；2023年水质最差，波幅明显，跨度时间长（Scor=0.10，Score<0的时段达45.78%），可能与2023年汛期降水量较同期偏多有关。
为验证因子分析评价结果，同步开展单因子评价，并绘制单因子评价法与因子分析法水质评价结果叠合图（图3）。根据单因子评价结果，研究断面水质变化幅度较大，水质类别在Ⅱ类至劣Ⅴ类范围内波动，各月份均有不同程度超标，超标时间主要集中在2—3月（冬春季）与7月前后（汛期），呈现季节性超标特征，其结果与因子分析综合得分结果基本吻合。单因子评价法的评价结果直观且易于理解，能够快速识别出超标因子和超标时段，但水质类别受限于单一指标的波动，评价结果严苛[31]，且受限于评价标准，无法全面反映水质的整体状况[32]。因子分析法综合了多个单因子评价法受限的指标（如EC、Tur、TN等），同时，综合考虑了参评指标的贡献，使得评价结果更具综合性，在描述水质波动形态、污染程度及污染持续时间方面展现出明显优势。
3 结论
（2）本次研究提取的3个主因子以64.29%的累计方差贡献率保留并解释原始评价指标，其中F1主因子为水质综合污染指标，反映水体中有机物及营养元素的污染状况，相关联的指标为TP、CODMn、Tur、NH3-N，是研究断面水质变化的重要因子。F2主因子为水质基本污染指标，相关联的评价指标为pH、DO。F3主因子为溶解性盐类污染指标，反映水体受溶解性盐类污染的程度，相关联的评价指标为TN、EC。
（3）研究断面水质波动明显（Score=-1.60~2.38，std=0.37），波峰主要集中在2—3月（冬春季）与6—10月（汛期），表明该时段内水质较差，呈现出明显的季节性特征；6—10月波幅明显增大（Score=-1.60~1.97，std=0.46），是全年水质较差的时段，该结论与单因子评价法的结果相吻合，验证了因子分析法的有效性和准确性。
（4）因子分析法能够涵盖如EC、Tur、TN等单因子评价法可能受限的指标，同时综合考虑参评指标的贡献，使得评价结果更具综合性，其在描述水质波动形态、污染程度以及污染持续时间等方面展现出明显优势。2021—2023年研究断面逐年的水质变化趋势总体相似，表明研究断面水质受到相似因素的影响。其中，2021年水质波动最大（std=0.38）；2022年水质波动最小（std=0.34），水质最好，跨度时间短（Score=-0.11，Score<0的时段达65.87%）；2023年水质最差，波幅明显，跨度时间长（Score=0.10，Score<0的时段达45.78%）。