基于GIS的凱江流域農田重金屬污染評價研究

侯沁言，張世熔*，馬小杰，潘小梅，王貴胤，鐘欽梅

（1.四川農業大學環境學院，成都 611130；2.四川省土壤環境保護重點實驗室，成都 611130）

近年來，由于人類活動強度的提高，農田土壤重金屬污染已成為全球關注的環境問題之一[1-2]。在污染土壤的重金屬元素中，Cd、Pb、Hg、As、Cr、Zn、Cu、Ni等8種元素的研究報道較多[3-6]。值得關注的是，即使在較低含量水平條件下，土壤Cd、Pb、Hg和As也不僅會影響農產品的產量和品質，還會通過食物鏈傳遞，對人體健康造成潛在威脅[5]。作為我國糧食、油料等作物主要生產區之一的四川盆地丘陵地區，其農田土壤重金屬的研究已有少量報道[7-9]。例如，在尚二萍等[7]的研究結果中，四川盆地耕地點位超標率為43.55%，而嚴明書等[8]和劉潔等[9]對重慶都市經濟圈土壤及沱江中游土壤質量的研究表明，其點位重金屬超標率分別為13.59%和0.6%。盡管樣點超標率可以描述土壤重金屬的污染狀況，但不同研究結果較大的空間異質性也增大了準確評估區域尺度重金屬污染分布特征的難度。因此，進一步開展該區域土壤重金屬空間分布及污染評價研究，為區域土壤環境質量管理提供參考依據是十分必要的。

關于區域土壤重金屬污染評價方法，傳統上多采用基于樣點的單因子污染指數法（Pi）、內梅羅綜合指數法（PN）、地累積指數法（Igeo）和潛在生態風險指數法（RI）等方法進行評價[10-12]。然而，重金屬在土壤中的移動性較差，基于樣點的指數法難于提供其污染空間分布的真實信息[13]。特別是在丘陵地區，重金屬在同一農田地塊內及不同地塊間的遷移易受地形的限制。隨著丘陵地區土壤重金屬污染防治進入規劃與工程實施階段，提供更準確的污染空間分布信息成為必然要求。因此，近年來一些研究者基于GIS、GPS和地統計學等空間分析方法研究了區域土壤重金屬的空間分布特征[9，14-16]。同時，也有研究者將這些空間分析法與指數法結合對其進行評價[17]。然而，迄今基于點位的污染超標率與基于空間分布的超標率之間的關系研究仍少見報道。此外，在相同的研究尺度（丘陵區）、采樣方案及軟件平臺下，系統性地對比不同空間插值方法對農田土壤重金屬超標率的影響仍鮮見報道。

凱江流域為四川盆地丘陵區涪江流域的一級支流，是四川盆地典型的糧油及蔬菜生產基地之一。流域內種植業和畜禽養殖業發達，人口密度高達500～700人·km-2。近年來，隨著農業集約化程度的提高和鄉鎮企業的發展，其農田土壤重金屬的含量水平及污染狀況尚不清楚。因此，本文根據我國《土壤環境質量農用地土壤污染風險管控標準（試行）》（GB 15618—2018）[18]，采用單因子指數法，在GIS和GPS軟件平臺上，系統對比研究反比距離權重法（Inverse distance weighted，IDW）、徑向基函數法（Radial basis function，RBF）、普通克里格法（Ordinary kriging，OK）、析取克里格法（Disjunctive kriging，DK）和概率克里格法（Probability kriging，PK）5種空間插值方法對該區域農田土壤重金屬污染評價結果的影響，旨在為區域土壤污染治理提供科學參考。

1 材料與方法

1.1 研究區域概況

研究區域地處東經 104°40′～105°6′和北緯 30°56′～31°4′，面積約 471 km2，屬亞熱帶濕潤性季風性氣候，年平均氣溫15℃，年均降水量800 mm，具有氣候溫和、四季分明、雨水充沛等特點。研究河段處于四川盆地丘陵低山地區，其土地利用方式主要為耕地，以種植水稻、小麥、棉花和油菜等為主。成土母質以侏羅系蓬萊組紫色砂泥巖、白堊系城墻巖群黃色砂巖、厚砂薄頁巖風化坡殘積物為主，凱江主流及支流沿岸有部分地區為第四系新沖積物。

