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

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

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

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

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

本文嘗試利用自行編制設計的數(shù)據(jù)集作為PMF模型的原始輸入，在已知數(shù)據(jù)集的源成分譜及貢獻率的條件下驗證PMF模型源解析結果的準確性，以便更好地理解PMF模型的工作原理以及擬合結果的現(xiàn)實意義，為土壤重金屬污染源解析工作下階段的進一步鋪開完善方法論。

1 研究方法

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

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

Xij=GikFkj+Eij

(1)

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

(2)

式中，en,m為殘差矩陣Eij中第n行第m列中的殘差數(shù)值；μn,m為與濃度數(shù)據(jù)矩陣Xij相對應的不確定度矩陣μ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，遠小于任一濃度數(shù)據(jù)xn,m值。

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

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

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

2 結果與討論

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

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

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

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

2.2 源解析矩陣與編制矩陣的對應關系

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

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

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

(4)

將編制矩陣Fkj*的行向量乘上編制矩陣Gik*列向量均值(簡記為Ave_G×Fkj*)，原始數(shù)值列于表2，轉化為百分比形式繪制于圖2c，可以發(fā)現(xiàn)該原始數(shù)值及百分比與PMF擬合數(shù)據(jù)均有較高的一致性。而圖2a(即Fkj*)與圖2b(即Fkj)百分比占比較接近,可能是源于編制情景1中Gik*列向量均值較為接近，因其所使用的生成隨機數(shù)范圍均為1～100，參見表1情景1，均值都接近50，參見表2Gik*列向量均值。這一推測可從編制情景2的擬合結果中得到驗證。情景2改變了Gik*的構造方式，使其列向量均值產(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)實意義的。當前利用PMF進行源解析的研究中，均把源成分譜Fkj的百分比占比理解為源貢獻率[16,17,18-20,22-24]，如圖3b PMF源成分譜擬合結果中Cd元素的F1占比為58%，則意味著F1所代表的污染源對Cd的平均貢獻率為58%。在明確了PMF擬合結果中的源成分譜數(shù)據(jù)Fkj實際上為Ave_G×Fkj*后，可以證明源成分譜Fkj的百分比占比的物理含義就是源貢獻率。基于編制數(shù)據(jù)嘗試計算F1對樣品1中Cd含量的貢獻率，樣品1中Cd的總濃度為x11*(第1個下標代表樣品編號i=1；第2個下標代表重金屬元素編號j=1)，而根據(jù)公式(1)或(4)可知x11*的計算遵循公式(5)。而根據(jù)源貢獻率的定義，F(xiàn)1對樣品1中Cd含量的貢獻率P(F1)的計算遵循公式(6)。

(5)

(7)

推而廣之，F(xiàn)1對100個樣品中Cd含量的平均源貢獻率P(F1)ave的計算應遵循公式(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 觀測誤差對源解析結果的影響

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

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

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

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

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

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

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

3 結論

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

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

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

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