基于地統計學的土壤重金屬分布與污染風險評價

楊奇勇 謝運球 羅為群 谷佳慧，2 曾紅春，2

(1.中國地質科學院巖溶地質研究所巖溶生態系統與石漠化治理重點實驗室， 桂林 541004；2.中國地質大學(北京)地球科學與資源學院， 北京 100083)

引言

土壤是農業生產的基礎，也是生態環境的重要組成部分[1]。它是環境中重金屬遷移、轉化的重要媒介，既是重金屬聚集的匯，也是向其他系統遷移的源[2]。土壤對污染物有一定的自凈能力，但是這種自凈能力是有一定限度的。隨著工農業發展、礦山開采等，人類向土壤排放的污染物越來越多，當進入土壤中的污染物超過土壤承受能力時，土壤就會向外界輸出污染物，從而危害人類健康[3-4]。因此，研究土壤重金屬元素含量的空間分布并對其污染風險進行定量評估有著重要的意義。

地統計學方法在環境污染物的空間分布和風險評估研究方面得到了越來越廣泛的應用[5-7]。普通克里格法(Ordinary Kriging, OK)能夠對未采樣區域實現最佳無偏性估計，是空間預測中最常用的方法。但是，在污染評估與風險管理的實際應用中，人們感興趣的可能并不是空間某一位置污染物含量的具體數值，而是其超出或低于某一閾值的風險概率分布。指示克里格法(Indicator Kriging，IK)能給出單一變量在一定閾值限制下的風險概率，多變量指示克里格法(Multiple variable indicator Kriging，MVIK)能夠將給定閾值的多個變量綜合成一個變量，并繪出綜合風險概率分布圖[8-11]，因而在環境風險評估中有較好的應用前景。如GOOVAERTS等[8]應用IK 法繪制了Swiss Jura地區土壤Cd、Cu和Pb等重金屬元素的概率分布圖，對該區域的土壤污染進行了風險評估; CHU等[9]利用臺灣中部的1082個土壤樣品，采用MVIK法對重金屬Cr、Cu、Ni和Zn進行土壤污染綜合風險評價，指出在給定的重金屬閾值下，調查區域土壤具有高的污染風險概率，風險概率高達0.82；LIN等[12]應用IK法對臺灣Chunghua縣的Cd、Cu、Ni、Zn等土壤重金屬元素進行單一重金屬污染風險評價，并用MVIK法對土壤重金屬污染進行了綜合風險評價。

本文選取云南省廣南縣廣南圖幅(1∶50 000)的野外調查數據，利用OK法進行土壤As、Cd、Cr、Cu、Pb、Hg、Ni、Zn等重金屬含量的空間變異研究，采用GB15618—1995《土壤環境質量標準》中的二級評價標準，利用IK法和MVIK法對各項土壤重金屬引起的土壤污染風險和污染綜合風險進行評估，以期為當地礦山開采、農業生產和環境保護等提供科學依據。

1 材料與方法

1.1 采樣區域概況

采樣區位于云南省廣南縣內，總面積約為402.91 km2，其中巖溶面積約為197.52 km2，占采樣區域總面積的49.02%(圖1)。采樣區屬于亞熱帶季風氣候，年平均氣溫為16.7℃，年平均降水量為1 042.1 mm。采樣區西南地勢低，為廣南盆地，其余地區為山地、丘陵。海拔處在971～1 701 m之間。采樣區內土壤類型主要是酸性母巖紅壤、酸性母巖黃壤、棕泥土和黃紫泥土。主要土地利用類型有林地、水田、旱地。林地主要分布在北部山地和東部的巖溶丘陵區，水田主要分布在西部，旱地主要分布在東部巖溶洼地。巖溶區內存在石漠化現象，且東部石漠化較嚴重。區域內有礦石開采與采煤等活動。

圖1 研究區位置與采樣點分布Fig.1 Location of study area and distribution of sample sites in study area

1.2 樣品采集與處理

按照DZ/T 0258—2014《多目標地球化學調查規范(1∶250 000)》要求進行間距為1 km的網格劃分，網格內預布點充分考慮土地利用的代表性，預布點布置在土壤易于匯集的平緩地、洼地、山間平壩地等部位，相鄰樣點距離大于500 m。采樣時間為2016年,采集0～20 cm深度的土柱，共獲取土壤樣品410件(圖1)。每個樣點由1個主樣點和2個分樣點組成，分樣點與主樣點距離超過25 m。利用手持GPS重新記錄主樣點位置作為采樣點坐標。