1.2 土壤樣品采集與分析方法

根據該區地形地貌、成土母質、土壤類型、土地利用等條件，按隨機網格法和分層抽樣法進行樣點的布設，共187個采樣點及20個檢查點（圖1）。在每個樣點按五點采樣法采集同一地塊0～20 cm表層土樣，混合均勻后用四分法保留土樣約1 kg；樣點位置采用GPS進行定位，并調查記錄地形、母質、土壤類型等自然環境條件，以及作物、灌溉、施肥等社會經濟信息。

圖1 樣點空間分布圖Figure 1 Distribution of soil sampling points

土壤樣品于實驗室自然風干后，去掉礫石、動植物殘骸等，進行研磨，再分別過20目、100目篩，待測。土壤pH值在水土比為2.5∶1條件下采用復合電極法進行測定；土壤Cd、Pb的總量用HNO3-HF-HClO4體系（V∶V∶V=1∶2∶2）消解，消解液冷卻后，加入5 mL 2 mol·L-1HCl溶解，并轉移至25 mL容量瓶，定容、過濾。待測液采用石墨爐原子吸收分光光度法測定；土壤Hg和As總量采用HCl與HNO3混合液體系（V∶V=1∶1）消解，原子熒光法測定[19]。在分析過程中所用試劑均為分析純，所有土壤樣品中每個重金屬全量測定時均設置3個重復，并設置空白樣，檢測過程中采用國家標準土壤樣品（GSS-3）進行質量控制，Cd、Hg、As和Pb的平均回收率分別為106.8%、92.5%、96.8%和102.3%。

1.3 單因子污染指數法

單因子污染指數法作為最簡單的評價方法之一，通過實際含量和評價標準的比值直接表現出重金屬的污染級別及危害程度。其計算公式[7]為：

式中：Pi為第i種土壤重金屬的環境質量指數；Ci為土樣中第i種重金屬的實測含量，mg·kg-1；Si為第i種重金屬的農用地土壤污染風險篩選值（試行）（GB 15618—2018）[18]，mg·kg-1。本文評價分為兩個等級：當 Pi＜1時，土樣重金屬含量低于篩選值，樣點基本未受到污染；當Pi≥1時，土壤樣點受到污染，指數值越高，污染的程度越高。

1.4空間插值方法

本研究基于樣點GPS坐標和ArcGIS 10.2地理信息系統軟件平臺，采用地統計學的普通克里格法進行空間插值獲得區域土壤4種重金屬含量的空間分布圖。同時，為系統對比研究不同插值方法對空間分布的影響，在相同的樣點GPS坐標和ArcGIS 10.2平臺上，分別采用反比距離權重法和徑向基函數法兩種確定性方法，以及普通克里格法和析取克里格法兩種地統計學方法對單因子指數Pi進行空間插值[20]，以研究不同空間插值方法對該區域農田土壤重金屬污染評價結果的影響。其中，基于地統計學的普通克里格法和析取克里格法均要求數據為正態分布。此外，本研究還采用概率克里格法[20]對單因子指數Pi進行了不同污染水平的風險分析。最后，根據以上5種空間插值方法獲得的信息，綜合討論該區域農田土壤重金屬污染超標狀況。

1.5 數據處理方法

本研究所有數據統計均在SPSS 19.0平臺上進行；樣點分布圖和各類空間插值圖利用ArcGIS 10.2軟件繪制。

2 結果與討論

2.1 土壤中重金屬含量統計

研究區土壤樣點pH范圍為4.90～8.25。土壤中Cd、Hg、As和Pb 4種重金屬均呈正態分布（表1），平均含量分別為 0.38、0.10、7.77 mg·kg-1和 22.54 mg·kg-1。其中，土壤Cd、Hg、As和Pb的平均含量分別為當地土壤環境背景值的1.58、1.43、1.27倍和1.11倍，表明該區域存在不同程度的重金屬累積情況。根據GB 15618—2018[18]，4種重金屬元素的平均值均低于其風險篩選值，但有部分土壤樣點的Cd和As含量超過其污染風險篩選值。變異系數（CV）相較于標準差更能反映出數據的離散性[21]，土壤Cd、Hg、As、Pb的變異系數分別為50%、40%、29%和13%，均處于中度變異程度（10%＜CV＜100%）[13，22]，表明這幾種元素數據的離散性較小，分布相對較為均勻。

