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

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

1.河南理工大學, 河南 焦作 454003 2.生態環境部土壤與農業農村生態環境監管技術中心, 北京 100012

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

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

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

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

1 材料與方法

1.1 研究區概況

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

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

1.2 樣品的采集與分析

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

1.3 數據處理

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

1.4 研究方法

1.4.1內梅羅綜合污染指數法

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

(1)

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

內梅羅指數評價法分級標準見表1.

表1 內梅羅綜合污染指數評價法分級標準

1.4.2指示克里格法

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

1.4.3PCA-MLR模型

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

2 結果與討論

2.1 土壤重金屬含量描述

由土壤重金屬含量的描述性統計分析結果(見表2)可知，研究區土壤中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種重金屬在土壤中的富集程度較高，可能受到的人為活動較劇烈. 依據GB 15618—2018《土壤環境質量 農用地土壤污染風險管控標準(試行)》可知，研究區土壤中Cd、Hg、As、Pb、Cr含量低于篩選值的點位比例分別為26.32%、100%、56.58%、39.47%、17.11%，超過管制值的點位比例分別為40.79%、0、5.26%、40.79%、2.63%，說明研究區Cd、Pb污染較為嚴重. 從變異系數看，Pb、Cd的變異系數最大，分別為148%、139%，屬于高程度變異；其次為As(95%)、Hg(76%)、Cr(55%)，為中等程度變異. 其中，Pb、Cd在研究區分布極不平衡，表明土壤中Pb、Cd的含量受到了較為明顯的人為影響.

表2 樣品描述性統計

2.2 內梅羅綜合污染指數評價分析

內梅羅綜合污染指數反映的研究區土壤點位重金屬的綜合污染情況如表3所示. 研究區綜合污染指數范圍為1.73～131.69，平均值為18.28，處于重度污染程度的土壤點位比例為90.79%，說明研究區土壤整體處于重度污染水平. 其中，土壤中Cd、Hg、As、Pb、Cr的單因子污染指數范圍分別為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%，表明研究區Cd、Pb和Cr污染比較嚴重.

表3 土壤重金屬內梅羅綜合污染指數

2.3 空間分布特征

選擇GB 15618—2018的篩選值作為5種重金屬元素的閾值，研究區土壤樣點中的Hg含量均低于篩選值，對農產品質量安全、農作物生長或土壤生態環境的風險低，因此不進行空間插值分析.

2.3.1半方差函數分析

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

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

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

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

2.3.2指示克里格分析

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

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

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

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

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

2.4 土壤重金屬的來源解析

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

2.4.1源類型識別結果

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

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

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

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

2.4.2源貢獻率結果

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

圖5 研究區整體、東岸和西岸表土中 重金屬各類源的相對貢獻率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

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

3 結論

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

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

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