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

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

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

2.生態環境部土壤與農業農村生態環境監管技術中心， 北京 100012

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

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

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

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

1 材料與方法

1.1 研究區概況

研究區域位于湖南省郴州市蘇仙區(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. 該區域地處郴州市西河流域，西河是郴江的一級支流、湘江的三級支流. 蘇仙區礦產資源豐富，素有“有色金屬之鄉”的美譽，有世界罕見的柿竹園多金屬礦，國有硚口鉛鋅礦，瑪瑙山錳礦，東波鉛鋅礦以及年產煤數萬噸的許家洞、街洞和棲鳳渡煤礦等. 發達的礦業生產給當地帶來了一系列環境問題.

1.2 數據來源與預處理

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

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

1.3 數據處理

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

1.4 研究方法

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

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

(1)

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

PMF模型基于加權最小二乘法進行限定和迭代計算，利用樣品的重金屬濃度和不確定度數據進行各樣點的加權計算，使得目標函數Q最小化. 目標函數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莫蘭指數

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

(3)

(4)

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

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

(5)

(6)

VarI=E(I2)-E2(I)

(7)

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

2 結果與討論

2.1 土壤重金屬含量描述

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

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

2.2 相關性分析

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

表2 各重金屬元素含量間的Pearson相關系數

2.3 土壤重金屬的來源解析

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

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

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

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

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

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

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

圖4 重點污染行業核密度分布

2.5 因子貢獻率與企業密度的空間相關性

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

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

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

表3 研究區重點污染企業及污染源貢獻率的單變量莫蘭指數

表4 研究區重點污染企業與污染源貢獻率間的雙變量莫蘭指數

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

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

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

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

3 結論

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

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

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