基于編制數(shù)據(jù)的土壤重金屬PMF源解析研究

李德安 鄧一榮 王俊 黃灶泉 谷培科

(廣東省環(huán)境科學(xué)研究院,廣東 廣州 510045)

土壤是經(jīng)濟(jì)社會可持續(xù)發(fā)展的物質(zhì)基礎(chǔ)[1]，而土壤重金屬污染問題正在引發(fā)越來越多的關(guān)注。我國于2005年4月—2013年12月開展了首次全國土壤污染狀況調(diào)查，調(diào)查結(jié)果顯示，我國土壤鎘、汞、砷、銅、鉛、鉻、鋅、鎳8種無機(jī)污染物點(diǎn)位超標(biāo)率分別為7.0%、1.6%、2.7%、2.1%、1.5%、1.1%、0.9%、4.8%，全國土壤環(huán)境狀況總體不容樂觀[2]。在掌握土壤環(huán)境質(zhì)量狀況的基礎(chǔ)上，有必要進(jìn)一步開展土壤污染溯源工作，為控源斷源精準(zhǔn)施策提供支撐。我國在“十三五”期間已在若干省份開展耕地土壤污染成因排查分析試點(diǎn)，預(yù)計該項工作將成為我國土壤污染防治“十四五”規(guī)劃的重要工作進(jìn)一步鋪開。在已發(fā)布的浙江省《土壤、地下水和農(nóng)業(yè)農(nóng)村污染防治“十四五”規(guī)劃》中，土壤污染溯源被列為重點(diǎn)任務(wù)[3]。土壤重金屬的來源往往包括自然源和多種人為源等，溯源工作較為復(fù)雜[4]?，F(xiàn)有的土壤污染源解析方法大致包括多元統(tǒng)計分析法[5-7]、同位素分析法[8,9]以及受體模型法[10,11]等。其中受體模型法因不需要分析污染物的遷移轉(zhuǎn)化路徑等優(yōu)點(diǎn)得到了廣泛應(yīng)用[12]。常用的受體模型法包括絕對主成分得分-多元線性回歸法(APCS-MLR)[13,14]、UNMIX模型[15]以及正定矩陣因子分析法(PMF)[16,17]等。

正定矩陣因子分析法(PMF)最早應(yīng)用于大氣污染物源解析，該方法利用用戶提供的不確定度作為各樣點(diǎn)濃度數(shù)據(jù)的權(quán)重，然后根據(jù)最小二乘法進(jìn)行迭代運(yùn)算。迭代計算結(jié)果可用于定量確定各污染源的成分譜(FactorProfile)及貢獻(xiàn)率(FactorContributions)[12]。目前土壤污染源解析領(lǐng)域利用PMF模型的主要目的在于確定污染源的貢獻(xiàn)占比[18]。由于無法確切得知土壤在成土過程及人為輸入影響中重金屬的真實來源及貢獻(xiàn)程度，因此試圖對模型擬合所得的貢獻(xiàn)占比的真假或優(yōu)劣作出評價是相當(dāng)困難的。目前對模型結(jié)果可信度的判斷主要借助軟件自帶的不確定度估計功能[19-21]或與其它受體模型結(jié)果進(jìn)行交叉對比[10,16,22-24]等手段。但仍然缺乏模型結(jié)果與真實結(jié)果的對比驗證。

本文嘗試?yán)米孕芯幹圃O(shè)計的數(shù)據(jù)集作為PMF模型的原始輸入，在已知數(shù)據(jù)集的源成分譜及貢獻(xiàn)率的條件下驗證PMF模型源解析結(jié)果的準(zhǔn)確性，以便更好地理解PMF模型的工作原理以及擬合結(jié)果的現(xiàn)實意義，為土壤重金屬污染源解析工作下階段的進(jìn)一步鋪開完善方法論。

1 研究方法

1.1 數(shù)據(jù)設(shè)計

PMF等受體模型將樣品濃度數(shù)據(jù)集看作i行j列的矩陣(Xij)，其中i為樣品數(shù)量，j為重金屬的數(shù)量，本研究樣品數(shù)量i設(shè)置為100，重金屬數(shù)量設(shè)定為Cd、Hg、As、Pb、Cr、Cu、Zn、Ni共8種。這些重金屬僅作為元素標(biāo)識，不代表實際意義。而PMF源解析的目標(biāo)在于將樣品濃度數(shù)據(jù)矩陣Xij分解成源貢獻(xiàn)率矩陣Gik和源成分譜矩陣Fkj，以及一個殘差矩陣Eij，如公式(1)所示[24]：

Xij=GikFkj+Eij

(1)

