基于PMF模型的湖北隨縣三里崗地區(qū)土壤重金屬源解析

董健彪，謝淑云，田 歡，楊文兵，李 華，周徐維，Diego Armando Pinzon Nunez

(1.中國地質(zhì)大學(xué)(武漢)地球科學(xué)學(xué)院，湖北 武漢 430074； 2.中國地質(zhì)大學(xué)(武漢)材料與化學(xué)學(xué)院，湖北 武漢 430074； 3.湖北省地質(zhì)局第八地質(zhì)大隊(duì)，湖北 襄陽 441002； 4.青海省地質(zhì)調(diào)查局，青海 西寧 810000)

1 引言

近年來，大量重金屬元素通過采礦、工業(yè)廢料排放、農(nóng)藥、化肥等途徑進(jìn)入土壤系統(tǒng)，同時由于其在環(huán)境中的持久性、生物累積性、難降解性及其與有益元素結(jié)合造成新的風(fēng)險(xiǎn)等特點(diǎn)，受到人們的廣泛關(guān)注(楊希 等，2021；Zhang et al.，2021)。2014年《全國土壤污染狀況調(diào)查公報(bào)》顯示，我國土壤重金屬污染情況研究也受到廣泛關(guān)注，點(diǎn)位超標(biāo)率達(dá)到了16.1%，其中As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八種重金屬是最主要的無機(jī)污染物。

土壤重金屬源解析工作是污染評價(jià)及治理工作的前提，具有重要的現(xiàn)實(shí)意義(Hu et al.，2018)。學(xué)界流行的源解析方法很多，如主成分分析法(PCA)、正定矩陣因子分解法(Positive Matrix Factorization)等(董騄睿 等，2015)。PMF模型是美國環(huán)保署推薦的源解析方法，起初多用于研究大氣以及水污染，該方法在識別污染源貢獻(xiàn)比例上表現(xiàn)良好，且在與GIS技術(shù)結(jié)合時可以更有效地識別污染源(Dong et al.，2019；Guan et al.，2018；Liang et al.，2017)，因此近年來越來越多的學(xué)者利用PMF模型識別土壤重金屬污染來源(Hu et al.，2018)。

三里崗鎮(zhèn)地處湖北省隨州市，是湖北省香菇主產(chǎn)區(qū)之一，香菇產(chǎn)品遠(yuǎn)銷海外，近年出口香菇中出現(xiàn)重金屬超標(biāo)產(chǎn)品，對于研究區(qū)農(nóng)產(chǎn)品銷量造成一定影響(楊定國 等，2014)，因此有必要對研究區(qū)展開土壤地球化學(xué)調(diào)查，探討土壤中重金屬元素的含量。本文運(yùn)用PMF模型結(jié)合GIS方法，以三里崗鎮(zhèn)為研究區(qū)，對土壤樣品中的As、Cd、Cr、Cu、Ni、Pb、Zn和Hg元素的含量和來源進(jìn)行評價(jià)和解析，并結(jié)合GIS方法分析不同污染源的源貢獻(xiàn)率分布情況，以期為土壤重金屬污染的治理和管控提供理論依據(jù)和技術(shù)支持。

2 材料與方法

2.1 區(qū)域地質(zhì)概況

圖1 三里崗鎮(zhèn)地質(zhì)圖(改自陳婕 等(2018))Fig.1 Geological map of Sanligang town(modified after Chen Jie et al.2018)1—第四紀(jì)；2—白堊紀(jì)；3—泥盆紀(jì)；4—志留紀(jì)；5—奧陶紀(jì)-志留紀(jì)；6—寒武紀(jì)；7— 震旦紀(jì)；8—元古代；9—一級構(gòu)造單元界限；10—二級構(gòu)造單元界限；11—三級構(gòu)造單元界限

2.2 數(shù)據(jù)來源與處理

本文數(shù)據(jù)來源于湖北省“金土地”工程-高標(biāo)準(zhǔn)基本農(nóng)田地球化學(xué)調(diào)查(項(xiàng)目編號HBJTD20140108)，該項(xiàng)目由湖北省地質(zhì)局第八地質(zhì)大隊(duì)主持完成，完成了對隨縣三里崗鎮(zhèn)的1∶50000表層土壤樣測量，平均采樣密度為 5.19 個點(diǎn)/km2。共采集表層土壤1696件，其中重復(fù)樣34件。

