基于PCA-SOM的北京市平谷区地下水污染溯源

林斯杰，齐永强，杨梦曦，杨 庆，杨梦凡, 刘 毅，胡 清*

1.南方科技大学环境科学与工程学院，广东 深圳 518055 2.清华大学环境学院，北京 100091 3.北京市地质矿产勘查院，北京 100195 4.北京环丁环保大数据研究院，北京 100083

地下水是北京市的重要供水水源，平谷区又是北京市的重要应急水源地. 随着平谷区工农业发展和人口增长，区域地下水水质呈恶化趋势[1]. 地下水具有自净能力弱、补给交互慢、污染隐蔽性强、难以逆转等特点，一旦污染将会给人类带来巨大危害. 因此，开展平谷区地下水污染源及污染特征研究，对及时控制污染源、采取针对性管控措施、保障北京应急水源地安全具有重要意义. 然而，目前对平谷区地下水的研究多集中在水位变化分析[2-4]、水化学特征分析[5]和水质评价[6]等方面，较少研究对平谷区地下水污染来源进行系统分析.

国内外学者采用多种方式对地下水污染溯源进行了研究，总体上可以分为模型法、实测法和统计法[7]. 模型法是根据含水层中实测污染物浓度，通过数值-机理模型反求污染源时间及空间分布[8]，或是由响应结果推求模型输入参数进行时间反演计算[9-10]. 实测法多采用氮、氧、碳等同位素作为示踪剂，基于化学质量平衡和多元统计分析测算示踪剂的时空分布，进而推算污染物来源[11-13]. 统计法建立在水质监测数据基础上，依靠图论[14]、相关性分析[15]、灰色关联分析法[16]、模糊数学法[17]、PCA (主成分分析法)[18-19]、结合GIS图像识别污染源[20-22]、回归分析方法[1]等对污染物进行分析. 这些方法虽然可以表征和识别污染指标的污染来源，但分析维度较为单一，未能体现指标间关联性，不能很好地反映污染物空间分布和来源特征.

SOM (Self-Organizing Map，自组织映射)是通过神经元间的竞争学习将相似的数据映射到相邻位置并保持拓扑结构不变的一种方法. 由于这种聚类过程是数据的全部维度映射，SOM本身具有多维数据融合的属性. 相较PCA (主成分分析法)，SOM通过对输入模式的反复学习，使权重向量空间能反映输入模式的统计特征. 该方法被应用于地下水污染评价和污染源溯源，其分析效果直观、结果较好[23-27]. 然而，SOM在训练前需要预先定义神经元的映射结构和映射大小，同时对结果的解读依赖于神经元激活强度特征图谱的目视比较，相比PCA无法进行定量分析[28]. 而PCA无法对主要因子之外的污染指标进行解释，因此将PCA和SOM两种方法相结合(简称“PCA-SOM法”)既能解释主要因子涉及的污染指标关联性，又能对所有污染指标来源进行解释.

该研究采用PCA-SOM法，通过PCA对研究区域地下水主要污染因子进行同源性分析，利用监测指标在神经元网络上映射特征图像对污染物的相关性进行分析，对平谷区地下水污染物来源及污染特征进行解析，以期为平谷区开展地下水污染防治工作提供参考依据.

1 研究方法

1.1 研究区域概况

选取北京市平谷区为研究区域. 平谷区位于39°30′N～40°30′N、115°45′E～117°10′E，其北、东、南三面环山，西南开口与华北平原相连接. 区内海拔由东北部的85 m向西南部逐渐降至20 m，坡降为1.8‰～2.3‰. 平谷区是北京市暴雨中心地区之一，据平谷区国家气象观测中心1959—2018年观测资料显示，多年平均降水量约635.5 mm. 降水量具有时间、空间分布不均匀的特征，其主要集中在5—10月，占全年降水量的80%以上.

平谷区是典型的农业区，近年来形成种植业稳步发展、林业小幅下降、牧渔业主动调减的发展格局；同时，平谷区也零星分布有工业园、高新技术产业园区、生活垃圾填埋场等潜在风险源. 平谷区内尚存有2个规模较大的集中式地下水饮用水水源地，主要地下水资源以第四系松散堆积物中孔隙水为主[29].