式中，Xij為i個樣品j個重金屬數(shù)據(jù)組成的i行j列矩陣；Gik為i個樣品k個源(或稱因子)的貢獻(xiàn)數(shù)據(jù)組成的i行k列矩陣，代表某個源在某個樣品中的貢獻(xiàn)大?。籉kj為k個源j個重金屬的濃度數(shù)據(jù)組成的k行j列矩陣，代表某個源中某種重金屬的含量高低；Eij為i行j列的殘差矩陣，代表某個樣品中未被模型解釋的某種重金屬含量的大小。PMF模型通過多次迭代計算獲得最佳分解矩陣Gik和Fkj，使得目標(biāo)函數(shù)Q最小，Q的定義如公式(2)所示：

(2)

式中，en,m為殘差矩陣Eij中第n行第m列中的殘差數(shù)值；μn,m為與濃度數(shù)據(jù)矩陣Xij相對應(yīng)的不確定度矩陣μij中第n行第m列中的不確定度數(shù)值。不確定度是PMF模型所要求的輸入數(shù)據(jù)集之一，代表原始觀測數(shù)據(jù)集Xij由于采樣和分析測試等誤差所帶來的不確定度，可由測試實驗室提供或用戶根據(jù)檢出限估算。根據(jù)檢出限估算不確定度如公式(3)所示：

(3)

式中，Urel為相對誤差，本研究取10%；xn,m為濃度數(shù)據(jù)矩陣Xij中第n行第m列中的濃度數(shù)值；MDL為檢出限，本研究取0.01，遠(yuǎn)小于任一濃度數(shù)據(jù)xn,m值。

本研究基于公式(1)進(jìn)行濃度數(shù)據(jù)矩陣Xij的反向構(gòu)造，即通過人為設(shè)定Gik和Fkj矩陣然后將其相乘得到濃度數(shù)據(jù)矩陣Xij；基于公式(3)或生成隨機(jī)數(shù)的方式構(gòu)造不確定度矩陣μij，共設(shè)計3種編制情景，每種情景的特征匯總見表1。所有情景樣品數(shù)i均為100，元素數(shù)j均為8，源數(shù)量k均為3。

1.2 數(shù)據(jù)處理與分析

本研究利用Microsoft Excel 2019進(jìn)行Gik中隨機(jī)數(shù)的生成、Gik和Fkj矩陣的乘法運(yùn)算、μij中隨機(jī)數(shù)的生成等；利用EPA PMF 5.0進(jìn)行PMF源解析；利用Origin 2018進(jìn)行圖表繪制。

2 結(jié)果與討論

2.1 源數(shù)量設(shè)定對源解析結(jié)果的影響

利用編制情景1中的濃度與不確定度數(shù)據(jù)進(jìn)行PMF源解析，源數(shù)量分別設(shè)定為2、3和4。不同源數(shù)量設(shè)定下模型對編制數(shù)據(jù)的擬合程度如圖1所示。從數(shù)據(jù)編制過程中源數(shù)量k值的設(shè)定可知(表1中5種編制情景k值均為3)，最符合編制數(shù)據(jù)的源數(shù)量應(yīng)為3。而從圖1可以看出，在PMF模型中源數(shù)量設(shè)定為2的情況下，模型對編制數(shù)據(jù)的擬合度不如3源和4源的設(shè)定情景，尤其是Cu、Cr、As 3種元素，相關(guān)系數(shù)r2值僅分別為0.20、0.62和0.64，而3源、4源條件下8種元素的r2值均達(dá)到1.0。三者目標(biāo)函數(shù)Q值的差異也較大，2源、3源和4源的Qture值分別為42855.6、0.00002和0.0002，其中2源的Qture值遠(yuǎn)大于理論值Q值(樣品數(shù)i與重金屬類型數(shù)量j的乘積[16]，本研究中為100×8=800)，而3源和4源均接近于0。

圖1 不同源數(shù)量設(shè)定下PMF模型對編制數(shù)據(jù)的擬合度對比

不同源數(shù)量的模擬結(jié)果表明，當(dāng)設(shè)定源數(shù)量小于最優(yōu)值時，數(shù)據(jù)擬合度隨設(shè)定源數(shù)量增加而顯著提高，Q值則隨設(shè)定源數(shù)量增加而顯著下降。當(dāng)設(shè)定源數(shù)量等于最優(yōu)值時，繼續(xù)增加源數(shù)量模型擬合度和Q值呈現(xiàn)基本穩(wěn)定的趨勢，因此源數(shù)量最優(yōu)值可看作模型擬合度和Q值變化曲線的拐點(diǎn)。該結(jié)論與文獻(xiàn)報道的源數(shù)量選取方式基本一致[25]。

2.2 源解析矩陣與編制矩陣的對應(yīng)關(guān)系

圖1b顯示了模型對編制數(shù)據(jù)Xij較為理想的擬合程度，但仍有必要進(jìn)一步考察模型擬合所得的源貢獻(xiàn)率矩陣Gik和源成分譜矩陣Fkj與編制矩陣的對應(yīng)關(guān)系，因為在污染源解析領(lǐng)域，這2個分解矩陣更具有實際應(yīng)用價值。