樣品加工、處理嚴格執行DZ/T 0258—2014《多目標地球化學調查規范(1∶250 000)》。土壤樣品風干后，過20目篩。按照4 km21個點進行樣品組合分析，共得到組合樣品102個，對組合樣品進行土壤重金屬含量測試，采用密碼插入法，按照分析樣品數插入一定數量的準確度和精密度一級標準控制物監控測試，確保樣品分析的準確度和精密度。重金屬元素和pH值的分析方法及檢出限見表1。

1.3 數據處理與分析

相關分析、經典統計分析在SPSS 16.0中進行。重金屬含量和重金屬指示變換值的半方差函數模擬在GS+7.0中完成，專題圖的制作在ArcGIS 9.2中完成。地統計分析的方法、原理詳見文獻[13]。

土壤污染風險指示克里格分析方法的主要步驟如下[14]：①根據問題的要求確定各項土壤重金屬污染評價的閾值。土壤重金屬污染閾值采用GB15618—1995《土壤環境質量標準》中的二級評價標準(表2)對采樣區8種重金屬進行污染評價。②確定土壤重金屬污染指標的指示函數。根據土壤重金屬污染閾值，利用指示函數對相應的采樣點數據進行二態指示變換，得到各樣點的指示變換值Zij(1或0)，用來評價相應采樣點上的土壤重金屬污染狀況。③在GS+7.0 中利用土壤重金屬指示變換值進行空間變異分析，擬合得到最佳變異函數模型。④將變異函數模型參數輸入到ArcGIS 9.2 中進行OK插值分析，得到滿足相應閾值的土壤重金屬元素污染風險的空間概率分布圖，計算平均概率。⑤將多個土壤重金屬指標合并為一個綜合指標

表1 樣品分析方法Tab.1 Testing methods for soil samples

(1)

式中k——重金屬指標個數

i——采樣點編號

⑥利用綜合指數重復步驟③、④，獲取土壤重金屬污染綜合指數空間分布圖，計算綜合指數污染風險平均概率。

表2 土壤重金屬污染閾值Tab.2 Threshold of pollution for soil heavy metal mg/kg

1.4 空間評價精度驗證

由于組合樣點數量有限，研究采用交叉驗證的方法對空間評價精度進行驗證[15]。具體的做法是：①從所有樣點中抽取一個樣點作為驗證樣點，利用剩余樣點進行空間預測。②將抽取的樣點依次替換剩余樣品中的其他樣點作為驗證樣點，每次替換后都用剩余的樣點進行一次空間預測，直到所有的樣點都替換完畢。③選取平均絕對誤差(MME)和均方根誤差(RRMSE)[16]對空間預測精度進行評價

(2)

(3)

式中O——驗證數據P——預測數據

n——驗證數據的數目

2 結果與分析

2.1 土壤重金屬的描述性統計分析

進行調查區土壤重金屬含量統計分析，發現巖溶地貌區土壤重金屬含量與非巖溶地貌區土壤重金屬含量存在很大的差異性(表3)。從表3中可以看出，巖溶區土壤重金屬含量最大值(除Ni元素外)、平均值要顯著高于非巖溶地貌區，巖溶區pH值平均值也高于非巖溶地貌區。這與已有的研究結果一致[17-18]。參考表2中土壤重金屬污染閾值標準，從平均值來看非巖溶地貌區只有Cd、Cu元素存在污染，但是巖溶地貌區除了As、Cr、Pb元素外其他重金屬元素都存在污染。從偏度系數來看，所有采樣點的重金屬元素含量都不呈正態分布，但經過對數轉換后，各項重金屬元素基本呈正態分布，適合利用OK法進行空間預測。

2.2 土壤重金屬元素間的相關性

土壤重金屬元素與pH值相關分析結果(表4)表明8種重金元素兩兩之間都存在顯著的正相關性(p<0.01)，相關系數在0.457～0.885之間，其中Zn元素和Cd元素的相關系數最高，As元素和Cd元素的相關系數最低。Zn元素與Cd、Cr、Pb、Hg、Ni，Cr元素與Pb、Ni，Pb元素與Hg的相關系數都在0.8以上，說明它們在來源上具有很大相似性，As與Cd、Cu的相關系數在0.5以下，說明三者的來源具有復合關系[19-20]。pH值與各項重金屬元素之間相關系數在0.365～0.673之間，相關性顯著(p<0.01)。