1.2 数据来源

注：极富水区，岩性为卵石含漂石； 富水区，岩性为砾卵石； 中等富水区，岩性为卵砾石夹粗砂； 弱富水区，岩性为砂、卵石夹中粗砂； 极弱富水区，岩性为粗中砂夹卵砾石.图1 平谷区浅层地下水污染监测井位置和岩层富水性分布[6]Fig.1 Location of shallow groundwater pollution monitoring wells and water abundance of rock stratum in Pinggu District[6]

研究区内共计33眼浅层地下水污染监测专井[6](见图1)，监测年份为2010—2018年，采样时间分别为每年的3月、6月、9月和12月. 水质采样方式为人工采集水样，检测指标为32项关键指标，分别为ρ(Al)、总碱度、ρ(As)、ρ(Ca2+)、ρ(Cd)、ρ(Cl-)、ρ(CN)(CN为氰化物)、ρ(CO2)、ρ(CO32-)、ρ(Cr6+)、EC(电导率)、ρ(Cu)、ρ(F-)、ρ〔Fe(Ⅱ)〕、ρ〔Fe(Ⅲ)〕、总硬度、ρ(HCO3-)、ρ(Hg)、ρ(K+)、ρ(Mg2+)、ρ(Mn)、ρ(Na+)、ρ(NH4+)、ρ(NO2)、ρ(NO3-)、ρ(DO)、pH、ρ(Se)、ρ(SO42-)、ρ(TDS)(TDS为溶解性总固体)、ρ(VP)(VP为挥发酚)、ρ(Zn). 分析指标依照GB/T 14848—2017《地下水质量标准》[30]要求执行.

1.3 数据分析方法

使用基于MATLAB开发的地下水污染预警工具包——GasHD 1.0软件,其内置的SOM算法针对全局数据进行补齐. 使用SPSS 22.0软件自带的PCA分析工具进行多元统计分析. 采用Pearson检验进行相关性分析，p<0.05作为变量间有显著相关性的判据. 采用回归法作为主因子分析中提取主因子后的旋转方法，并对提取的各主因子意义进行分析. 在GasHD 1.0软件中使用SOM映射全部维度数据并进行聚类分析，结合区域地表人为活动情况、区域水文地质情况分析污染物来源和空间分布特征. SOM聚类数量指标通过平谷区水文地质条件划分和聚类站点映射图得到.

2 结果与讨论

2.1 地下水监测指标类型及统计特征

平谷区地下水化学成分如表1所示. 平谷区浅层地下水水质状况较好，pH范围为7.0～8.3，整体呈弱碱性，部分区域地下水水质超过GBT 14848—2017《地下水质量标准》[30]Ⅲ类标准限值，仅适用于农业和部分工业用水，不可直接作为生活饮用水. 与2018年相比，2011年枯、丰水期阳离子浓度平均值变化较小，呈小幅波状起伏；丰水期ρ(Na+)、ρ(K+)最高值均有上升，ρ(Ca2+)、ρ(Mg2+)均呈下降趋势. 由表1可见，ρ(NH4+)、ρ(SO42-)、ρ(Cd)、ρ〔Fe(Ⅱ)〕、ρ(NO2)变异系数较大，分别为3.402、2.967、5.021、2.700和5.836，5个指标空间分布离散性和差异性较大，存在局部富集现象.

对监测指标间进行Pearson检验表明，地下水所有化学元素中EC与ρ(Ca2+)(p=0.936)、总碱度与ρ(HCO32-)(p=0.981)、ρ(Mg2+)与总硬度(p=0.944)指标之间呈显著相关，表明平谷区地下水水化学类型主要以HCO3-Ca型为主，其次为HCO3-Mg型. 姜体胜等[5]通过Piper三线图、Schoeller图、Gibbs图、离子比率图等方法得出了相同结论. 由表1可见，ρ(Ca2+)、ρ(Mg2+)、ρ(HCO3-)的变异系数较小. 由于平谷区地下水系统围岩岩性以奥陶系白云岩为主，白云岩主要以Ca、Mg元素为主，因此平谷地区地下水中Ca2+、Mg2+、HCO3-为自然起源.