對源解析文獻(xiàn)中常用的源成分譜百分比數(shù)據(jù)進(jìn)行考察，圖2中a為編制情景1實際所構(gòu)造的Fkj矩陣，b則為PMF模型擬合所得的Fkj矩陣，二者均轉(zhuǎn)化為百分比形式。表2則分別列出了二者未轉(zhuǎn)化為百分比的原始數(shù)據(jù)。

可以發(fā)現(xiàn)，圖2中的百分比堆積柱狀圖a與b一致程度較高，但表2中非百分比的原始數(shù)據(jù)a與b存在數(shù)量級級別的差異。進(jìn)一步查詢PMF擬合結(jié)果中的Gik矩陣，發(fā)現(xiàn)軟件對該矩陣進(jìn)行了歸一化處理(即Gik矩陣任一列向量的均值均為1)。但PMF模型分解出的Gik與Fkj矩陣相乘所得的Xij仍然與編制矩陣一致，見圖1b，說明在PMF擬合過程中，Gik任一列向量的均值被轉(zhuǎn)移至Fkj矩陣的對應(yīng)行向量中。該過程用矩陣形式表示如下(為便于區(qū)別，上標(biāo)帶*號的表示編制矩陣，不帶*號的表示PMF擬合矩陣，Ave_G表示Gik某一列向量的均值)：

(4)

將編制矩陣Fkj*的行向量乘上編制矩陣Gik*列向量均值(簡記為Ave_G×Fkj*)，原始數(shù)值列于表2，轉(zhuǎn)化為百分比形式繪制于圖2c，可以發(fā)現(xiàn)該原始數(shù)值及百分比與PMF擬合數(shù)據(jù)均有較高的一致性。而圖2a(即Fkj*)與圖2b(即Fkj)百分比占比較接近,可能是源于編制情景1中Gik*列向量均值較為接近，因其所使用的生成隨機(jī)數(shù)范圍均為1～100，參見表1情景1，均值都接近50，參見表2Gik*列向量均值。這一推測可從編制情景2的擬合結(jié)果中得到驗證。情景2改變了Gik*的構(gòu)造方式，使其列向量均值產(chǎn)生明顯大小差異，參見表1情景2。PMF擬合的Fkj原始數(shù)據(jù)與百分比形式分別列于表3和圖3?？梢钥闯?該情景下不管是原始數(shù)據(jù)還是百分比占比，F(xiàn)kj都與編制Fkj*差異較大，而與Ave_G×Fkj*相接近。

表2 編制情景1源成分譜對比

圖2 編制情景1源成分譜百分比對比度對比(F1～F3分別指代3個源Factor1～Factor3) 圖3 編制情景2源成分譜百分比對比度對比

實際上這種處理方式是有現(xiàn)實意義的。當(dāng)前利用PMF進(jìn)行源解析的研究中，均把源成分譜Fkj的百分比占比理解為源貢獻(xiàn)率[16,17,18-20,22-24]，如圖3b PMF源成分譜擬合結(jié)果中Cd元素的F1占比為58%，則意味著F1所代表的污染源對Cd的平均貢獻(xiàn)率為58%。在明確了PMF擬合結(jié)果中的源成分譜數(shù)據(jù)Fkj實際上為Ave_G×Fkj*后，可以證明源成分譜Fkj的百分比占比的物理含義就是源貢獻(xiàn)率?；诰幹茢?shù)據(jù)嘗試計算F1對樣品1中Cd含量的貢獻(xiàn)率，樣品1中Cd的總濃度為x11*(第1個下標(biāo)代表樣品編號i=1；第2個下標(biāo)代表重金屬元素編號j=1)，而根據(jù)公式(1)或(4)可知x11*的計算遵循公式(5)。而根據(jù)源貢獻(xiàn)率的定義，F(xiàn)1對樣品1中Cd含量的貢獻(xiàn)率P(F1)的計算遵循公式(6)。

(5)

(7)

推而廣之，F(xiàn)1對100個樣品中Cd含量的平均源貢獻(xiàn)率P(F1)ave的計算應(yīng)遵循公式(7)。從公式(7)可知，P(F1)ave實際上等于PMF源成分譜Fkj的百分比占比。從公式(7)還可以發(fā)現(xiàn)，F(xiàn)kj列向量的加和即為所有樣品中第j個元素的平均值，因此從表3的2行合計項中可以發(fā)現(xiàn)二者幾乎一致。

表3 編制情景2源成分譜對比

2.3 觀測誤差對源解析結(jié)果的影響