表1 土壤重金屬含量統計分析Table 1 Statistical analysis of soil heavy metals

2.2 土壤重金屬含量空間特征分析

由于該區域土壤Cd、Hg、As和Pb含量均符合正態分布，因此本研究在ArcGIS 10.2平臺上選擇地統計學模塊中的普通克里格法進行空間插值，獲得它們含量的空間分布圖（圖2）。研究區土壤Cd（圖2a）含量在西北角和凱河一帶有兩個高值中心（＞0.8 mg·kg-1），同時呈現向四周遞減的趨勢。研究區西北角靠近縣城，工業活動頻繁；同時，在凱河附近曾有一養豬場，其養殖廢棄物可能是造成污染的主要原因。土壤Hg的空間分布（圖2b）呈現以南華北部為高值區（＞0.17 mg·kg-1），向東、向北逐漸減少到 0.06～0.10 mg·kg-1的低值區，并在東南部最低（＜0.06 mg·kg-1）。由于研究區土壤Hg平均含量較低，因此該分布可能與土壤母質因素有較大關系。土壤As的空間分布（圖2c）呈現以南華和萬安的兩個高值中心（＞11 mg·kg-1），向四周遞減的趨勢。研究區內土壤Pb（圖2d）的高值中心（＞25 mg·kg-1）位于南華、凱河處，并向周圍呈現遞減趨勢，其東北及西北部屬于低值區（＜20 mg·kg-1）。有研究表明，土壤中As和Pb主要源于交通污染和大氣沉降[23-24]。研究區土壤As和Pb的高值區均處于交通較好的區域，同時在南華附近存在一個垃圾焚燒廠，這也可能是造成高值的原因。

圖2 凱江流域土壤重金屬空間分布圖Figure 2 Spatial distribution of soil heavy metal in the Kaijiang Valley

2.3 土壤重金屬污染評價

2.3.1 基于樣點的單因子污染指數評價

研究區內土壤Cd、Hg、As和Pb的單項污染指數的平均值分別為0.487、0.111、0.381和0.102（表2），其中僅土壤Cd和As存在Pi≥1的超過風險篩選值的樣點，其占比分別為7.5%和0.5%。因此，Cd是研究區農田土壤中主要的重金屬污染物。

2.3.2 基于空間的單因子污染指數評價

為進一步探究研究區重金屬污染指數的空間分布特征，分別采用反比距離權重法、徑向基函數法、普通克里格法和析取克里格法對4種重金屬的Pi進行空間插值，并分級統計其面積占總面積的百分比（表3）。

總體上4種插值方法的平均誤差（ME）、標準化平均誤差（MSE）和20個校驗點的檢驗誤差（MEP）均較小（表3），表明它們用于插值是可行的，但普通克里格法的誤差相對較小。從土壤重金屬種類分析，Hg和Pb的4種空間插值分析方法的統計結果完全相同，全區域均為Pi＜0.4等級的面積分布。因此，這兩種土壤重金屬在研究區內無污染風險。然而，不同空間插值方法所得的土壤Cd和As的面積統計結果卻呈現一定差異（表3）。

土壤 As僅存在 Pi＜0.4、0.4≤Pi＜0.6 和 0.6≤Pi＜0.8這3個污染指數等級的空間分布。因此，盡管這4種空間插值方法所得的分布等級相同，但在前3個等級各自分布的面積卻呈現一定差異。例如，在0.6≤Pi＜0.8等級的分布中，占區域面積百分數由小到大依次為析取克里格法面積（僅0.07%）＜徑向基函數法（0.10%）＜普通克里格法（0.67%）＜反比距離權重法（0.89%），與基于樣點污染指數的評價結果（2.7%）不同。此外，4種空間插值方法皆因僅有1個超標樣點而無法形成超標等級的面積分布。

表2 基于樣點的單因子指數評價Table 2 Assessment of single factor indexes based on sampling points

表3 基于不同插值方法的單因子指數分級面積比Table 3 Area ratios of single factor index based on the different interpolation methods

土壤Cd污染指數的空間分布則因插值方法不同呈現明顯差異。但無論哪一種空間插值方法，其超標區域占研究區總面積的百分比（0.14%～1.06%）均明顯低于基于樣點統計的百分比（7.5%）。

