基于受體模型與地統(tǒng)計的耕地土壤重金屬污染源解析

呂柏楠, 王 超, 師華定, 李明秋

1.河南理工大學, 河南 焦作 454003 2.生態(tài)環(huán)境部土壤與農(nóng)業(yè)農(nóng)村生態(tài)環(huán)境監(jiān)管技術中心, 北京 100012

土壤是自然生態(tài)系統(tǒng)的重要組成部分，是人類賴以生存和農(nóng)業(yè)生產(chǎn)的重要物質(zhì)基礎[1]. 近年來，我國環(huán)境問題日益嚴重，工礦、農(nóng)業(yè)土壤環(huán)境狀況不容樂觀[2]. 2014年全國土壤污染狀況調(diào)查顯示，我國耕地土壤污染點位超標率達19.4%，主要污染物為重金屬，其中鎘污染點位超標率達到7.0%，其次分別為鎳、砷、銅、汞、鉛、鉻等[3]；據(jù)2016年原農(nóng)業(yè)部監(jiān)測系統(tǒng)調(diào)查數(shù)據(jù)，我國每年因為土壤重金屬污染而使糧食產(chǎn)量降低1.0×107t以上[4]. 重金屬污染物在土壤中不斷富集，不僅會造成耕地質(zhì)量下降、農(nóng)產(chǎn)品產(chǎn)量降低等負面影響，還會通過食物鏈進入人體，威脅人類健康[5]. 因此，明晰重金屬污染物的來源，是進行土壤重金屬污染防治的重要前提.

空間插值技術可以通過合理布設土壤點位，實現(xiàn)研究區(qū)土壤污染數(shù)據(jù)的全覆蓋. 因此，如何進行空間模擬預測及方法的選用是準確預測區(qū)域土壤屬性空間分布特征的關鍵[6]. 克里格(Kriging)插值為土壤屬性的空間預測提供了一種無偏最優(yōu)估值方法，是目前地統(tǒng)計學中應用最廣的方法之一[7-8]. 作為克里格的衍生方法，指示克里格是一種不依賴于正態(tài)分布假設條件的非參數(shù)估值方法，可以處理高程度的空間變異現(xiàn)象而不必削減重要的高值數(shù)據(jù)，適合用于解決存在多特異值的土壤污染數(shù)據(jù)[9-10]. 段藝芳等[11]采用指示克里格法，發(fā)現(xiàn)陜北地區(qū)8種土壤重金屬生態(tài)風險在空間上均表現(xiàn)出由南向北、由東南向西北逐漸增加的趨勢；楊奇勇等[12]采用指示克里格法，發(fā)現(xiàn)云南省廣南縣除Cu元素外，其他重金屬元素高概率污染區(qū)域都集中在巖溶區(qū).

污染源解析技術是針對污染物的來源進行定性(源識別)或定量(源解析)判斷的方法[13]. 目前，受體模型法是當前土壤污染物源解析研究中最主要的技術方法. 受體模型種類較多，常用的有化學質(zhì)量平衡(CMB)模型[14]、主成分分析/因子分析-多元線性回歸(PCA/FA-MLR)模型[15-16]、正定矩陣因子分解(PMF)模型[17]和UNMIX模型[18]等. 其中，主成分分析/因子分析-多元線性回歸(PCA/FA-MLR)模型在土壤重金屬源解析領域應用較早，它不需要準確的源成分譜數(shù)據(jù)，只需要知道受體信息及排放源組成，使用較為廣泛[19]. Li等[20]采用FA-MLR模型，揭示了重慶4個典型農(nóng)業(yè)縣的工業(yè)排放、燃料消耗和當?shù)剞r(nóng)業(yè)對8種重金屬的主要人為貢獻為64%～90%；Yang等[21]利用PCA-MLR模型，揭示了湖北省青山區(qū)東部地區(qū)2010—2014年期間增加的重金屬濃度有82.5%來自工業(yè)/農(nóng)業(yè)活動、17.5%來自交通源.