本文使用Excel 2019和SPSS 21進(jìn)行數(shù)據(jù)整理、計(jì)算和分析，利用GS+(Version7)進(jìn)行半變異函數(shù)分析，通過EPA PMF 5.0進(jìn)行重金屬源解析，利用ArcGIS 10.2進(jìn)行空間分析和圖件制作。

2.3 污染評價(jià)方法

本文利用單因子指數(shù)法和內(nèi)梅羅指數(shù)法對研究區(qū)重金屬污染情況進(jìn)行評價(jià)(楊濤毅 等，2011；Meng et al.，2011)。單因子指數(shù)法是將土壤中某污染元素的含量值Ci除以評價(jià)標(biāo)準(zhǔn)值Si，根據(jù)結(jié)果Pi對土壤中單一污染元素進(jìn)行污染評價(jià)。公式如下：

Pi=Ci/Si

(1)

內(nèi)梅羅綜合污染指數(shù)法利用多種元素的單因子污染指數(shù)的平均值(Pivag)以及最大值(Pimax)進(jìn)行計(jì)算，兼顧了多元素的綜合特征以及單一因素的極端影響。公式如下：

(2)

2.4 半變異函數(shù)

半變異函數(shù)是地質(zhì)統(tǒng)計(jì)學(xué)中用來描述區(qū)域化變量性質(zhì)的一種常用方法，通過半變異函數(shù)，可以對區(qū)域化變量的隨機(jī)性特征以及結(jié)構(gòu)性特征進(jìn)行分析，得到研究變量的空間分布規(guī)律(耿治鵬 等，2022)。半變異函數(shù)公式為

(3)

式中，y(h)為半變異函數(shù)值，N(h)表示距離向量為h時的采樣點(diǎn)數(shù)目，h表示空間距離，Z(x)為區(qū)域化變量，一般表示為

Z(x)=(xw，xv，xw)

(4)

2.5 PMF模型

PMF全稱正定因子矩陣分解法(Positive Matrix Factorization)，是由Paatero et al.(1994)基于因子分析原理所提出的一種分析模型，利用樣品的濃度和不確定度，通過最小二乘法計(jì)算因子數(shù)(徐源，2021)。

PMF模型將重金屬元素含量視作i×j階的矩陣，將該矩陣拆分為三個因子矩陣，分別是因子貢獻(xiàn)矩陣G(i×k)，因子成分矩陣F(k×j)以及殘差矩陣E(i×j)，公式表示為

(5)

其中Xij為原始矩陣第i個樣品第j個化學(xué)元素，Gik為第k個因子對于第i個樣品的因子貢獻(xiàn)度，F(xiàn)jk是第k個因子對第j個元素的貢獻(xiàn)濃度，p則表示因子數(shù)(薛建龍，2014)。

PMF模型定義了一個目標(biāo)函數(shù)Q：

(6)

式中，c表示樣品序數(shù)，d表示元素序數(shù)，k表示因子的序數(shù)，因此zcd表示第c個樣品第d個元素的含量值，gck表示第k個因子對第c個樣品的貢獻(xiàn)度，dk表示第k個因子對第d個元素的貢獻(xiàn)濃度，ucd表示第c個樣品第d個元素的不確定度，ecd為殘差。

不確定度計(jì)算公式為

(7)

當(dāng)元素濃度小于或等于檢出限(MDL)時，的不確定度計(jì)算公式為

(8)

3 結(jié)果與討論

3.1 土壤重金屬元素參數(shù)特征

表1 隨縣土壤重金屬元素描述性統(tǒng)計(jì)表①Table 1 Descriptive statistical table of heavy metal elements in soil of Sui county

依據(jù)GB15618-2018《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)(試行)》(以下簡稱《標(biāo)準(zhǔn)》)中的風(fēng)險(xiǎn)篩選值為評價(jià)標(biāo)準(zhǔn)值對研究區(qū)水田和其他用地樣品分別統(tǒng)計(jì)分析，結(jié)果見表2。表2顯示，兩類用地中，Cd元素污染樣品最多，水田中超標(biāo)樣品數(shù)為32，占比為5.6%，其他用地中超標(biāo)樣品數(shù)為379，占比為38.9%，其次是As元素，在兩類用地中超標(biāo)樣品數(shù)均超過了5%，而Cr、Ni、Zn元素超標(biāo)樣品數(shù)均在2%左右， Cu、Pb、Hg元素?zé)o超標(biāo)樣品。整體來看，研究區(qū)Cd元素超標(biāo)率遠(yuǎn)遠(yuǎn)超過《全國污染調(diào)查公報(bào)》中的全國點(diǎn)位平均超標(biāo)率7%，存在一定的風(fēng)險(xiǎn)(呂悅風(fēng) 等，2019)。

