湖南省郴州市蘇仙區(qū)重點污染企業(yè)影響區(qū)的土壤重金屬污染源解析

徐 源， 師華定,*， 王 超， 費 楊， 于靖靖， 舒 密

1.中國環(huán)境科學研究院土壤與固體廢物環(huán)境研究所， 北京 100012

2.生態(tài)環(huán)境部土壤與農(nóng)業(yè)農(nóng)村生態(tài)環(huán)境監(jiān)管技術(shù)中心， 北京 100012

3.北京信息科技大學， 北京 100101

近年來，隨著城市化進程的加快，重金屬通過污水灌溉、大氣干濕沉降、工業(yè)廢渣等途徑進入土壤[1-2]，重金屬污染已成為當今土壤污染中污染面積最廣、危害最大的環(huán)境問題之一[3-4]. Cd、Hg、As、Pb、Cr作為重金屬“五毒”元素，在土壤中具有移動性大、毒性高、無法降解等特點，在生產(chǎn)活動中容易被作物吸收富集，不僅嚴重影響作物的產(chǎn)量和品質(zhì)，還可以通過食物鏈在人體積累，危害人體健康[5-7]. 研究土壤“五毒”元素污染的來源以及污染源的主要影響區(qū)域，對維持區(qū)域環(huán)境健康具有重要意義[8-9].

目前，土壤重金屬污染來源解析主要利用基于重金屬總量的受體模型，常用的受體模型有化學質(zhì)量平衡(CMB)模型、主成分分析/因子分析-多元線性回歸(PCA/FA-MLR)模型、正定矩陣因子分解(PMF)模型、UNMIX模型等[10-11]. 其中，PMF模型具有適用性廣、不需要測量源成分譜、能對因子貢獻作非負約束等優(yōu)點，得到廣泛的應用. 例如，黃華斌等[12]利用PMF模型解析出研究區(qū)土壤重金屬的來源；董騄睿等[13]通過PMF模型解析結(jié)果與Pb穩(wěn)定同位素比值結(jié)果進行對比分析，指出PMF模型可以較好地應用于土壤重金屬源解析研究.

然而，單一的PMF模型的解析結(jié)果往往比較籠統(tǒng)，對污染源的判別中主觀性較強，且缺乏直觀視覺效果，因此較多學者將其與地統(tǒng)計學方法結(jié)合運用，進一步得到污染源在空間上的分布狀況[14-16]. 這類研究雖對解析結(jié)果有了很大改進，但在源成分譜未知的情況下，對于污染物的具體來源依然無法準確識別. 另外，關(guān)于PMF模型對土壤重金屬的源解析研究中，目前大多僅停留在污染源大類上，如自然源、工業(yè)源、農(nóng)業(yè)源等，其中對工業(yè)源的探討并未涉及具體的重點污染行業(yè). 雙變量莫蘭指數(shù)方法(bivariate Moran′s I)通過指數(shù)的計算對空間內(nèi)兩種要素作比較，表征兩種要素間的空間相關(guān)性. 周俊馳[17]采用莫蘭指數(shù)對耕地土壤重金屬進行分析，識別出研究區(qū)內(nèi)Cd、As、Hg這3種重金屬在耕地土壤中的富集與工礦開采排放關(guān)系密切；于靖靖等[18]利用莫蘭指數(shù)和廣義加性模型探討了重點污染企業(yè)對重金屬Cd的影響，表明依據(jù)空間自相關(guān)理論可以對特定的污染源做很好的識別. 鑒于此，該研究將地統(tǒng)計學方法與PMF模型結(jié)合起來，在利用傳統(tǒng)PMF模型進行污染源解析的基礎(chǔ)上，利用莫蘭指數(shù)表征不同行業(yè)與污染源的空間關(guān)系，明確工業(yè)源下的不同重點污染企業(yè)對土壤重金屬的影響，以輔助解析驗證PMF模型的有效性，以期為區(qū)域企業(yè)管理和土壤污染治理提供參考.

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)域位于湖南省郴州市蘇仙區(qū)(112°53′55″E～113°16′22″E、25°30′21″N～26°03′29″N)，海拔150～160 m，以山地為主，崗地、水面較少，年降雨量 1 452.1 mm，屬中亞熱帶季風濕潤氣候，土地總面積為 1 342.27 km2. 該區(qū)域地處郴州市西河流域，西河是郴江的一級支流、湘江的三級支流. 蘇仙區(qū)礦產(chǎn)資源豐富，素有“有色金屬之鄉(xiāng)”的美譽，有世界罕見的柿竹園多金屬礦，國有硚口鉛鋅礦，瑪瑙山錳礦，東波鉛鋅礦以及年產(chǎn)煤數(shù)萬噸的許家洞、街洞和棲鳳渡煤礦等. 發(fā)達的礦業(yè)生產(chǎn)給當?shù)貛砹艘幌盗协h(huán)境問題.

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