2.3 土壤重金屬元素含量的空間分布

表5為土壤各項重金屬元素含量對數轉換后的半方差函數理論模型和參數。從表5可以看出，采樣區土壤重金屬元素含量主要選用球狀、高斯和指數3種模型進行半方差函數擬合，各模型決定系數都大于0.8，處在0.896～0.979之間，說明所選取的模型能較好地反映各項土壤重金屬元素含量的空間結構特[21]。塊金值與基臺值之比C0/(C0+C)通常作為衡量變量空間相關程度的尺度，若其值小于0.25，則為空間強相關；處在0.25～0.75為空間中等強度相關；大于0.75為弱相關[13]。從表5可以看出，除Cu元素含量(塊基比為0.322)呈現中等偏強的空間相關性外，其他重金屬元素含量都具有強烈的空間相關性。說明采樣區域土壤重金屬含量的空間異質性主要是由巖性、成土母質等內在因素引起的，而人類活動等隨機性因素影響較小。利用獲得的模型及其參數(表5)，在ArcGIS 9.2中分別對8種重金屬元素含量進行空間插值，得到它們的空間分布圖(圖2)。

表3 研究區域不同地貌類型土壤重金屬含量統計特征Tab.3 Characteristics of soil heavy metals contents for different landforms in study area

注：下標1、2、3分別表示非巖溶區、巖溶區和全部樣點。

表4 土壤重金屬元素間的相關系數Tab.4 Correlation coefficients of soil heavy metals

注：所有的相關性在p<0.01水平顯著(雙尾)。

表5 土壤重金屬含量半方差函數的理論模型和參數Tab.5 Theoretical models and parameters of soil heavymetals contents

從圖2可以看出，8種重金屬元素含量總體上都是東南部、北部巖溶區含量高，西南部非巖溶地貌區的盆地中重金屬元素含量低；但在細節上，8種重金屬含量在西北部、北部空間分布有顯著差異。As、Pb和Hg含量空間分布相似，Cd和Cr含量空間分布相似，Cu、Ni和Zn含量空間分布相似等。重金屬含量空間分布的相似程度與它們在統計上的相關性強弱是一致的。

2.4 土壤重金屬污染風險評價

利用表2中土壤重金屬污染閾值對各項土壤重金屬元素進行指示變換，發現采樣區的Pb元素含量的指示變換值都為0，說明采樣區土壤沒有重金屬Pb元素污染；而Cr元素的指示變換值也只有1個為1，其余都為0，說明采樣區基本可以忽略土壤重金屬Cr元素的污染；As、Cd、Cu、Hg、Ni和Zn等重金屬元素指示變換值為1的樣點數量依次為33、84、17、19、27、27個，也就是說這6種重金屬元素分別有32.35%、82.35%、16.67%、18.63%、26.47%和26.47%的樣點土壤有污染。因此，采樣區的土壤重金屬污染評價只考慮這6種重金屬元素。利用這6種重金屬元素的指示變換值分別進行OK空間插值，得到各項重金屬元素土壤污染概率空間分布圖(圖3)。

圖2 重金屬含量空間分布Fig.2 Spatial distribution maps for soil heavy metal contents

圖3 土壤重金屬污染概率空間分布Fig.3 Probability pollution distribution maps of each soil heavy metal

從圖3可以看出，除Cu元素外，其他重金屬元素高概率污染區域都集中在巖溶區，其中As和Hg污染的高概率區域基本與巖溶地貌區分布一致(圖1和圖3)。將Cu元素的概率分布圖與Cu元素含量分布圖進行比較，發現Cu元素含量高值分布的東南部并沒有出現大面積的土壤污染高概率區域。這一現象與采樣區域東南部的高pH值分布有關，與朱禮學[22]的研究具有相似的結論。Cd元素污染概率分布主要集中在0.7～1.0的高概率區間，而其余5種重金屬污染區域主要集中在0.1～0.3的低概率區域。統計分析結果表明As、Cd、Cu、Hg、Ni和Zn重金屬污染的平均概率依次為0.326、0.805、0.185、0.192、0.267、0.270。因此，在本研究采用的重金屬污染閾值下，調查區域土壤具有較高的Cd和As重金屬污染風險。

利用綜合指數變換值Zij進行空間插值，得到土壤重金屬污染綜合風險概率分布圖(圖4)。