2.2 污染源解析

剔除数据分布不呈正态分布的指标后还剩29个指标. 采用KMO和Bartlett球形检验对研究区域地下水剩余的29个指标的标准化数据进行相关矩阵检验(见表2)，其中KMO值为0.685，接近0.7，适合进行因子分析；Bartlett球形检验P值接近于0，满足P<0.05的置信度，表明各变量之间具有较高的相关性. 该研究共提取8个公因子，累计方差贡献率约71.550%，能反映29项监测指标的信息.

表1 平谷区地下水化学成分统计表

注：加粗体为变异系数大于2的指标.

为了对公因子典型指标进行更好的解释，使用凯撒正态化最大方差法将荷载因子矩阵进行正交旋转(见表3). 针对公因子中高得分和包含污染物指标的公因子进行分析. 公因子F1主要包括ρ(Ca2+)、ρ(Mg2+)、ρ(SO42-)、ρ(TDS)、总硬度和EC六项监测指标. 从地质构造角度考虑，平谷区第四系松散层孔隙水主要储存在冲洪积作用形成的含水层中，这些分布区域的地下水径流条件相对较好，含钙、镁化合物的岩层与地下水交互作用明显；同时，地下水在流经硫酸盐沉积物时发生溶滤过程，使得SO42-普遍存在，因此公因子F1代表了溶滤-富集作用.

公因子F2主要包括ρ(Na+)、ρ(NH4+)、ρ(Cl-)、ρ(Mn)和ρ(NO2)五项监测指标. 由于平谷区农业和养殖业密集程度较高，生产活动过程中含氮化肥、含氯农药和养殖废水下渗造成污染，因此Na+、Cl-、NH4+和NO2反映了地表农业和养殖业对地下水的影响. Mn指标受区域成土母质、土壤类型和人为活动影响较大[5]，但研究区内Mn主要分布在平谷区平原区西南部的前芮营垃圾填埋场以及平谷区东鹿角，并且向四周扩散，呈现出人为活动的表征. 因此，公因子F2代表了人为活动影响中的农业、养殖业和填埋场等人为活动产生的污染作用.

公因子F3主要包括的指标为ρ(CO2). 从水文地质化学上来看，浅层地下水中CO2可能来自大气降水渗透、地下微生物呼吸作用，以及含水层碳酸盐岩水文化学作用，但该研究中CO2的来源还需要结合SOM结果进行进一步分析.

表2 主要因子的特征值和累计方差贡献率

表3 旋转因子荷载矩阵

注：加粗体为荷载因子矩阵进行正交旋转后得分较高的指标.

图2 主要监测指标映射到SOM上的特征图谱Fig.2 Result of main components in SOM picture

图3 监测井映射到SOM上的特征图谱Fig.3 SOM pictures of all wells

公因子F4主要包括的指标为ρ(HCO3-)、ρ(F-)和ρ(As)，As主要集中分布在研究区西部的马昌营薄各庄村和柏店村，ρ(As)最高值为0.183 mg/L；F-的分布地点主要有研究区南西方向的马坊镇、马昌营薄各庄村以及前芮营村泃河附近，ρ(F-)整体上呈西部平原区大于东部近山区的特征，其最高值达4.60 mg/L，表明泃河可能已经遭受F-污染. 因此，公因子F4代表了人为活动影响中的工业污染作用.

公因子F5主要包括的指标为ρ〔Fe(Ⅱ)〕和浊度，表明浊度与Fe(Ⅱ)密切相关，反映了浊度主要来自Fe(Ⅱ). 公因子F6主要包括的指标为ρ〔Fe(Ⅲ)〕，反映了铁元素在水文地质化学中的变化情况. 公因子F7主要包括的指标为ρ(VP)，其主要存在于平谷区中心城区，表明人为活动中的工业源或生活源污染作用. 公因子F8主要包括的指标为ρ(Hg)，反映垃圾填埋场影响地下水质量.

2.3 污染物空间分布

