CN105468926A - 地下水型饮用水水源地污染源解析方法 - Google Patents

Abstract

本发明提供一种地下水型饮用水水源地污染源解析方法，包括：对采样数据进行预处理；计算预处理后的采样数据的特征值和特征向量；利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵；对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转；使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵。对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，避免矩阵分解的结果中出现负值，使得因子载荷和因子得分具有可解释性和明确的物理意义，同时确定污染源和贡献，而不需要事先知道源成分谱，运用简便，从而有效消除测量检测的不确定性。

Description

地下水型饮用水水源地污染源解析方法

技术领域

本发明涉及水体检测技术领域，特别是指一种地下水型饮用水水源地污染源解析方法。

背景技术

近年来，就全国范围而言，我国地下水质量下降趋势明显，不宜作为饮用水水源的比例逐渐增大，即使存在地下水水源地饮用水指标并不一定受到直接污染而超标现象，但其水质可能已经不同程度承受污染风险，且水质同时受地下水运动演化影响，因此考虑选择反映水质演化的地下水水化学类型常规指标，进行污染来源解析。地下水水质指标在取样与分析过程中，由于取样和样品处理、试剂和水纯度、仪器量度和仪器洁净、采用的分析方法、测定过程以及数据处理等过程均会产生测量误差(系统误差，随机误差，过失误差)。从取样到分析结果计算误差都绝对存在，虽然在各个过程中进行质量控制，但无法完全消除不确定性的影响。

发明内容

本发明要解决的技术问题是提供一种地下水型饮用水水源地污染源解析方法，能够有效消除测量检测的消除不确定性。

为解决上述技术问题，本发明的实施例提供一种地下水型饮用水水源地污染源解析方法，所述地下水型饮用水水源地污染源解析方法包括：

对采样数据进行预处理，所述预处理包括缺失值处理、不确定性处理和数据合理性分析中的至少一种；

计算预处理后的采样数据的特征值和特征向量，得到特征值矩阵和特征向量矩阵；

利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵；

对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转；

使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，所述因子载荷矩阵对应污染源指纹图谱，所述因子得分矩阵对应主要污染源贡献。

优选的，对采样数据进行缺失值处理包括：

丢弃采样数据中含缺失数据的记录；

用全局变量或属性的平均值代替缺失数据；

先根据欧式距离或相关分析来确定距离缺失数据样本最近的K个样本，将这K个值加权平均来估计该样本的缺失数据，其中k为大于零的整数；

采用神经网络模型预测每一个缺失数据；

对低于数据检测限的数据用数据检测限值或1/2检测限代替；

优选的，对采样数据进行不确定性处理包括：

利用下式计算采样数据的不确定性：

$U=\sqrt{{(s\×c)}^2+l^2}$

其中，s为误差百分数；c为指标浓度值；l为因子数据检出限。

优选的，对采样数据进行数据合理性分析，包括：

采用下式对采样数据进行筛选：

$S/N=\sqrt{\frac{\underset{i=1}{\overset{n}{\Σ}}{(x_{ij}-s_{ij})}^2}{\underset{i=1}{\overset{n}{\Σ}}{s_{ij}}^2}}$

其中，S/N为信噪比，x_ij表示第i采样点第j个样品的浓度，s_ij表示第i采样点第j个样品的标准偏差。

优选的，所述计算预处理后采样数据的特征值和特征向量，得到特征值矩阵和特征向量矩阵，包括：

对预处理后的采样数据进行无量纲化；

计算无量纲化后采样数据的协方差矩阵；

利用雅各布方法计算协方差矩阵的特征值和特征向量，得到特征值矩阵和特征向量矩阵；

优选的，所述计算无量纲化后采样数据的协方差矩阵，包括：

利用下式计算无量纲化后采样数据的协方差矩阵：

Z＝DD′

其中，Z为协方差矩阵，D为无量纲化后的样品数据矩阵，D′为D的转置；

利用雅各布方法计算协方差矩阵的特征值和特征向量，得到特征值矩阵和特征向量矩阵，协方差矩阵表示为：

Z＝QEQ′

其中，E为特征值矩阵，Q为特征向量矩阵，Q′为Q的转置。

优选的，所述利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵之前，包括：

利用累计方差贡献率提取显著性因子；

所述利用累计方差贡献率提取显著性因子的提取条件为：

$\mathrm{CV}=\frac{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_j}{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j}\×100\%$ ＞第一限定值

其中：n为显著性因子个数；m为污染物个数；λ为特征值。

优选的，所述利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵，所述因子载荷矩阵表示为：

S＝QE^1/2

因子得分矩阵可表示为：

C＝(S′S)^-1S′D