云南是礦產(chǎn)資源大省，素有“有色金屬王國”之稱，礦產(chǎn)資源礦種齊全、儲量豐富，金子河流域銅礦、鉛鋅礦、金礦等礦產(chǎn)資源豐富，開發(fā)歷史悠久，黑色金屬礦和有色金屬礦采選企業(yè)數(shù)量較多. 因此，該研究以位于金子河流域中游的銅廠鄉(xiāng)耕地為研究區(qū)域，運用地統(tǒng)計學方法和主成分分析-多元線性回歸(PCA-MLR)模型對研究區(qū)土壤中的Cd、Hg、As、Pb、Cr 5種重金屬的空間分布及來源進行解析，探尋金子河兩岸農(nóng)用地土壤重金屬污染成因，以期為該區(qū)域土壤污染防控提供科學依據(jù).

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于云南省紅河哈尼族彝族自治州金平縣西部(22°43′N～22°52′N、102°55′E～103°03′E)，如圖1所示，地勢南高北低，最高海拔 1 191 m，最低海拔413 m，年平均氣溫13.6 ℃，年均降雨量800～1 000 mm，主導風向為西南風；土壤類型以黃棕壤為主，成土母巖以灰?guī)r和玄武巖為主；研究區(qū)西岸區(qū)域以旱地為主，東岸區(qū)域以水田為主；主要重金屬污染企業(yè)有黑色金屬(鐵、銅)礦采選廠、有色金屬(金)礦采選廠以及有色金屬礦山尾礦庫，主要分布在研究區(qū)的北部和西部區(qū)域.

圖1 金子河流域礦區(qū)分布及研究區(qū)采樣點分布Fig.1 Distribution of mining areas in the Jinzi River Basin and distribution of sampling points in the study area

1.2 樣品的采集與分析

土壤采樣點采用網(wǎng)格布點法進行布設[22]，網(wǎng)格大小為250 m×250 m，采樣時用GPS進行精確定位，采集表層(0～20 cm)土壤樣品76個，并選取其中8個點位采集深層(100～150 cm)土壤樣品. 采用梅花形布點法采集4～6個子樣，每個子樣均使用竹鏟垂直采集，去除雜草、礫石、肥料塊等雜質(zhì)，等份組合成一個混合樣. 室內(nèi)將樣品進一步用球磨機碾碎過200目尼龍篩，土壤重金屬全量采用HNO3-HCl-HF-HClO4溶樣，采用電感耦合等離子體原子發(fā)射光譜法(ICP-AES)測定土壤中Cd、Cr、Pb的含量，采用原子熒光法測定土壤中的As、Hg的含量[23].

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

研究區(qū)土壤中5種重金屬元素的描述性統(tǒng)計分析與樣品污染源解析均借助SPSS 19.0軟件完成；半方差函數(shù)借助GS+9.0軟件完成；指示克里格插值借助ArcGIS 10.2軟件完成.

1.4 研究方法

1.4.1內(nèi)梅羅綜合污染指數(shù)法

內(nèi)梅羅綜合污染指數(shù)(Nemerow integrated pollution index，)由美國學者Nemerow[24]提出，是在綜合各因子平均污染狀況的基礎上，兼顧極值和均值，有效評價研究區(qū)各采樣點重金屬的綜合污染狀況，是目前進行綜合污染指數(shù)計算最常用的方法之一[25]，計算公式：

(1)

式中，PN為內(nèi)梅羅綜合污染指數(shù)，Pimax為采樣點重金屬i污染物單項污染指數(shù)中的最大值，Piave為采樣點重金屬污染物單項污染指數(shù)的平均值.

內(nèi)梅羅指數(shù)評價法分級標準見表1.

表1 內(nèi)梅羅綜合污染指數(shù)評價法分級標準

1.4.2指示克里格法

指示克里格(indicator kriging)方法是由Journel提出的，被廣泛應用于具有極端值和高偏離度特性的土壤污染研究[26-29]. 指示克里格預測值的大小代表大于或小于特定閾值的概率，即從指示數(shù)據(jù)得出的非抽樣地點的期望值等于變量的累計分布[30]. 指示克里格具體計算方法[31]：①根據(jù)指標的閾值確定指示器；②計算最佳指示半方差函數(shù)；③利用半方差函數(shù)模型參數(shù)在ArcGIS 10.2軟件中進行指示克里格插值，得到土壤重金屬污染概率的空間分布圖.

1.4.3PCA-MLR模型

主成分分析-多元線性回歸法(PCA-MLR)首先利用主成分分析法定性識別污染源，在獲得污染因子數(shù)量后，對各示蹤元素的濃度進行多元線性回歸(MLR)，就可得到重金屬濃度與各示蹤元素濃度的回歸式，回歸系數(shù)用于計算各源的相對貢獻率[32-33].