利用ArcGIS 10.3在研究區(qū)范圍內(nèi)按照2 km×2 km布設(shè)302個網(wǎng)格，用GPS對302個采樣點精確定位. 采樣時采用梅花形布點法對土壤樣品進行采集，用木鏟采集每個網(wǎng)格表層(0～20 cm)的土壤樣品至少6個，樣品混合均勻后按四分法獲取每個采樣點采集的土壤不少于2 kg，裝入樣品袋中帶回，并記錄現(xiàn)場土樣的相關(guān)信息及周邊土地利用現(xiàn)狀. 將采集的土壤樣品棄去雜草、礫石、動植物殘體等雜物，混勻風干后用研缽磨碎過100目(孔徑相當于0.15 mm)篩后保存. 采用HCl-HNO3-HF微波密閉消解技術(shù)進行土壤樣品消解，采用電感耦合等離子體原子發(fā)射光譜法(ICP-AES)測定土壤中w(Cd)、w(Cr)、w(Pb)，采用原子熒光法測定土壤中w(As)、w(Hg). 研究區(qū)采樣點及主要潛在污染源分布見圖1.

圖1 蘇仙區(qū)采樣點及主要潛在污染源分布概況

1.3 數(shù)據(jù)處理

研究區(qū)土壤中5種重金屬元素的描述性統(tǒng)計分析、正態(tài)分布檢驗和轉(zhuǎn)換、Pearson相關(guān)性分析，借助SPSS 19.0軟件完成；樣品污染源解析借助PMF 5.0完成；空間自相關(guān)分析借助GeoDa軟件完成；克里金插值借助ArcCIS 10.3完成.

1.4 研究方法

1.4.1正定矩陣因子分解(PMF)模型

PMF模型是由Paatero等[19]提出的一種因子分析受體模型. 該模型將采樣數(shù)據(jù)矩陣(X)分解成因子貢獻矩陣(G)、因子成分矩陣(F)以及殘差矩陣(E)[20-21]：

(1)

式中，a為受體樣品個數(shù)，b為所測的化學物質(zhì)種類，p為主因子數(shù)(即主要源個數(shù)).

PMF模型基于加權(quán)最小二乘法進行限定和迭代計算，利用樣品的重金屬濃度和不確定度數(shù)據(jù)進行各樣點的加權(quán)計算，使得目標函數(shù)Q最小化. 目標函數(shù)Q定義[20]如下：

式中，zcd為第c(c=1,2,…,a)個樣品中第d(d=1,2,…,b)個元素的含量，gck為源k對第c個樣品的相對貢獻，ucd為第c個樣品中第d個元素含量的不確定性大小，fkd為源k中第d個元素的含量，ecd為殘差.

1.4.2莫蘭指數(shù)

該研究利用莫蘭指數(shù)作為測度指標，探討屬性值之間是否具有特殊的空間形態(tài)[22-23]，分為單變量莫蘭指數(shù)和雙變量莫蘭指數(shù). 其中，單變量莫蘭指數(shù)可以指出區(qū)域同一屬性值的分布是聚集、離散或者隨機模式[24-25]，雙變量莫蘭指數(shù)揭示了空間中某一要素的一個指標與其相鄰位置要素的另一個指標的依賴關(guān)系[18]. 二者計算方法[26-27]分別見式(3)(4).

(3)

(4)

