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

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

土壤是經(jīng)濟(jì)社會(huì)可持續(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ù)計(jì)該項(xiàng)工作將成為我國土壤污染防治“十四五”規(guī)劃的重要工作進(jìn)一步鋪開。在已發(fā)布的浙江省《土壤、地下水和農(nóng)業(yè)農(nóng)村污染防治“十四五”規(guī)劃》中，土壤污染溯源被列為重點(diǎn)任務(wù)[3]。土壤重金屬的來源往往包括自然源和多種人為源等，溯源工作較為復(fù)雜[4]。現(xiàn)有的土壤污染源解析方法大致包括多元統(tǒng)計(jì)分析法[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)算。迭代計(jì)算結(jié)果可用于定量確定各污染源的成分譜(FactorProfile)及貢獻(xiàn)率(FactorContributions)[12]。目前土壤污染源解析領(lǐng)域利用PMF模型的主要目的在于確定污染源的貢獻(xiàn)占比[18]。由于無法確切得知土壤在成土過程及人為輸入影響中重金屬的真實(shí)來源及貢獻(xiàn)程度，因此試圖對模型擬合所得的貢獻(xiàn)占比的真假或優(yōu)劣作出評價(jià)是相當(dāng)困難的。目前對模型結(jié)果可信度的判斷主要借助軟件自帶的不確定度估計(jì)功能[19-21]或與其它受體模型結(jié)果進(jìn)行交叉對比[10,16,22-24]等手段。但仍然缺乏模型結(jié)果與真實(shí)結(jié)果的對比驗(yàn)證。

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

1 研究方法

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

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)識，不代表實(shí)際意義。而PMF源解析的目標(biāo)在于將樣品濃度數(shù)據(jù)矩陣Xij分解成源貢獻(xiàn)率矩陣Gik和源成分譜矩陣Fkj，以及一個(gè)殘差矩陣Eij，如公式(1)所示[24]：

Xij=GikFkj+Eij

(1)

式中，Xij為i個(gè)樣品j個(gè)重金屬數(shù)據(jù)組成的i行j列矩陣；Gik為i個(gè)樣品k個(gè)源(或稱因子)的貢獻(xiàn)數(shù)據(jù)組成的i行k列矩陣，代表某個(gè)源在某個(gè)樣品中的貢獻(xiàn)大?。籉kj為k個(gè)源j個(gè)重金屬的濃度數(shù)據(jù)組成的k行j列矩陣，代表某個(gè)源中某種重金屬的含量高低；Eij為i行j列的殘差矩陣，代表某個(gè)樣品中未被模型解釋的某種重金屬含量的大小。PMF模型通過多次迭代計(jì)算獲得最佳分解矩陣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由于采樣和分析測試等誤差所帶來的不確定度，可由測試實(shí)驗(yàn)室提供或用戶根據(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è)計(jì)3種編制情景，每種情景的特征匯總見表1。所有情景樣品數(shù)i均為100，元素?cái)?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í)，數(shù)據(jù)擬合度隨設(shè)定源數(shù)量增加而顯著提高，Q值則隨設(shè)定源數(shù)量增加而顯著下降。當(dāng)設(shè)定源數(shù)量等于最優(yōu)值時(shí)，繼續(xù)增加源數(shù)量模型擬合度和Q值呈現(xiàn)基本穩(wěn)定的趨勢，因此源數(shù)量最優(yōu)值可看作模型擬合度和Q值變化曲線的拐點(diǎn)。該結(jié)論與文獻(xiàn)報(bào)道的源數(shù)量選取方式基本一致[25]。

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

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

對源解析文獻(xiàn)中常用的源成分譜百分比數(shù)據(jù)進(jìn)行考察，圖2中a為編制情景1實(shí)際所構(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é)果中得到驗(yàn)證。情景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個(gè)源Factor1～Factor3) 圖3 編制情景2源成分譜百分比對比度對比

實(shí)際上這種處理方式是有現(xiàn)實(shí)意義的。當(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實(shí)際上為Ave_G×Fkj*后，可以證明源成分譜Fkj的百分比占比的物理含義就是源貢獻(xiàn)率?；诰幹茢?shù)據(jù)嘗試計(jì)算F1對樣品1中Cd含量的貢獻(xiàn)率，樣品1中Cd的總濃度為x11*(第1個(gè)下標(biāo)代表樣品編號i=1；第2個(gè)下標(biāo)代表重金屬元素編號j=1)，而根據(jù)公式(1)或(4)可知x11*的計(jì)算遵循公式(5)。而根據(jù)源貢獻(xiàn)率的定義，F(xiàn)1對樣品1中Cd含量的貢獻(xiàn)率P(F1)的計(jì)算遵循公式(6)。

(5)

(7)

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

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

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

在編制情景1和2中，濃度數(shù)據(jù)矩陣Xij*直接等于Gik*與Fkj*的乘積，并未引入觀測誤差。而在實(shí)際應(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è)計(jì)，因此每個(gè)樣品的不確定度均與濃度成比例，且相應(yīng)的信噪比S/N均高達(dá)10.0。而情景3引入隨機(jī)誤差后，信噪比明顯下降，更符合實(shí)際應(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，對于實(shí)際應(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可以看出，盡管因?yàn)橛^測誤差的存在導(dǎo)致PMF模型對原始數(shù)據(jù)的擬合度有一定下降，但從源成分譜的百分比占比來看，PMF擬合結(jié)果與原始編制數(shù)據(jù)所表征的源貢獻(xiàn)率是基本一致的。擬合所得的各元素的主導(dǎo)貢獻(xiàn)源均未發(fā)生偏離。而從表4可以看出，與表3未引入誤差的情景2類似，PMF擬合Fkj的合計(jì)項(xiàng)與Ave_G×Fkj*的合計(jì)項(xiàng)仍保持高度一致。這表明在原始濃度數(shù)據(jù)帶有誤差的條件下，PMF模型仍能取得較為理想的擬合結(jié)果以及源貢獻(xiàn)率信息。

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

3 結(jié)論

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

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

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

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