李德安 邓一荣 王俊 黄灶泉 谷培科

(广东省环境科学研究院,广东 广州 510045)

土壤是经济社会可持续发展的物质基础[1]，而土壤重金属污染问题正在引发越来越多的关注。我国于2005年4月—2013年12月开展了首次全国土壤污染状况调查，调查结果显示，我国土壤镉、汞、砷、铜、铅、铬、锌、镍8种无机污染物点位超标率分别为7.0%、1.6%、2.7%、2.1%、1.5%、1.1%、0.9%、4.8%，全国土壤环境状况总体不容乐观[2]。在掌握土壤环境质量状况的基础上，有必要进一步开展土壤污染溯源工作，为控源断源精准施策提供支撑。我国在“十三五”期间已在若干省份开展耕地土壤污染成因排查分析试点，预计该项工作将成为我国土壤污染防治“十四五”规划的重要工作进一步铺开。在已发布的浙江省《土壤、地下水和农业农村污染防治“十四五”规划》中，土壤污染溯源被列为重点任务[3]。土壤重金属的来源往往包括自然源和多种人为源等，溯源工作较为复杂[4]。现有的土壤污染源解析方法大致包括多元统计分析法[5-7]、同位素分析法[8,9]以及受体模型法[10,11]等。其中受体模型法因不需要分析污染物的迁移转化路径等优点得到了广泛应用[12]。常用的受体模型法包括绝对主成分得分-多元线性回归法(APCS-MLR)[13,14]、UNMIX模型[15]以及正定矩阵因子分析法(PMF)[16,17]等。

正定矩阵因子分析法(PMF)最早应用于大气污染物源解析，该方法利用用户提供的不确定度作为各样点浓度数据的权重，然后根据最小二乘法进行迭代运算。迭代计算结果可用于定量确定各污染源的成分谱(FactorProfile)及贡献率(FactorContributions)[12]。目前土壤污染源解析领域利用PMF模型的主要目的在于确定污染源的贡献占比[18]。由于无法确切得知土壤在成土过程及人为输入影响中重金属的真实来源及贡献程度，因此试图对模型拟合所得的贡献占比的真假或优劣作出评价是相当困难的。目前对模型结果可信度的判断主要借助软件自带的不确定度估计功能[19-21]或与其它受体模型结果进行交叉对比[10,16,22-24]等手段。但仍然缺乏模型结果与真实结果的对比验证。

本文尝试利用自行编制设计的数据集作为PMF模型的原始输入，在已知数据集的源成分谱及贡献率的条件下验证PMF模型源解析结果的准确性，以便更好地理解PMF模型的工作原理以及拟合结果的现实意义，为土壤重金属污染源解析工作下阶段的进一步铺开完善方法论。

1 研究方法

1.1 数据设计

PMF等受体模型将样品浓度数据集看作i行j列的矩阵(Xij)，其中i为样品数量，j为重金属的数量，本研究样品数量i设置为100，重金属数量设定为Cd、Hg、As、Pb、Cr、Cu、Zn、Ni共8种。这些重金属仅作为元素标识，不代表实际意义。而PMF源解析的目标在于将样品浓度数据矩阵Xij分解成源贡献率矩阵Gik和源成分谱矩阵Fkj，以及一个残差矩阵Eij，如公式(1)所示[24]：

Xij=GikFkj+Eij

(1)

式中，Xij为i个样品j个重金属数据组成的i行j列矩阵；Gik为i个样品k个源(或称因子)的贡献数据组成的i行k列矩阵，代表某个源在某个样品中的贡献大小；Fkj为k个源j个重金属的浓度数据组成的k行j列矩阵，代表某个源中某种重金属的含量高低；Eij为i行j列的残差矩阵，代表某个样品中未被模型解释的某种重金属含量的大小。PMF模型通过多次迭代计算获得最佳分解矩阵Gik和Fkj，使得目标函数Q最小，Q的定义如公式(2)所示：

(2)

式中，en,m为残差矩阵Eij中第n行第m列中的残差数值；μn,m为与浓度数据矩阵Xij相对应的不确定度矩阵μij中第n行第m列中的不确定度数值。不确定度是PMF模型所要求的输入数据集之一，代表原始观测数据集Xij由于采样和分析测试等误差所带来的不确定度，可由测试实验室提供或用户根据检出限估算。根据检出限估算不确定度如公式(3)所示：

(3)

式中，Urel为相对误差，本研究取10%；xn,m为浓度数据矩阵Xij中第n行第m列中的浓度数值；MDL为检出限，本研究取0.01，远小于任一浓度数据xn,m值。