圖4 土壤重金屬綜合污染概率空間分布圖Fig.4 Probability pollution distribution map of soil heavy metals

采樣區域重金屬污染綜合風險概率處在0.09～0.78之間(圖4)，其中0.09～0.3之間的低概率污染區域面積約為193.36 km2，占整個采樣區的47.99%，主要分布在采樣區域的西南部盆地到東北部丘陵一帶。重金屬污染綜合風險在0.4～0.78之間的高概率污染區域面積約為69.95 km2，占整個采樣區域的17.36%，主要分布在采樣區域的東南部巖溶區。采樣區域重金屬污染綜合風險平均概率為0.335，在本研究采用的重金屬污染閾值下，具有較高的土壤污染風險。大量的研究表明巖溶區土壤重金屬具有較高的背景值[17-18]，但是由于巖溶區富鈣偏堿的環境，卻很少有出現重金屬誘發的地方病報道[22]。因此，針對巖溶區高背景的重金屬值和富鈣偏堿的巖溶環境，研究確定巖溶地貌區土壤重金屬污染閾值有著實際的意義。

2.5 預測精度評估

經過交叉驗證法檢驗，各項重金屬元素含量的OK法預測、土壤污染風險概率的IK法及土壤重金屬污染綜合風險概率的MVIK法預測都具有較低的MME值和RRMSE值(表6)，表明各項預測都有較高的精度，因而評價結果是可信的。

表6 預測結果精度驗證Tab.6 Precision evaluation for spatial prediction

3 結論

(1)土壤重金屬之間均呈現極其顯著的相關關系，其中Zn元素與Cd、Cr、Pb、Ni，Cr元素與Pb、Ni，Pb元素與Hg的相關系數都在0.8以上，具有較高的同源性，而As與Cd、Cu的相關系數在0.5以下，說明三者在來源上具有復合性。

(2)除Cu元素具有中等偏強的空間結構特征外，其余各項重金屬元素都具有強烈的空間結構特征，其空間異質性主要是由巖性、成土母質等內在因素引起的。重金屬含量高值主要分布在東南部和北部的巖溶地貌區，土壤重金屬含量低值主要分布在西南部的非巖溶地貌區。

(3)指示克里格分析表明，采樣區域不存在Pb和Cr重金屬污染，而As、Cd、Cu、Hg、Ni和Zn等重金屬元素存在不同程度的污染風險，其中Cd污染風險最高，達到0.805。土壤重金屬綜合污染風險為0.335，高風險區域主要在東南部的巖溶地貌區，低風險區域主要分布在西南部的非巖溶地貌區。

1 鄭袁明，陳同斌，陳煌，等．北京市不同土地利用方式下土壤鉛的積累[J]. 地理學報，2005,60(5)：791-797.

ZHENG Y M, CHEN T B, CHEN H, et al. Lead accumulation in soils under different land use types in Beijing City[J]. Acta Geographica Sinica, 2005, 60(5): 791-797. (in Chinese)

2 周萍，文安邦，史忠林，等. 三峽庫區不同土地利用土壤重金屬分布特征與污染評價[J/OL]. 農業機械學報, 2017, 48(7)：207-213. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20170726&journal_id =jcsam. DOI:10.6041/j.issn.1000-1298.2017.07.026.