2 結果與討論

2.1 土壤重金屬含量描述

由土壤重金屬含量的描述性統(tǒng)計分析結果(見表2)可知，研究區(qū)土壤中Cd、Hg、As、Pb、Cr的平均含量分別為3.17、0.16、44.34、779.57和358.85 mg/kg，均高于云南省土壤背景值. 其中，Pb的平均含量(779.57 mg/kg)遠高于云南省土壤背景值(40.60 mg/kg)，前者約是后者的19.20倍；其次，Cd的平均含量(3.17 mg/kg)約是云南土壤背景值(0.22 mg/kg)的14.55倍，表明這2種重金屬在土壤中的富集程度較高，可能受到的人為活動較劇烈. 依據(jù)GB 15618—2018《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風險管控標準(試行)》可知，研究區(qū)土壤中Cd、Hg、As、Pb、Cr含量低于篩選值的點位比例分別為26.32%、100%、56.58%、39.47%、17.11%，超過管制值的點位比例分別為40.79%、0、5.26%、40.79%、2.63%，說明研究區(qū)Cd、Pb污染較為嚴重. 從變異系數(shù)看，Pb、Cd的變異系數(shù)最大，分別為148%、139%，屬于高程度變異；其次為As(95%)、Hg(76%)、Cr(55%)，為中等程度變異. 其中，Pb、Cd在研究區(qū)分布極不平衡，表明土壤中Pb、Cd的含量受到了較為明顯的人為影響.

表2 樣品描述性統(tǒng)計

2.2 內(nèi)梅羅綜合污染指數(shù)評價分析

內(nèi)梅羅綜合污染指數(shù)反映的研究區(qū)土壤點位重金屬的綜合污染情況如表3所示. 研究區(qū)綜合污染指數(shù)范圍為1.73～131.69，平均值為18.28，處于重度污染程度的土壤點位比例為90.79%，說明研究區(qū)土壤整體處于重度污染水平. 其中，土壤中Cd、Hg、As、Pb、Cr的單因子污染指數(shù)范圍分別為0.14～96.75、0.62～13.50、0.07～7.07、0.23～178.77、1.45～19.59，而土壤點位中Cd、Pb、Cr屬于重度污染的比例分別為63.16%、59.21%和78.95%，表明研究區(qū)Cd、Pb和Cr污染比較嚴重.

表3 土壤重金屬內(nèi)梅羅綜合污染指數(shù)

2.3 空間分布特征

選擇GB 15618—2018的篩選值作為5種重金屬元素的閾值，研究區(qū)土壤樣點中的Hg含量均低于篩選值，對農(nóng)產(chǎn)品質(zhì)量安全、農(nóng)作物生長或土壤生態(tài)環(huán)境的風險低，因此不進行空間插值分析.

2.3.1半方差函數(shù)分析

對各元素變量進行半方差函數(shù)模型擬合，分析局部地區(qū)變量的空間變異及相關程度. 擬合結果顯示：Cd、As、Pb這3種重金屬的決定系數(shù)均大于0.9，所選擬合模型均符合地統(tǒng)計要求，而Cr的決定系數(shù)僅為0.09，說明Cr元素在研究區(qū)范圍內(nèi)不具有空間相關性，不能進行半方差函數(shù)分析，不適宜空間局部插值. 塊基比能夠反映隨機因素導致的空間變異占總體變異的比例[35]，當塊基比值小于0.25時表示變量具有顯著的空間自相關，塊基比值范圍在0.25～0.75時表示變量具有中等自相關，塊基比值大于0.75時表示變量自相關程度不強[36]. 變程大小反映了當?shù)刂亟饘倏臻g異質(zhì)性的尺度或空間自相關尺度[37].

基于GS+ 9.0軟件，獲得研究區(qū)Cd、As、Pb的最佳指示半方差函數(shù)模型(見圖2)及參數(shù)(見表4). Cd、As的半方差函數(shù)符合球形模型，Pb符合高斯模型. Cd、As、Pb的塊基比分別為0.192、0.113、0.066，均在0.25以下，都表現(xiàn)為顯著的空間自相關性，說明在當前閾值條件下，Cd、As、Pb條件概率的空間分布是由結構性因素(如土壤母質(zhì)、地形、土壤類型等)共同作用引起的. As的變程為 2 523.88 m，Cd的變程為 1 723.00 m，Pb的變程為 1 190.04 m. 從變程數(shù)據(jù)來看，研究區(qū)域土壤中Cd、As空間分布的均一性相比Pb要強，而Pb在土壤中的分布存在小范圍的變異現(xiàn)象[38].