式中，I為單變量莫蘭指數(shù)，為雙變量莫蘭指數(shù)，xi、xj分別為要素i、j的屬性值，為屬性值的平均值，是第二要素的平均值，wi,j分別為要素i和j之間的空間權(quán)重，n為要素總數(shù).

對計算得到的莫蘭指數(shù)，利用Z分布進行顯著性檢驗，檢驗公式：

(5)

(6)

VarI=E(I2)-E2(I)

(7)

式中，E(I)為莫蘭指數(shù)的期望值，E(I2)為莫蘭指數(shù)方差的期望值. 當Z得分大于1.96或小于-1.96時，說明該要素在95%置信區(qū)間內(nèi)呈現(xiàn)明顯的聚集或離散特征；當Z得分介于-1.96～1.96之間時，則說明該要素在95%置信區(qū)間內(nèi)呈隨機分布.

2 結(jié)果與討論

2.1 土壤重金屬含量描述

研究區(qū)土壤中5種重金屬元素含量的描述性統(tǒng)計結(jié)果(見表1)顯示，w(Cd)、w(Hg)、w(As)、w(Pb)的平均值均超過了湖南省土壤環(huán)境背景值，分別是湖南省土壤環(huán)境背景值的15.0、3.2、4.6、7.3倍，表明這4種重金屬元素含量的積累受人為活動的影響，其中，w(Cd)、w(As)、w(Pb)的平均值均超過GB 15618—2018《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風險管控標準(試行)》農(nóng)用地土壤污染風險篩選值，說明可能存在土壤污染風險；w(Cr)的平均值和中值均未超過湖南省土壤環(huán)境背景值，表明Cr未受到人為活動影響或影響較小. 該研究區(qū)5種重金屬元素含量的變異系數(shù)差異較大，數(shù)據(jù)離散程度較高，表明5種重金屬元素含量的空間差異比較明顯.

表1 研究區(qū)土壤中各重金屬元素含量

2.2 相關(guān)性分析

將5種重金屬元素含量正態(tài)化處理后，采用SPSS 19.0對蘇仙區(qū)土壤中5種重金屬元素含量進行Pearson相關(guān)性分析，結(jié)果(見表2)顯示，w(Cd)、w(As)、w(Pb)兩兩之間均呈顯著正相關(guān)(P<0.001)，且Pearson相關(guān)系數(shù)均在0.7以上，表明這3種元素間同源性很強.w(Hg)、w(Pb)之間呈顯著正相關(guān)(P<0.001)，但Pearson相關(guān)系數(shù)較小，說明這兩種元素亦存在部分同源性.w(As)、w(Cr)之間呈顯著負相關(guān)(P<0.05)，表明兩種元素間不存在同源性.

表2 各重金屬元素含量間的Pearson相關(guān)系數(shù)

2.3 土壤重金屬的來源解析

利用PMF 5.0對302個土壤樣品中5種重金屬元素的來源進行解析，并給出各因子的貢獻率. 初步將各重金屬元素數(shù)據(jù)載入后信噪比(S/N)均大于3，定義為“strong”. 運行模型后，多次調(diào)試重金屬元素的“strong”“weak”以及改變因子個數(shù)進行多次迭代計算，最終確定4個因子數(shù)，重金屬元素含量實測值與模擬值的擬合系數(shù)均大于0.75，實測Q值與理論Q值的偏差小于10%，解析結(jié)果較好.

圖2 研究區(qū)土壤重金屬污染源貢獻率