本研究基于公式(1)进行浓度数据矩阵Xij的反向构造，即通过人为设定Gik和Fkj矩阵然后将其相乘得到浓度数据矩阵Xij；基于公式(3)或生成随机数的方式构造不确定度矩阵μij，共设计3种编制情景，每种情景的特征汇总见表1。所有情景样品数i均为100，元素数j均为8，源数量k均为3。

1.2 数据处理与分析

本研究利用Microsoft Excel 2019进行Gik中随机数的生成、Gik和Fkj矩阵的乘法运算、μij中随机数的生成等；利用EPA PMF 5.0进行PMF源解析；利用Origin 2018进行图表绘制。

2 结果与讨论

2.1 源数量设定对源解析结果的影响

利用编制情景1中的浓度与不确定度数据进行PMF源解析，源数量分别设定为2、3和4。不同源数量设定下模型对编制数据的拟合程度如图1所示。从数据编制过程中源数量k值的设定可知(表1中5种编制情景k值均为3)，最符合编制数据的源数量应为3。而从图1可以看出，在PMF模型中源数量设定为2的情况下，模型对编制数据的拟合度不如3源和4源的设定情景，尤其是Cu、Cr、As 3种元素，相关系数r2值仅分别为0.20、0.62和0.64，而3源、4源条件下8种元素的r2值均达到1.0。三者目标函数Q值的差异也较大，2源、3源和4源的Qture值分别为42855.6、0.00002和0.0002，其中2源的Qture值远大于理论值Q值(样品数i与重金属类型数量j的乘积[16]，本研究中为100×8=800)，而3源和4源均接近于0。

图1 不同源数量设定下PMF模型对编制数据的拟合度对比

不同源数量的模拟结果表明，当设定源数量小于最优值时，数据拟合度随设定源数量增加而显著提高，Q值则随设定源数量增加而显著下降。当设定源数量等于最优值时，继续增加源数量模型拟合度和Q值呈现基本稳定的趋势，因此源数量最优值可看作模型拟合度和Q值变化曲线的拐点。该结论与文献报道的源数量选取方式基本一致[25]。

2.2 源解析矩阵与编制矩阵的对应关系

图1b显示了模型对编制数据Xij较为理想的拟合程度，但仍有必要进一步考察模型拟合所得的源贡献率矩阵Gik和源成分谱矩阵Fkj与编制矩阵的对应关系，因为在污染源解析领域，这2个分解矩阵更具有实际应用价值。

对源解析文献中常用的源成分谱百分比数据进行考察，图2中a为编制情景1实际所构造的Fkj矩阵，b则为PMF模型拟合所得的Fkj矩阵，二者均转化为百分比形式。表2则分别列出了二者未转化为百分比的原始数据。

可以发现，图2中的百分比堆积柱状图a与b一致程度较高，但表2中非百分比的原始数据a与b存在数量级级别的差异。进一步查询PMF拟合结果中的Gik矩阵，发现软件对该矩阵进行了归一化处理(即Gik矩阵任一列向量的均值均为1)。但PMF模型分解出的Gik与Fkj矩阵相乘所得的Xij仍然与编制矩阵一致，见图1b，说明在PMF拟合过程中，Gik任一列向量的均值被转移至Fkj矩阵的对应行向量中。该过程用矩阵形式表示如下(为便于区别，上标带*号的表示编制矩阵，不带*号的表示PMF拟合矩阵，Ave_G表示Gik某一列向量的均值)：

(4)

将编制矩阵Fkj*的行向量乘上编制矩阵Gik*列向量均值(简记为Ave_G×Fkj*)，原始数值列于表2，转化为百分比形式绘制于图2c，可以发现该原始数值及百分比与PMF拟合数据均有较高的一致性。而图2a(即Fkj*)与图2b(即Fkj)百分比占比较接近,可能是源于编制情景1中Gik*列向量均值较为接近，因其所使用的生成随机数范围均为1～100，参见表1情景1，均值都接近50，参见表2Gik*列向量均值。这一推测可从编制情景2的拟合结果中得到验证。情景2改变了Gik*的构造方式，使其列向量均值产生明显大小差异，参见表1情景2。PMF拟合的Fkj原始数据与百分比形式分别列于表3和图3。可以看出,该情景下不管是原始数据还是百分比占比，Fkj都与编制Fkj*差异较大，而与Ave_G×Fkj*相接近。

表2 编制情景1源成分谱对比

图2 编制情景1源成分谱百分比对比度对比(F1～F3分别指代3个源Factor1～Factor3) 图3 编制情景2源成分谱百分比对比度对比