表4 土壤重金屬含量的理論模型及相關參數(shù)

圖2 指示半方差函數(shù)模型Fig.2 Indicating semivariance function model

2.3.2指示克里格分析

根據(jù)Cd、As、Pb的指示半方差函數(shù)模型，利用ArcGIS 10.2軟件進行克里格插值，得到研究區(qū)耕地土壤Cd含量大于0.3 mg/kg、As含量大于40 mg/kg、Pb含量大于90 mg/kg的污染概率分布結果(見圖3).

從圖3可以看出，研究區(qū)西岸土壤中Cd、As、Pb的污染風險較高，特別是Cd和Pb，在絕大部分區(qū)域的污染風險概率高于0.7，且Cd和Pb的污染風險概率分布相似，可能存在相同的污染源；而As僅在北部區(qū)域的污染風險概率低于0.3，在中南部絕大部分區(qū)域的污染風險概率高于0.7，而且與Cd和Pb的污染風險概率分布相似，三者可能存在相同的污染源. 而西岸地區(qū)屬于海拔為500～1 000 m的山地區(qū)域，富含礦產(chǎn)資源，且附近存在多家金屬礦選企業(yè). 因此，推斷西岸重金屬的主要污染來自礦產(chǎn)資源的開發(fā)及金屬冶煉所產(chǎn)生的廢水、廢氣、廢渣[39-41].

圖3 表層土壤Cd、As和Pb污染概率分布Fig.3 Pollution probability distribution of Cd, As and Pb in top layer soils

研究區(qū)的東岸地區(qū)，土壤Cd、As、Pb的污染風險較低，其中，Cd的污染風險概率僅在北部區(qū)域高于0.7，中南部絕大部分區(qū)域低于0.3；As的污染風險概率在東岸絕大部分區(qū)域低于0.3；Pb僅在北部小部分區(qū)域存在概率高于0.7的污染風險，其余大部分區(qū)域的污染風險概率低于0.3. 而東岸附近無企業(yè)，因此，推斷東岸重金屬的主要污染來自土壤母質(zhì)或土壤類型等自然源.

總體而言，Cd、As、Pb污染的高概率區(qū)域主要分布在研究區(qū)的西部與西南部地區(qū)，北部地區(qū)是Cd、Pb污染的高概率區(qū)域，而Cd、As、Pb污染的低概率區(qū)域主要分布在研究區(qū)的東部以及東南部地區(qū)，Cd、As、Pb污染風險分布具有顯著的地域性差異.

2.4 土壤重金屬的來源解析

根據(jù)研究區(qū)重金屬元素的空間分布特征，將研究區(qū)的東西兩岸(金子河為界)及研究區(qū)整體分別利用PCA-MLR模型進行來源解析，探究局部(東西兩岸)與整體(研究區(qū))污染源的關系.

2.4.1源類型識別結果

由表5可見，對于研究區(qū)整體而言，PCA-MLR源1主要載荷元素為Pb、Cd、As和Hg. 金屬礦山的開采冶煉、工業(yè)廢渣和礦渣堆放等是土壤中Pb、Cd污染的主要來源[42-43]；As常見于多種金屬礦床中，尤其多伴生于金礦床[44]；金屬冶煉工藝是環(huán)境中汞元素的主要來源[45]. 而依據(jù)實地調(diào)查可知，研究區(qū)的北部地區(qū)存在有色金屬礦(金礦)采選廠及尾礦庫，西部地區(qū)存在黑色金屬礦(鐵礦、銅礦)采選廠和有色金屬(金礦、鉛鋅礦)采選廠. 因此，將因子1識別為以選冶為主的工業(yè)源. 對于東岸和西岸而言，PCA-MLR源1主要載荷元素均為Cd、As和Pb. 其中，東岸的源成分譜中As占比較大，因此，將東岸因子1識別為以金礦選冶為主的工業(yè)源；西岸的源成分譜中As所占比重較小，因此，將西岸因子1識別為以黑色金屬礦(鐵礦、銅礦)和鉛鋅礦選冶為主的工業(yè)源.