图4 研究区域NH4+、NO3-和NO2污染物空间分布Fig.4 Spatial distribution of NH4+, NO3- and NO2 pollutants in the study area

通过SOM将全部指标以及监测井聚类结果映射到神经元上，根据神经元形成的特征图像比对指标间的关联性(见图2、3). 根据平谷区的水文地质情况映射到3个维度，分别代表了平谷区山前地带、过渡地带和平原地带[6]. 将特征图像以图片正中心为原点建立直角坐标系分析，神经元激活位置一致则说明指标正相关，反之则说明指标负相关. 由图2、3可见：①Na+响应区域位于第四象限，并且与平原地带特征图像吻合，表明了平原地带人为活动的影响；NH4+响应区域位于第三、四象限交界处，反映了NH4+兼受山前地带林业、种植业以及平原地带农业养殖业的双重影响；CO2分布与山前地带监测井特征图谱相似，均位于第二、第三象限，说明由于山前地带多为极弱富水区或弱富水区，降水过程中下渗极快，CO2组分不易逸失，据此判断PCA中的公因子F3代表了雨水下渗作用. ②NH4+、NO2、NO3-的响应区域分别位于不同象限，反映了3种污染物的空间分布特征不同，来源不一致. 平谷盆地山前地带含水层结构较为单一，以粗中砂夹卵砾石，砂、卵石夹中粗砂为主，而第四系沉积物主要分布在冲洪积平原，沉积厚度从数米至数百米不等，岩性为卵砾石夹粗砂. 在此条件下，NH4+作为阳离子较容易与土壤中阳离子发生交换，从而导致山前地带NH4+浓度低于平原地带；同时，山前地带水文地质条件也提供了良好的氧化氛围，影响了NO3-的分布(见图4). 从NH4+、NO2、NO3-三者的特征图谱响应强度来看，NO3-远超出NH4+和NO2，揭示NO3-污染可能是历史累积因素所致. Fe(Ⅱ)与Fe(Ⅲ)的响应区域也位于不同象限，Fe(Ⅱ) 更接近分布在山前地带，而Fe(Ⅲ)更接近分布在山前地带、平原地带，其原因也与山前地带比平原地带氧化氛围更好有关. ③从重金属分布上看，Cr6+、Cd、Al、Hg、Zn分布在山前地带，Cu、Mn、As分布在平原地带，Se在山前地带、平原地带和过渡地带均有分布，整体上保持自然背景水平[31]. 从有机物分布上，VP、CN分别分布在山前地带、过渡地带，由于此区域存在一定的村落和工业园区，因此有机物来源于人类活动；同时，Cd、Al、CN均在相同或相邻位置激活神经元，反映二者在一定程度上具有同源性.

3 结论

a) 平谷区地下水电导率与ρ(Ca2+)(p=0.936)、总碱度与ρ(HCO32-)(p=0.981)、ρ(Mg2+)与总硬度(p=0.944)指标之间呈显著相关，地下水化学类型主要以HCO3-Ca型为主，其次为HCO3-Mg型. 地下水中Ca2+、Mg2+、HCO3-均为自然起源. 平谷区地下水中NH4+、SO42-、Cd、Fe(Ⅱ)、NO2指标的空间分布离散性和差异性较大，说明上述指标存在局部富集现象.

b) PCA结果表明，影响平谷区地下水水质的8个公因子中，首要影响因子为溶滤-富集作用(贡献率为22.398%)、次要影响因子为农业、养殖业和填埋场等人为活动作用(贡献率为16.533%)，雨水下渗作用(贡献率为8.035%)、工业源人为活动(贡献率为7.466%)对地下水也有一定影响. As分布于平原地带，来源于地质环境背景因子. 西部平原区ρ(F-)大于东部近山区部分，反映该区域存在潜在F-污染.

c) SOM补充解释了PCA方法中非公因子的污染指标关联性和来源，反映了NH4+兼受山前地带林业、种植业和平原地带农业养殖业的双重影响；同时，NH4+、NO3-、NO2三者之间及Fe(Ⅱ)与Fe(Ⅲ)之间来源不同，Cd、Al、CN三者具有同一来源.