其中,S为因子载荷矩阵，C为因子得分矩阵，S’为S的转置。

优选的，所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，包括：

利用下式对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转：

S₁＝ST₁ ^-1

C₁＝T₁C

其中，S₁为旋转后的因子载荷矩阵；C₁为旋转后的因子得分矩阵；T₁为转换矩阵，且T₁＝(C^*C′)(C′C^*)^-1，C^*为把C中的负值替换为零后的因子得分矩阵；

利用因子载荷矩阵S₁构建转换矩阵T₂对因子载荷矩阵S₁和因子得分矩阵C₁继续旋转：

S₂＝S₁T₂

C₂＝T₂ ^-1C₁

其中：S₂为二次旋转后的因子载荷矩阵；C₂为二次旋转后的因子得分矩阵；T₂为二次转换矩阵，且T₂＝(S₁′S₁)^-1(S₁′S₁ ^*)，S₁ ^*为S₁中的负值换为零后的因子载荷矩阵。

优选的，所述使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，包括：

重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，因子载荷中负值的平方和小于预先设定的第二限定值时，输出因子载荷矩阵和因子得分矩阵。

本发明的上述技术方案的有益效果如下：

上述方案中，对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，避免矩阵分解的结果中出现负值，使得因子载荷和因子得分具有可解释性和明确的物理意义，使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，能够同时确定污染源和贡献，而不需要事先知道源成分谱，运用简便，从而有效消除测量检测的不确定性。

附图说明

图1为本发明实施例的地下水型饮用水水源地污染源解析方法流程图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

如图1所示，本发明的实施例一种地下水型饮用水水源地污染源解析方法，所述地下水型饮用水水源地污染源解析方法包括：

步骤101：对采样数据进行预处理，所述预处理包括缺失值处理、不确定性处理和数据合理性分析中的至少一种。

步骤102：计算预处理后的采样数据的特征值和特征向量，得到特征值矩阵和特征向量矩阵。

步骤103：利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵。

步骤104：对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转。

步骤105：使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，所述因子载荷矩阵对应污染源指纹图谱，所述因子得分矩阵对应主要污染源贡献。

本发明实施例的地下水型饮用水水源地污染源解析方法，对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，避免矩阵分解的结果中出现负值，使得因子载荷和因子得分具有可解释性和明确的物理意义，使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，能够同时确定污染源和贡献，而不需要事先知道源成分谱，运用简便，从而有效消除测量检测的不确定性。

正定矩阵因子分析是基于多元统计的分析方法，对数据有效性具有一定的要求，因此在进行分析之前首先对数据进行预处理。

优选的，对采样数据进行缺失值处理包括：

丢弃采样数据中含缺失数据的记录；

用全局变量或属性的平均值代替缺失数据；

先根据欧式距离或相关分析来确定距离缺失数据样本最近的K个样本，将这K个值加权平均来估计该样本的缺失数据，其中k为大于零的整数；

采用神经网络模型预测每一个缺失数据；

对低于数据检测限的数据用数据检测限值或1/2检测限代替；

优选的，对采样数据进行不确定性处理包括：

利用下式计算采样数据的不确定性：

$U=\sqrt{{(s\×c)}^2+l^2}$

其中，s为误差百分数；c为指标浓度值；l为因子数据检出限。

优选的，对采样数据进行数据合理性分析，包括：

采用下式对采样数据进行筛选：

$S/N=\sqrt{\frac{\underset{i=1}{\overset{n}{\Σ}}{(x_{ij}-s_{ij})}^2}{\underset{i=1}{\overset{n}{\Σ}}{s_{ij}}^2}}$

其中，S/N为信噪比，x_ij表示第i采样点第j个样品的浓度，s_ij表示第i采样点第j个样品的标准偏差。

信噪比小，说明样品的噪声大，信噪比越大则表示样品检出的可能性越大，越适合模型。

优选的，所述计算预处理后采样数据的特征值和特征向量，得到特征值矩阵和特征向量矩阵，包括：

对预处理后的采样数据进行无量纲化；

计算无量纲化后采样数据的协方差矩阵；

利用雅各布方法计算协方差矩阵的特征值和特征向量，得到特征值矩阵和特征向量矩阵；

优选的，所述计算无量纲化后采样数据的协方差矩阵，包括：

利用下式计算无量纲化后采样数据的协方差矩阵：

Z＝DD′

其中，Z为协方差矩阵，D为无量纲化后的样品数据矩阵，D′为D的转置；

利用雅各布方法计算协方差矩阵的特征值和特征向量，得到特征值矩阵和特征向量矩阵，协方差矩阵表示为：

Z＝QEQ′

其中，E为特征值矩阵，Q为特征向量矩阵，Q′为Q的转置。