表5 土壤樣品的主成分因子載荷

對于研究區(qū)整體、東岸及西岸而言，PCA-MLR源2主要載荷元素均為Cr、Hg. 該研究區(qū)Cr的變異系數(shù)較小，受人為活動影響較小. 由于深層土壤中的重金屬含量高于表層土壤，可能受土壤成土母質(zhì)的影響[46]. 由圖4可見，研究區(qū)大部分表層土壤中Cr的含量低于深層土壤，因此，將因子2識別為自然源，受成土母質(zhì)影響.

圖4 研究區(qū)Cr元素表層和深層土濃度對比Fig.4 Contrast of Cr element surface and deep soil concentration in the study area

2.4.2源貢獻率結果

由圖5可見，研究區(qū)整體、東岸和西岸的土壤重金屬來源均為2種，分別是工業(yè)源和自然源. 對于研究區(qū)整體而言，工業(yè)源的相對貢獻率最大，為87.21%；對于東岸而言，自然源的相對貢獻率最大，為93.46%；對于西岸而言，工業(yè)源的相對貢獻率最大，為91.02%. 因此，成土母質(zhì)等自然因素是東岸土壤重金屬含量的主要來源；工業(yè)礦選冶煉是西岸土壤重金屬含量的主要來源，也是整個研究區(qū)土壤重金屬含量的主要來源.

圖5 研究區(qū)整體、東岸和西岸表土中 重金屬各類源的相對貢獻率Fig.5 The relative contribution rate of various sources of heavy metals in the surface soils of the whole study area, the east bank and the west bank

造成東西岸區(qū)域土壤重金屬污染主要來源不同的原因，主要是西岸區(qū)域周圍較多的礦產(chǎn)資源開采冶煉過程中產(chǎn)生的廢水、煙粉塵及冶煉渣、尾礦庫尾砂中的殘余重金屬淋濾遷移，以及關閉或搬遷企業(yè)歷史遺留重金屬廢渣常年隨雨水沖刷、淋溶，進入周邊土壤. Xu等[47]研究發(fā)現(xiàn)，礦業(yè)對Cd、Pb、As的主要影響發(fā)生在周邊地區(qū)5 km范圍內(nèi)，其中，西岸區(qū)域距離礦產(chǎn)開采冶煉場所基本在7 km范圍內(nèi)，而東岸區(qū)域距離礦產(chǎn)開采冶煉場所基本超過7 km. 尾礦庫作為礦山開采主要污染區(qū)域，尾礦庫中的尾礦砂隨風力遷移，及廢棄礦山構筑物淋濾水力作用驅(qū)動，而影響周圍環(huán)境[48]. 西岸區(qū)域海拔高且落差大，尾礦砂中的重金屬受降水等水力作用影響范圍較廣，而東岸距離尾礦庫較遠，地勢較為平緩，受到來自尾礦庫的影響較小. 因此，要控制研究區(qū)耕地周圍重點污染企業(yè)“三廢”的排放，加強對關停或搬遷的遺留工業(yè)場地周邊污染的監(jiān)測，才能有效防范工業(yè)產(chǎn)生的重金屬污染對周邊耕地的影響.

3 結論

a) 金子河流域土壤中Cd、Hg、As、Pb、Cr含量低于GB 15618—2018《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風險管控標準(試行)》的篩選值的點位比例分別為26.32%、100%、56.58%、39.47%、17.11%，含量超過管制值的點位比例分別為40.79%、0%、5.26%、40.79%、2.63%；研究區(qū)域內(nèi)重金屬綜合污染指數(shù)較高，處于重度污染程度的土壤點位比例為90.79%，且Cd、Pb和Cr的污染比較嚴重.

b) 從空間分布特征可以看出，元素Cd、As、Pb污染的高概率區(qū)域主要分布在金子河流域西部與西南部地區(qū)，Cd、Pb污染的高概率區(qū)域主要分布在研究區(qū)的北部地區(qū)，Cd、As、Pb污染風險分布具有顯著的地域性差異.

c) 金子河流域東西兩岸表土中重金屬污染的主要來源不同. 成土母質(zhì)等自然因素是東岸土壤重金屬含量的主要來源；工業(yè)礦選冶煉是西岸土壤重金屬含量的主要來源，也是整個研究區(qū)土壤重金屬含量的主要來源.