由圖2可見，因子1對5種重金屬元素的貢獻率各不相同，其中對Cr的貢獻率最高，達到88.2%. 該研究區(qū)w(Cr)平均值低于湖南省土壤環(huán)境背景值，受人為活動影響較小，因此，將因子1識別為自然源，受成土母質(zhì)影響. 因子2對Hg的貢獻率最高，為80.6%，其次是Pb，貢獻率為24.2%，二者含量相關(guān)性較強，表明二者存在部分同源性. 有色金屬冶煉、燃煤、金礦和汞礦活動是我國最主要的汞排放源[29]. Hg和Pb的環(huán)境遷徙性較強，多項研究指出大氣沉降是土壤中Hg、Pb元素的主要來源[30-32]，因此，將因子2識別為工業(yè)源中的大氣干濕沉降源. 因子3對Cd和Pb兩種元素的貢獻率較高，分別為63.2%和65.9%；Cd和Pb含量相關(guān)性較強，顯示同源性較強. 金屬礦山的開采、冶煉、重金屬尾礦、冶煉廢渣和礦渣堆放等是土壤中Cd、Pb污染的主要來源[33-35]. 因此，將因子3識別為工業(yè)生產(chǎn)過程中直接排放的工業(yè)廢棄物. 因子4對As的貢獻率最大，為75.1%. 人為As來源眾多，且與Cd、Pb、Hg等元素來源相似，含砷礦石的開采、運輸、加工等各環(huán)節(jié)都有損耗；另外，砷化合物作為原料的玻璃、顏料、原藥、紙張的生產(chǎn)以及煤的燃燒等過程，都可產(chǎn)生含砷廢水、廢氣和廢渣[36-37]，因此，將因子4識別為工業(yè)混合源.

綜上，該研究區(qū)表層土壤中5種重金屬元素的來源為自然源和3類工業(yè)源，其中，自然源的綜合貢獻率為41.60%，因子2、3、4三類工業(yè)源的綜合貢獻率分別為22.05%、20.05%和16.30%.

圖3 主要污染源貢獻率的空間分布

2.4 污染源貢獻率的空間分布

為探究各污染源的主要影響區(qū)域，將PMF模型解析得到的各源對單個樣品的貢獻率在ArcGIS 10.3下用普通克里金插值，得到各源貢獻率的空間分布情況(見圖3). 因子1代表自然源，該源貢獻率的高值點主要分布在蘇仙區(qū)北部，源貢獻率均超過52.75%，低值點主要分布在蘇仙區(qū)中部，說明蘇仙區(qū)北部受人為活動影響較小，中部受人為活動影響較大，污染程度較高. 整體來看，因子1對蘇仙區(qū)大部分區(qū)域的土壤樣品都有著較大的貢獻，貢獻率大多在32.39%以上，說明因子1作為自然源對土壤樣品的影響具有普遍性. 因子2代表工業(yè)源中的大氣干濕沉降源，該源貢獻率的高值點主要分布在蘇仙區(qū)中東部，貢獻率大多超過36.84%. 該高值區(qū)主要地形為山地，西部和北部鄰近區(qū)域分布大量的礦山開采和重金屬冶煉企業(yè)，地理位置特殊；整體來看，該污染源對大部分地區(qū)的影響較小，貢獻率大多低于24.62%，說明因子2作為遷徙源具有區(qū)域聚集性. 因子3代表直接排放的工業(yè)廢棄物，該源貢獻率的高值點主要分布在蘇仙區(qū)中部及南部少部分地區(qū)，說明這些區(qū)域受工礦企業(yè)影響較大，人為污染較為嚴重. 因子4貢獻率的高值點主要分布在蘇仙區(qū)中部及北部，整體來看，該源對大部分區(qū)域影響較小，貢獻率較低，大多在26.43%以下.

圖4 重點污染行業(yè)核密度分布

2.5 因子貢獻率與企業(yè)密度的空間相關(guān)性

按照GB/T 4754—2017《國民經(jīng)濟行業(yè)分類》將研究區(qū)內(nèi)342家污染企業(yè)分為四大類，其中采選業(yè)(08黑色金屬礦采選業(yè)、09有色金屬礦采選業(yè))215個，冶煉和壓延加工業(yè)(31黑色金屬冶煉和壓延加工業(yè)、32有色金屬冶煉和壓延加工業(yè))105個，化學原料和化學制品制造業(yè)(26化學原料和化學制品制造業(yè))12個，其他行業(yè)(廢棄資源綜合利用業(yè)、倉儲業(yè)、生態(tài)保護和環(huán)境治理業(yè)等)10個. 利用ArcGIS 10.3繪制全部企業(yè)及不同行業(yè)的核密度圖，分析各污染源與不同行業(yè)間的空間關(guān)系，明晰3類工業(yè)源的具體類別. 由圖4可見，全部企業(yè)密度高值區(qū)主要分布在蘇仙區(qū)中部及西南少部分地區(qū)，其中采選業(yè)主要聚集在蘇仙區(qū)中部，冶煉和壓延加工業(yè)主要聚集在蘇仙區(qū)中部、東北部及西南少部分地區(qū)，化學原料和化學制品制造業(yè)主要聚集在中北部，其他行業(yè)主要聚集在蘇仙區(qū)中部少部分地區(qū).