优选的，所述利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵之前，包括：

利用累计方差贡献率提取显著性因子；

所述利用累计方差贡献率提取显著性因子的提取条件为：

$\mathrm{CV}=\frac{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_j}{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j}\×100\%$ ＞第一限定值

其中：n为显著性因子个数；m为污染物个数；λ为特征值。

其中，第一限定值可以设置为85％。

本实施例中，利用累计方差贡献率提取显著性因子，能够使高维变量空间降维后能尽可能保留原来指标信息。

优选的，所述利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵，所述因子载荷矩阵表示为：

S＝QE^1/2

因子得分矩阵可表示为：

C＝(S′S)^-1S′D

其中,S为因子载荷矩阵，C为因子得分矩阵，S’为S的转置。

优选的，所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，包括：

利用下式对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转：

S₁＝ST₁ ^-1

C₁＝T₁C

其中，S₁为旋转后的因子载荷矩阵；C₁为旋转后的因子得分矩阵；T₁为转换矩阵，且T₁＝(C^*C′)(C′C^*)^-1，C^*为把C中的负值替换为零后的因子得分矩阵；

利用因子载荷矩阵S₁构建转换矩阵T₂对因子载荷矩阵S₁和因子得分矩阵C₁继续旋转：

S₂＝S₁T₂

C₂＝T₂ ^-1C₁

其中：S₂为二次旋转后的因子载荷矩阵；C₂为二次旋转后的因子得分矩阵；T₂为二次转换矩阵，且T₂＝(S₁′S₁)^-1(S₁′S₁ ^*)，S₁ ^*为S₁中的负值换为零后的因子载荷矩阵。

本实施例中，考虑了原始数据的不确定性，对数据偏差进行了校正，避免矩阵分解的结果中出现负值，使得因子载荷和因子得分具有可解释性和明确的物理意义和科学的解释。

优选的，所述使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，包括：

重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，因子载荷中负值的平方和小于预先设定的第二限定值时，输出因子载荷矩阵和因子得分矩阵。

其中，本领域技术人员可以根据需要设置第二限定值。

本实施例中，用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，得到的污染源不需要转换就可以直接与原始数据矩阵作比较，PMF方法能够同时确定污染源和贡献，而不需要事先知道源成分谱。适用于水文地质条件简单，观测数据量较大，污染源和污染种类相对较少的地区，运用简便，可应用分析软件进行计算。

模型基本方程可以为：

X_nm＝G_npF_pm+E

其中，n为取样点数，m为各取样点测试的成分数量，p为污染源个数；X_nm为取样点各成分含量；G_np为主要源的贡献率；F_pm为源指纹图谱。

上述方案中，对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，避免矩阵分解的结果中出现负值，使得因子载荷和因子得分具有可解释性和明确的物理意义，使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，能够同时确定污染源和贡献，而不需要事先知道源成分谱，运用简便，从而有效消除测量检测的不确定性。

本次选取反映地下水水化学类型的11项指标(Cl^-、SO₄ ^2-、NO₃ ^-、HCO₃ ^-、Mg²⁺、Ca²⁺、Na⁺、EC、TDS、T-Hard、NH₄-N)进行因子分析。为提高数据集质量，首先对数据缺失值进行处理、shapiro-wilk正态性检验、Box-Cox变换。

2)不确定性计算

考虑浓度值15％的数据误差，s为误差百分数，c为指标浓度值，l为因子数据检出限，计算数据不确定性

3)根据公式计算信噪比

本研究所用数据在放入模型前以信噪比S/N(signaltonoise)作为标准进行筛选，信噪比S/N为:式中，x_ij表示第i采样点第j个样品的浓度，s_ij表示第i采样点第j个样品的标准偏差。

4)模型的输入

根据PMF软件对数据格式的要求，将数据以矩阵的形式输入模型。

5)模型运行与结果

通过修改模型迭代次数和分解因子数，对模型残差进行分析，选择误差最小的参数组合进行正定矩阵因子分解，运行模型，查看运行结果。

本次模拟将水质影响因素划分为3个主要指标，因子F1以氨氮为主，F1反映地下水受到工业或农业污染。主因子F2中以Cl^-、SO₄ ^2-、Mg²⁺、NO₃ ^-、Ca²⁺为主，溶滤作用；因子F3以电导率、TDS、总硬度、Na⁺、HCO₃ ^-为主，受蒸发浓缩作用影响。

综上所述，根据以上基于多元统计分析的方法进行水质指标及因子的分类，分别受溶滤作用、蒸发浓缩作用以及人为活动影响。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (10)

1.一种地下水型饮用水水源地污染源解析方法，其特征在于，所述地下水型饮用水水源地污染源解析方法包括：