在編制情景1和2中，濃度數(shù)據(jù)矩陣Xij*直接等于Gik*與Fkj*的乘積，并未引入觀測誤差。而在實際應(yīng)用中誤差是不可避免的。因此在第3種編制情景中，利用Excel生成隨機(jī)數(shù)的方式為濃度數(shù)據(jù)矩陣Xij*引入隨機(jī)誤差，詳見表1編制情景3，并將該誤差絕對值的2倍作為輸入PMF模型的不確定度數(shù)據(jù)矩陣ij*，考察帶有誤差條件下PMF的擬合情況。

情景3與情景2的樣品濃度與不確定度的對比如圖4所示。在情景2中，由于不確定度數(shù)據(jù)μij*均按照公式(3)進(jìn)行設(shè)計，因此每個樣品的不確定度均與濃度成比例，且相應(yīng)的信噪比S/N均高達(dá)10.0。而情景3引入隨機(jī)誤差后，信噪比明顯下降，更符合實際應(yīng)用情況。但信噪比S/N仍滿足大于1.0的條件，表明樣品濃度數(shù)據(jù)的信號強(qiáng)度滿足模型擬合要求。

圖4 樣品濃度與不確定度對比

將情景3的模型擬合結(jié)果與該情景編制濃度數(shù)據(jù)進(jìn)行對比，結(jié)果如圖5所示。在圖1b中，對于未引入誤差的情景1，PMF模型取得了很好的擬合效果，3源情況下各元素的r2值均達(dá)到1.0。引入隨機(jī)誤差后，PMF模型對情景3的擬合優(yōu)度有所下降，但仍達(dá)到了相對較高的水平，各元素的r2值在0.87～0.99，對于實際應(yīng)用而言是可以接受的。因誤差的引入，PMF擬合Qture值相比情景2也有一定增長，為118.02，但仍遠(yuǎn)小于理論Q值800，表明模型擬合結(jié)果較好。

圖5 PMF模型對情景3編制數(shù)據(jù)的擬合情況 圖6 編制情景3源成分譜百分比對比度對比

進(jìn)一步考察在土壤重金屬PMF源解析領(lǐng)域中最為關(guān)注的源貢獻(xiàn)率的情況，源貢獻(xiàn)率的百分比占比和原始數(shù)值分別參見圖6和表4。從圖6可以看出，盡管因為觀測誤差的存在導(dǎo)致PMF模型對原始數(shù)據(jù)的擬合度有一定下降，但從源成分譜的百分比占比來看，PMF擬合結(jié)果與原始編制數(shù)據(jù)所表征的源貢獻(xiàn)率是基本一致的。擬合所得的各元素的主導(dǎo)貢獻(xiàn)源均未發(fā)生偏離。而從表4可以看出，與表3未引入誤差的情景2類似，PMF擬合Fkj的合計項與Ave_G×Fkj*的合計項仍保持高度一致。這表明在原始濃度數(shù)據(jù)帶有誤差的條件下，PMF模型仍能取得較為理想的擬合結(jié)果以及源貢獻(xiàn)率信息。

表4 編制情景2源成分譜對比

3 結(jié)論

本文利用自行編制設(shè)計的數(shù)據(jù)集作為PMF模型的原始輸入，在已知數(shù)據(jù)集的源成分譜及貢獻(xiàn)率的條件下驗證了PMF模型的源解析結(jié)果，主要結(jié)論如下。

在PMF模型擬合過程中需要設(shè)定源數(shù)量，源數(shù)量最優(yōu)值可近似按照模型擬合度和Q值變化曲線的拐點(diǎn)確定。當(dāng)設(shè)定源數(shù)量小于最優(yōu)值時，數(shù)據(jù)擬合度隨設(shè)定源數(shù)量增加而顯著提高，Q值則隨設(shè)定源數(shù)量增加而顯著下降。當(dāng)設(shè)定源數(shù)量等于最優(yōu)值時，繼續(xù)增加源數(shù)量模型擬合度和Q值呈現(xiàn)基本穩(wěn)定的趨勢。

在不帶觀測誤差和帶有觀測誤差2種條件下，PMF模型均能較好地擬合原始數(shù)據(jù)，其中不帶誤差條件下擬合r2值可達(dá)1.0，帶有誤差條件下擬合度有所下降，但仍高于0.87。

PMF擬合所得的源成分譜Fkj實際為編制源貢獻(xiàn)矩陣列向量均值與編制源成分譜矩陣行向量的乘積，其物理含義即為某個源對某種重金屬的平均貢獻(xiàn)含量，因此源成分譜的百分比占比的物理含義即為源貢獻(xiàn)率。在不帶觀測誤差和帶有觀測誤差2種條件下，PMF擬合所得源貢獻(xiàn)率與原始編制數(shù)據(jù)所表征的源貢獻(xiàn)率一致性較好，表明在滿足模型假設(shè)的條件下，PMF用于土壤重金屬源解析可取得較為理想的結(jié)果。