基于單因子指數(shù)法與內(nèi)梅羅指數(shù)法對研究區(qū)水田與其他用地土壤中重金屬進(jìn)行評價(jià)，評價(jià)標(biāo)準(zhǔn)值表2已給出，評價(jià)結(jié)果見表3，污染評價(jià)標(biāo)準(zhǔn)見表4(鄧通德 等，2022)。從表3中可以看出研究區(qū)8種重金屬元素的單因子指數(shù)均值皆小于1，研究區(qū)8種重金屬在全區(qū)的平均污染程度不高，除Cu、Pb、Hg外其余元素均有污染樣品，但只有Cd元素污染樣品數(shù)超過了20%，形成一定規(guī)模。綜合污染指數(shù)P綜均值為0.73，根據(jù)表4中評價(jià)標(biāo)準(zhǔn)，研究區(qū)整體污染程度處于尚清潔等級，有253件P綜>1的污染樣品，占比達(dá)到16.39%，說明研究區(qū)污染情況值得重視。

表2 重金屬超標(biāo)樣品數(shù)目統(tǒng)計(jì)②Table 2 Statistics of the number of heavy metal samples exceeding the standard

表3 重金屬污染指數(shù)描述統(tǒng)計(jì)量③Table 3 Descriptive statistics of heavy metal pollution index

續(xù)表

表4 污染評價(jià)標(biāo)準(zhǔn)Table 4 Standards for pollution assessment

3.2 土壤半變異函數(shù)分析

利用GS+軟件對數(shù)據(jù)進(jìn)行半變異函數(shù)分析，并根據(jù)決定系數(shù)R2越趨于1、殘差平方和RSS越趨于0則模型擬合度越好的原則選取合適的模型，半變異函數(shù)分析要求數(shù)據(jù)符合正態(tài)分布或近正態(tài)分布(吳敏 等，2016；蔡大為 等，2020)，而8種元素原數(shù)據(jù)都不符合正態(tài)分布特征，進(jìn)行半變異函數(shù)分析前對數(shù)據(jù)進(jìn)行轉(zhuǎn)換，轉(zhuǎn)換前后的偏峰系數(shù)對比見表5。半變異函數(shù)分析結(jié)果見表6，8種元素的決定系數(shù)R2均大于0.7，殘差平方和都比較小，說明模型擬合較好。

表5 偏峰系數(shù)檢驗(yàn)Table 5 Test of skewness coefficient and kurtosis coefficient

表6 半變異函數(shù)分析Table 6 Semi-variogram analysis

3.3 土壤重金屬PMF源解析

利用EPA PMF 5.0對8種重金屬元素進(jìn)行PMF源解析，根據(jù)Qrob/Qexp比值快速下降處因子數(shù)最佳(李嬌 等，2019)將因子數(shù)確定為4。表7列出了4個因子對8種元素的源成分譜和源貢獻(xiàn)率。圖2是各因子對各元素的源貢獻(xiàn)率，利用ArcGIS對PMF源貢獻(xiàn)率進(jìn)行克里格插值，繪制出4個因子的源貢獻(xiàn)率分布圖(圖3)。

表7 PMF源解析出的各源成分譜及其貢獻(xiàn)率Table 7 Component spectrum and contribution rate of each source resolved by PMF source

圖2 各元素因子源貢獻(xiàn)率Fig.2 Factor source contribution rate of each element

圖3 因子源貢獻(xiàn)率分布圖Fig.3 Distribution of factor source contribution rates

根據(jù)表7及圖2，因子1主要代表As元素的來源，對As元素的源貢獻(xiàn)率高達(dá)71.78%，對Pb元素也有一定的源貢獻(xiàn)率(22.45%)。由于我國煤資源中As元素含量很高(白向飛，2003)，高砷煤燃燒過程中As便會進(jìn)入環(huán)境中(Finkelman et al.，2018)，但煤燃燒同樣會導(dǎo)致Hg污染，而因子1對于Hg的源貢獻(xiàn)率卻很低，因此因子1應(yīng)代表其他來源。通過觀察圖1和圖3中因子1源貢獻(xiàn)率分布圖可以看出，因子1高源貢獻(xiàn)率與震旦系(Z)以及白堊系(K)地層重合度很高，因此推測因子1可能與地層有較大關(guān)聯(lián)。表8統(tǒng)計(jì)了1544組樣品、不同年代地層中8種重金屬元素的含量，有35組樣品屬于巖漿巖，由于占比少，未參與統(tǒng)計(jì)。從表6中可以看出，震旦系以及白堊系地層中的As含量(15.13 mg/kg和14.79 mg/kg)和Pb(25.70 mg/kg和26.66 mg/kg)含量相比其他地層更高，同時也高于研究區(qū)含量均值(As：11.52 mg/kg；Pb：24.42 mg/kg)。因此因子1代表的來源主要是震旦系和白堊系地層母巖。