对采样数据进行预处理，所述预处理包括缺失值处理、不确定性处理和数据合理性分析中的至少一种；

计算预处理后的采样数据的特征值和特征向量，得到特征值矩阵和特征向量矩阵；

利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵；

对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转；

使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，所述因子载荷矩阵对应污染源指纹图谱，所述因子得分矩阵对应主要污染源贡献。

2.根据权利要求1所述的地下水型饮用水水源地污染源解析方法，其特征在于，对采样数据进行缺失值处理包括：

丢弃采样数据中含缺失数据的记录；

用全局变量或属性的平均值代替缺失数据；

先根据欧式距离或相关分析来确定距离缺失数据样本最近的K个样本，将这K个值加权平均来估计该样本的缺失数据，其中k为大于零的整数；

采用神经网络模型预测每一个缺失数据；

对低于数据检测限的数据用数据检测限值或1/2检测限代替。

3.根据权利要求1所述的地下水型饮用水水源地污染源解析方法，其特征在于，对采样数据进行不确定性处理包括：

利用下式计算采样数据的不确定性：

其中，s为误差百分数；c为指标浓度值；l为因子数据检出限。

4.根据权利要求1所述的地下水型饮用水水源地污染源解析方法，其特征在于，对采样数据进行数据合理性分析，包括：

采用下式对采样数据进行筛选：

其中，S/N为信噪比，x_ij表示第i采样点第j个样品的浓度，s_ij表示第i采样点第j个样品的标准偏差。

5.根据权利要求1所述的地下水型饮用水水源地污染源解析方法，其特征在于，所述计算预处理后采样数据的特征值和特征向量，得到特征值矩阵和特征向量矩阵，包括：

对预处理后的采样数据进行无量纲化；

计算无量纲化后采样数据的协方差矩阵；

利用雅各布方法计算协方差矩阵的特征值和特征向量，得到特征值矩阵和特征向量矩阵。

6.根据权利要求5所述的地下水型饮用水水源地污染源解析方法，其特征在于，所述计算无量纲化后采样数据的协方差矩阵，包括：

利用下式计算无量纲化后采样数据的协方差矩阵：

Z＝DD′

其中，Z为协方差矩阵，D为无量纲化后的样品数据矩阵，D′为D的转置；

利用雅各布方法计算协方差矩阵的特征值和特征向量，得到特征值矩阵和特征向量矩阵，协方差矩阵表示为：

Z＝QEQ′

其中，E为特征值矩阵，Q为特征向量矩阵，Q′为Q的转置。

7.根据权利要求1所述的地下水型饮用水水源地污染源解析方法，其特征在于，所述利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵之前，包括：

利用累计方差贡献率提取显著性因子；

所述利用累计方差贡献率提取显著性因子的提取条件为：

其中：n为显著性因子个数；m为污染物个数；λ为特征值。

8.根据权利要求6所述的地下水型饮用水水源地污染源解析方法，其特征在于，所述利用所述特征值矩阵和特征向量矩阵求得因子载荷矩阵和因子得分矩阵，所述因子载荷矩阵表示为：

S＝QE^1/2

因子得分矩阵可表示为：

C＝(S′S)^-1S′D

其中,S为因子载荷矩阵，C为因子得分矩阵，S’为S的转置。

9.根据权利要求8所述的地下水型饮用水水源地污染源解析方法，其特征在于，所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，包括：

利用下式对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转：

S₁＝ST₁ ^-1

C₁＝T₁C

其中，S₁为旋转后的因子载荷矩阵；C₁为旋转后的因子得分矩阵；T₁为转换矩阵，且T₁＝(C^*C′)(C′C^*)^-1，C^*为把C中的负值替换为零后的因子得分矩阵；

利用因子载荷矩阵S₁构建转换矩阵T₂对因子载荷矩阵S₁和因子得分矩阵C₁继续旋转：

S₂＝S₁T₂

C₂＝T₂ ^-1C₁

其中：S₂为二次旋转后的因子载荷矩阵；C₂为二次旋转后的因子得分矩阵；T₂为二次转换矩阵，且T₂＝(S₁′S₁)^-1(S₁′S₁ ^*)，S₁ ^*为S₁中的负值换为零后的因子载荷矩阵。

10.根据权利要求9所述的地下水型饮用水水源地污染源解析方法，其特征在于，所述使用最小二乘方法进行迭代运算重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，确定因子载荷矩阵和因子得分矩阵，包括：

重复所述对所述因子载荷矩阵和因子得分矩阵进行非负约束的旋转，因子载荷中负值的平方和小于预先设定的第二限定值时，输出因子载荷矩阵和因子得分矩阵。