实际上这种处理方式是有现实意义的。当前利用PMF进行源解析的研究中，均把源成分谱Fkj的百分比占比理解为源贡献率[16,17,18-20,22-24]，如图3b PMF源成分谱拟合结果中Cd元素的F1占比为58%，则意味着F1所代表的污染源对Cd的平均贡献率为58%。在明确了PMF拟合结果中的源成分谱数据Fkj实际上为Ave_G×Fkj*后，可以证明源成分谱Fkj的百分比占比的物理含义就是源贡献率。基于编制数据尝试计算F1对样品1中Cd含量的贡献率，样品1中Cd的总浓度为x11*(第1个下标代表样品编号i=1；第2个下标代表重金属元素编号j=1)，而根据公式(1)或(4)可知x11*的计算遵循公式(5)。而根据源贡献率的定义，F1对样品1中Cd含量的贡献率P(F1)的计算遵循公式(6)。

(5)

(7)

推而广之，F1对100个样品中Cd含量的平均源贡献率P(F1)ave的计算应遵循公式(7)。从公式(7)可知，P(F1)ave实际上等于PMF源成分谱Fkj的百分比占比。从公式(7)还可以发现，Fkj列向量的加和即为所有样品中第j个元素的平均值，因此从表3的2行合计项中可以发现二者几乎一致。

表3 编制情景2源成分谱对比

2.3 观测误差对源解析结果的影响

在编制情景1和2中，浓度数据矩阵Xij*直接等于Gik*与Fkj*的乘积，并未引入观测误差。而在实际应用中误差是不可避免的。因此在第3种编制情景中，利用Excel生成随机数的方式为浓度数据矩阵Xij*引入随机误差，详见表1编制情景3，并将该误差绝对值的2倍作为输入PMF模型的不确定度数据矩阵ij*，考察带有误差条件下PMF的拟合情况。

情景3与情景2的样品浓度与不确定度的对比如图4所示。在情景2中，由于不确定度数据μij*均按照公式(3)进行设计，因此每个样品的不确定度均与浓度成比例，且相应的信噪比S/N均高达10.0。而情景3引入随机误差后，信噪比明显下降，更符合实际应用情况。但信噪比S/N仍满足大于1.0的条件，表明样品浓度数据的信号强度满足模型拟合要求。

图4 样品浓度与不确定度对比

将情景3的模型拟合结果与该情景编制浓度数据进行对比，结果如图5所示。在图1b中，对于未引入误差的情景1，PMF模型取得了很好的拟合效果，3源情况下各元素的r2值均达到1.0。引入随机误差后，PMF模型对情景3的拟合优度有所下降，但仍达到了相对较高的水平，各元素的r2值在0.87～0.99，对于实际应用而言是可以接受的。因误差的引入，PMF拟合Qture值相比情景2也有一定增长，为118.02，但仍远小于理论Q值800，表明模型拟合结果较好。

图5 PMF模型对情景3编制数据的拟合情况 图6 编制情景3源成分谱百分比对比度对比

进一步考察在土壤重金属PMF源解析领域中最为关注的源贡献率的情况，源贡献率的百分比占比和原始数值分别参见图6和表4。从图6可以看出，尽管因为观测误差的存在导致PMF模型对原始数据的拟合度有一定下降，但从源成分谱的百分比占比来看，PMF拟合结果与原始编制数据所表征的源贡献率是基本一致的。拟合所得的各元素的主导贡献源均未发生偏离。而从表4可以看出，与表3未引入误差的情景2类似，PMF拟合Fkj的合计项与Ave_G×Fkj*的合计项仍保持高度一致。这表明在原始浓度数据带有误差的条件下，PMF模型仍能取得较为理想的拟合结果以及源贡献率信息。

表4 编制情景2源成分谱对比

3 结论

本文利用自行编制设计的数据集作为PMF模型的原始输入，在已知数据集的源成分谱及贡献率的条件下验证了PMF模型的源解析结果，主要结论如下。

在PMF模型拟合过程中需要设定源数量，源数量最优值可近似按照模型拟合度和Q值变化曲线的拐点确定。当设定源数量小于最优值时，数据拟合度随设定源数量增加而显著提高，Q值则随设定源数量增加而显著下降。当设定源数量等于最优值时，继续增加源数量模型拟合度和Q值呈现基本稳定的趋势。

在不带观测误差和带有观测误差2种条件下，PMF模型均能较好地拟合原始数据，其中不带误差条件下拟合r2值可达1.0，带有误差条件下拟合度有所下降，但仍高于0.87。

PMF拟合所得的源成分谱Fkj实际为编制源贡献矩阵列向量均值与编制源成分谱矩阵行向量的乘积，其物理含义即为某个源对某种重金属的平均贡献含量，因此源成分谱的百分比占比的物理含义即为源贡献率。在不带观测误差和带有观测误差2种条件下，PMF拟合所得源贡献率与原始编制数据所表征的源贡献率一致性较好，表明在满足模型假设的条件下，PMF用于土壤重金属源解析可取得较为理想的结果。