ZHOU P, WEN A B, SHI Z L, et al. Distribution characteristics and pollution evaluation of soil heavy metals of different land use types in Three Gorges Reservior Region[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(7)：207-213. (in Chinese)

3 廖海軍．北京市密云水庫上游土壤重金屬污染調查評價[J]. 中國環保產業，2007, 2(8): 23-26.

LIAO H J. Pollution evaluation of soil heavy metals in the upper Miyun reservoir of Beijing City[J]. Environmental Protection Industry of Protection Industry of China, 2007, 2(8): 23-26. (in Chinese)

4 WU G, KANG H, ZHANG X, et al. A critical review on the bio-removal of hazardous heavy metals from contaminated soils: issues, progress, eco-environmental concerns and opportunities[J]. Journal of Hazardous Materials, 2010, 174(1-3):1-8.

5 司涵，張展羽，呂夢醒，等. 小流域土壤氮磷空間變異特征分析[J/OL]. 農業機械學報, 2014, 45(3):90-96. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140316&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2014.03.016.

SI H, ZHANG Z Y, Lü M X, et al. Spatial variability of soil nitrogen and phosphorus in small water shed[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(3):90-96. (in Chinese)

6 GHADERMAZI J, SAYYAD G, MOHAMMADI J, et al. Spatial prediction of nitrate concentration in drinking water using pH as auxiliary Co-kriging variable[J]. Procedia Environmental Sciences, 2011(3):130-135.

7 JANG C S, LIU C W, LU K L, et al. Delimitation of arsenic-contaminated groundwater using risk-based indicator approaches around blackfoot disease hyperendemic areas of southern Taiwan[J]. Environmental Monitoring and Assessment, 2007, 134(1-3): 293-304.

8 GOOVAERTS P，WEBSTER R，DUBOIS J P. Assessing the risk of soil contamination in the Swiss Jura using indicator geostatistics[J]. Enviromental and Ecological Statistics，1997，4(1): 31-48.

9 CHU H J, LIN Y P, JANG C S, et al. Delineating the hazard zone of multiple soil pollutants by multivariate indicator Kriging and conditioned Latin hypercube sampling[J]. Geoderma, 2010, 158(3-4): 242-251.

10 JUANG K W, LEE D Y. Comparison of three non-parametric Kriging methods for delineating heavy-metal contaminated soils[J]. Journal of Environmental Quality, 2000, 29(1): 197-205.

11 BELKHIRI L, MOUNI L, NARANY T S, et al. Evaluation of potential health risk of heavy metals in groundwater using the integration of indicator Kriging and multivariate statistical methods[J]. Groundwater for Sustainable Development, 2017, 4:12-22.

12 LIN Y P, CHENG B Y, GUEYSHIN S, et al. Combining a finite mixture distribution model with indicator Kriging to delineate and map the spatial patterns of soil heavy metal pollution in Chunghua County, central Taiwan[J]. Environmental Pollution, 2010, 158(1):235-244.

13 王政權. 地統計學及在生態學中的應用[M]. 北京: 科學出版社，1999: 65-132.

14 楊奇勇，楊勁松，余世鵬. 禹城市耕地土壤鹽分與有機質的指示克里格分析[J]. 生態學報，2011，31(8): 2196-2202.

YANG Q Y, YANG J S, YU S P. Evaluation on spatial distribution of soil salinity and soil organic matter by indicator Kriging in Yucheng City[J]. Acta Ecologica Sinica，2011，31(8): 2196-2202. (in Chinese)

15 MYERS J C. Geostatistical error management[M]. New York: Van Nostrand Reinhold, 1997.

16 WU C F, WU J P, LUO Y M, et al. Spatial estimation of soil total nitrogen using Cokriging with predicted soil organic matter content[J]. Soilence Society of America Journal, 2009, 73(73):1676-1681.

17 覃朝科，易鷂，劉靜靜，等. 廣西某鉛鋅礦區廢水匯集洼地土壤重金屬污染調查與評價[J]. 中國巖溶，2013, 32(3):318-324.

QIN C K, YI Y, LIU J J, et al. Investigation and evaluation on heavy metal pollution in the waste water collecting depression in a lead-zinc mine of karst area[J]. Carsologica Sinica, 2013, 32(3):318-324. (in Chinese)

18 鄧琴，吳迪，秦樊鑫，等. 巖溶鉛鋅礦區土壤重金屬污染特征[J]. 中國巖溶, 2017, 36(2):248-254.

DENG Q, WU D, QIN F X, et al. Pollution characteristics of heavy metals in soil of lead-zinc mining in karst areas[J]. Carsologica Sinica, 2017, 36(2):248-254. (in Chinese)

19 CHEN T, YAN B. Fixation and partitioning of heavy metals in slag after incineration of sewage sludge[J]. Waste Management, 2012, 32(5):957-964.

20 QUAN S X, YAN B, LEI Q, et al. Distribution of heavy metal pollution in sediments from an acid leaching site of e-waste[J]. Science of the Total Environment, 2014, 499(1): 349-355.

21 DUFFERA M, WHITE J G, WEISZ R. Spatial variability of Southeastern U.S. Coastal Plain soil physical properties: implications for site-specific management[J]. Geoderma, 2007, 137(3-4):327-339.

22 朱禮學. 土壤pH值在CaCO3在多目標地球化學調查中的研究意義[J]. 物探化探計算技術，2001, 21(2): 140-143.

ZHU L X. A study to the multi-aim geochemistry agrology survey[J]. Geophysical and Geochemical Exploration，2001, 21(2): 140-143. (in Chinese)