運用GeoDa軟件分別計算重點污染企業(yè)及各因子貢獻率的單變量莫蘭指數(shù)，結(jié)果(見表3)顯示，在P<0.001顯著性水平下，研究區(qū)重點污染企業(yè)及各因子貢獻率的單變量莫蘭指數(shù)均為正，且數(shù)值較大，Z得分均大于1.96，說明不同行業(yè)及各因子貢獻率在區(qū)域內(nèi)的分布具有聚集性特征. 單變量莫蘭指數(shù)分析結(jié)果均通過顯著性檢驗，具有統(tǒng)計學意義，可為后續(xù)雙變量檢驗提供理論支撐.

為探討各源貢獻率與企業(yè)間的空間依賴關(guān)系，將各源貢獻率的空間分布結(jié)果與各行業(yè)的核密度結(jié)果進行雙變量莫蘭指數(shù)分析，結(jié)果(見表4)顯示，因子1在P<0.001顯著性水平下，與全部企業(yè)、采選業(yè)、化學原料和化學制品制造業(yè)的雙變量莫蘭指數(shù)為負，表明三者對因子1的貢獻存在負影響. 在重點污染企業(yè)分布密集區(qū)域，工業(yè)源貢獻率較高，自然源貢獻率較低. 因子1與冶煉和壓延加工業(yè)、其他行業(yè)間未通過顯著性水平(P>0.1)，因其各自單變量莫蘭指數(shù)均通過顯著性檢驗，說明要素本身在空間上呈現(xiàn)聚集特征，但因子1與二者間卻不存在空間上的依賴關(guān)系，呈空間隨機性. 所得結(jié)果進一步驗證了因子1作為自然源的結(jié)論.

表3 研究區(qū)重點污染企業(yè)及污染源貢獻率的單變量莫蘭指數(shù)

表4 研究區(qū)重點污染企業(yè)與污染源貢獻率間的雙變量莫蘭指數(shù)

因子2在P<0.001顯著性水平下，與全部企業(yè)和4類重點污染行業(yè)間的雙變量莫蘭指數(shù)均為負，表明在4類企業(yè)分布較密集的區(qū)域內(nèi)，因子2的貢獻率多為低值，說明因子2為非直接工業(yè)源. 因子2對Hg元素的貢獻最大，氣態(tài)單質(zhì)汞(GEM)、氣態(tài)氧化汞(GOM)和顆粒態(tài)汞(PBM)可以在大氣中傳輸，尤其氣態(tài)單質(zhì)汞可隨大氣運動進行跨界傳輸[38-39]. 宋正城[40]指出，有色金屬冶煉活動對周邊表層土壤中汞的影響最嚴重的區(qū)域并非位于冶煉廠區(qū)，而是位于冶煉廠數(shù)公里外，說明汞的積累具有空間遷移性. 所得結(jié)果進一步驗證了因子2作為非直接工業(yè)源中的大氣干濕沉降源的結(jié)論. 對于遷徙源，污染源與其影響區(qū)域不存在對應的空間位置關(guān)系，難以識別到具體行業(yè)中，其詳細成因還需進一步確認.

因子3在P<0.001顯著性水平下，與全部企業(yè)和采選業(yè)的雙變量莫蘭指數(shù)較大，表明因子3與重點污染企業(yè)中的采選業(yè)在空間上存在顯著正相關(guān)，在采選業(yè)分布密集的區(qū)域內(nèi)，因子3的貢獻率亦多為高值；因子3與其余3類行業(yè)均不具有正相關(guān)關(guān)系或相關(guān)性不明顯，所得結(jié)果進一步明確了因子3為采選業(yè)排放的工業(yè)廢棄物，該污染源與其貢獻率的高值區(qū)存在對應的空間位置關(guān)系，因此該類廢棄物指在空間上未經(jīng)遠距離主動或被動運輸轉(zhuǎn)移，直接對排放口周圍的土壤造成污染.