表8 不同年代地層8種重金屬元素含量④Table 8 Contents of eight heavy metals in strata of different ages

從表7及圖2中看出，因子2對Hg、Pb兩種元素的源貢獻(xiàn)率高達(dá)76.96%、63.28%，對As、Cd、Zn等元素也表現(xiàn)出了一定的影響，源貢獻(xiàn)率分別為19.00%、16.24%、25.80%。有研究表明，交通排放是一個重要的重金屬(包括Hg、Pb、Cd等元素)來源(Liu et al.，2012；Yang et al.，2017)，而三里崗鄉(xiāng)村道路分布全區(qū)，從圖3也可以看出，因子2源貢獻(xiàn)率超過24.20%的區(qū)域在全區(qū)分布廣泛，因此因子2應(yīng)該受到了交通排放所帶來的影響，但因?yàn)檠芯繀^(qū)為鄉(xiāng)鎮(zhèn)，車輛流通量小，所以交通排放到土壤中的重金屬含量并不高，所以并未有樣品的Hg、Pb含量超過風(fēng)險(xiǎn)篩選值。其次，圖3中因子2高源貢獻(xiàn)率區(qū)域(>34.34%)分布在北部、中部以及南部，與圖1元古界(Pt)地層重合程度較高，說明因子2可能也受到了來自元古界地層母巖的影響。半變異函數(shù)分析表明Pb受到自然因素和人為因素共同影響，因此因子2解釋為元古界地層母巖與道路交通排放的共同影響較為合理。

從表7以及圖2中看出，因子4對于Cr、Cu、Ni的影響最高，分別達(dá)到了88.40%、82.14%和88.77%，同時對Zn元素的貢獻(xiàn)率也達(dá)到了42.46%。從圖1及圖3中因子4源貢獻(xiàn)率分布可以看出，因子4高源貢獻(xiàn)率(>39.87%)區(qū)域與奧陶系到志留系(O-S)地層高度重合，從表8中也可以看出，Cr、Cu、Ni在寒武系到志留系地層中的含量遠(yuǎn)高于其余地層中的含量，Zn元素含量也高于研究區(qū)背景值，所以因子4可以解釋為奧陶系到志留系地層母巖的影響。

4 結(jié)論

本文利用半變異函數(shù)及PMF模型對三里崗鎮(zhèn)8種重金屬元素的含量及來源進(jìn)行了分析，結(jié)果表明：

(1)8種重金屬元素除Hg、Pb外的均值都高于全國背景值，依據(jù)《標(biāo)準(zhǔn)》中所給出的風(fēng)險(xiǎn)篩選值進(jìn)行評價(jià)，Cd元素超標(biāo)率超過20%，內(nèi)梅羅綜合污染指數(shù)法顯示三里崗鎮(zhèn)污染樣品占比為16.39%，平均污染指數(shù)為0.73，處于尚清潔水平。

(2)半變異函數(shù)分析結(jié)果表明，8種重金屬元素的塊金效應(yīng)均<50%，表明研究區(qū)重金屬元素含量主要受自然因素影響。

(3)PMF源解析得到了4個來源，依據(jù)不同地層中重金屬元素的含量以及因子源貢獻(xiàn)率分布情況將因子1解釋為震旦系和白堊系地層母巖的影響、因子2解釋為元古界地層母巖與鄉(xiāng)村道路交通排放的混合源、因子3解釋為震旦系和寒武系地層母巖的影響、因子4解釋為奧陶系到志留系地層母巖的影響。

致謝：本研究數(shù)據(jù)主要來自湖北省地質(zhì)局第八地質(zhì)大隊(duì)和“金土地”工程-高標(biāo)準(zhǔn)基本農(nóng)田地球化學(xué)調(diào)查課題(項(xiàng)目編號HBJTD20140108)，受中國地質(zhì)大學(xué)(武漢)中央高校基金生態(tài)地球化學(xué)團(tuán)隊(duì)項(xiàng)目(CUG170104)聯(lián)合資助。