為進一步觀察4種空間插值方法對其污染指數，特別是超標指數空間分布的影響，本文選取在土壤質量環境評價研究中采用最頻繁的反比距離權重法[25]和普通克里格法[9]作圖，進行輔助分析（圖3和圖4）。對比圖3、圖4和表3可知，反比距離權重法和徑向基函數法中土壤Cd的Pi≥1超標等級面積（1.06%和0.67%）明顯高于普通克里格法（0.42%）和析取克里格法（0.14%），這是由于算法公式上的差異所致[26-27]。前兩種方法屬確定性插值方法，均突出高值點對其鄰近空間區域影響，其權重更大，所以插值土壤Cd的超標面積偏大（圖3a）；反之，基于地統計學的普通克里格法和析取克里格法并不強調高值點對其鄰近空間區域的影響權重，且不同含量樣點取等權重，因而插值結果具有一定的平滑效應[20，28]。特別是析取克里格法，其平滑效應明顯，導致其超標面積僅為反比距離權重法和普通克里格法的13.2%和33.3%。

圖3 基于反比距離權重法的單因子指數空間分布Figure 3 Spatial distribution of single factor index for soil heavy metals by inverse distance weighting method

2.4 單因子指數的風險概率評價

與前面4種研究單因子污染指數值空間分布特征的插值方法不同，概率克里格法[20]能基于污染指數的某一閾值水平來觀察研究區超過該閾值水平概率的空間分布特征[29]。毛竹等[30]利用該方法分析了漢源富泉鉛鋅礦區土壤重金屬的概率空間分布，預測出研究區各區域重金屬的超標概率。本研究基于土壤Cd、Hg、As、Pb的含量水平及單因子污染指數統計結果（表2），通過預測概率分析，將Cd和As的閾值設定為1.0，而Hg和Pb的最大閾值分別為0.3和0.2，以超過這些閾值為條件，用概率克里格法依次進行空間插值獲得它們的風險概率分布圖（圖5）。其中，Hg在Pi＞0.3及Pb在Pi＞0.2條件下均無高風險概率，因此整個研究區農田土壤受Hg和Pb污染的概率很低，該高值區域評價結果與前面基于樣點或空間的單因子指數評價結果均一致。

土壤Cd與As污染超標（Pi＞1）的風險概率區域占整個研究區面積的比例則與前面幾種評價不完全一致。其中，在萬安附近，有一小塊農田土壤As的超標概率達到0.8～1.0的水平（圖5c，占研究區面積的1.05%），即該區域土壤有80%的概率受到As污染。其評價結果與基于樣點的單因子污染指數評價相似，而未像基于空間的前4種插值方法一樣漏掉有1個樣點土壤As超標的信息。同時，土壤Cd存在兩個概率為1.0的高值區域，其空間分布與Pi值的超標區域空間分布特征相似，但面積約為區域面積的0.16%，更接近基于空間的析取克里格法的評價結果。

3 結論

（1）研究區土壤Cd、Hg、As和Pb均呈正態分布，平均含量依次是當地土壤環境背景值的1.58、1.43、1.27倍和1.11倍，存在不同程度的重金屬累積情況，其中有部分土壤點位Cd和As含量超過GB 15618—2018中的風險篩選值。

（2）采用普通克里格法獲得區域土壤中4種重金屬含量的空間分布圖。研究區內土壤Cd、As和Pb的含量總體上呈現以兩個高值中心向四周遞減的趨勢；土壤Hg則總體呈現西部高、東部低的趨勢。

（3）基于樣點的單因子指數評價得出土壤污染程度為 Cd＞As＞Hg＞Pb。基于空間的面積統計結果表明，研究區僅土壤Cd存在超標面積。4種方法中，普通克里格插值方法的誤差最小，更接近于真實情況。

（4）概率克里格空間插值表明研究區僅有土壤Cd和As存在高污染風險概率區域，其空間分布特征與Pi的高值區域空間分布特征基本一致。

圖4 基于普通克里格法的土壤重金屬單因子指數空間分布Figure 4 Spatial distribution of single factor index for soil heavy metals by ordinary kriging method

圖5 基于概率克里格法的土壤重金屬單因子指數概率空間分布Figure 5 Probability distribution of single factor index for soil heavy metals by probability