因子4在P<0.001顯著性水平下，與4類重點污染行業(yè)間的雙變量莫蘭指數(shù)均為正，其中與化學原料和化學制品制造業(yè)的雙變量莫蘭指數(shù)最大，表明4類重點污染行業(yè)對因子4的貢獻率均有正向影響，在重點污染企業(yè)分布密集的區(qū)域內(nèi)，因子4的貢獻率亦多為高值，進一步說明因子4為各行業(yè)的工業(yè)混合源，且以化學原料和化學制品制造業(yè)排放的工業(yè)廢棄物為主.

該研究首先對研究區(qū)土壤中5種重金屬元素進行污染源大類的定性識別，并計算出每個污染源的貢獻率，然后依據(jù)污染源貢獻率與企業(yè)密度間的空間關(guān)系，將工業(yè)源細分到具體的重點污染行業(yè)上，對于存在相似貢獻的工業(yè)源，亦有很好地表征和識別. 秦治恒等[41]利用空間自相關(guān)方法，發(fā)現(xiàn)湘江流域中土壤Cd含量分布與污染企業(yè)的位置存在很強的相關(guān)性. 于靖靖等[18]基于莫蘭指數(shù)和廣義加性模型發(fā)現(xiàn)包含蘇仙區(qū)在內(nèi)的湘江子流域重點污染企業(yè)影響區(qū)中，采選業(yè)對土壤中Cd的影響最大. 這兩項研究都利用了空間自相關(guān)理論對特定污染源的定性識別，但尚未對污染源的貢獻量進行探討，陳雅麗等[42]統(tǒng)計了近10年關(guān)于土壤重金屬污染源解析的研究成果，統(tǒng)計結(jié)果指出，對湖南省而言，工業(yè)活動及其產(chǎn)生的大氣沉降是土壤中重金屬累積的主要來源. 其中，在礦產(chǎn)業(yè)發(fā)達的湖南省，Cd主要受采礦業(yè)的影響，Hg主要受工業(yè)活動產(chǎn)生的大氣沉降影響，Pb一方面來自工業(yè)活動，另一方面是含鉛灰塵的沉降進入表層土，As在工業(yè)發(fā)達地區(qū)受土壤母質(zhì)和工業(yè)活動的影響，Cr主要受土壤母質(zhì)影響. 筆者所得結(jié)果與上述研究結(jié)果基本一致. 對比可知，PMF模型與莫蘭指數(shù)結(jié)合的方法，一方面可以定量解析出各污染源的影響程度及主要影響區(qū)域；另一方面可以達到明確相關(guān)重點污染企業(yè)對土壤中5種重金屬元素的影響，從而有針對性地提出區(qū)域企業(yè)管理、重點污染源監(jiān)管及土壤環(huán)境治理策略.

3 結(jié)論

a) 研究區(qū)5種重金屬元素受4類源影響. 其中，自然源是Cr的主要來源，對Cr的貢獻率為88.2%；大氣干濕沉降是Hg的主要來源，對Hg的貢獻率為80.6%；采選業(yè)排放的工業(yè)廢棄物是Cd和Pb的主要來源，對Cd和Pb的貢獻率分別為63.2%、65.9%；以化學原料和化學制品制造業(yè)排放的廢棄物為主的工業(yè)混合源是As的主要來源，對As的貢獻率為75.1%.

b) 研究區(qū)中部受重點污染企業(yè)影響嚴重，且主要以采選業(yè)排放的廢棄物為主；此外，中部偏北地區(qū)主要以化學原料和化學制品制造業(yè)廢棄物為主；中部偏東地區(qū)受大氣干濕沉降影響嚴重，該地區(qū)連接3條支流上游，應警惕河流輸送和污水灌溉等活動引起的面源污染；西北部地區(qū)受重點污染企業(yè)影響較小.

c) PMF模型與莫蘭指數(shù)相結(jié)合的方法，使解析結(jié)果從工業(yè)源大類明晰到不同重點污染行業(yè)上，不僅能更好地識別有相似貢獻的工業(yè)源，而且可以輔助優(yōu)化驗證PMF的有效性.