基于時空不確定性分析的北京市農田土壤重金屬鎘含量等級劃分

李曉嵐，高秉博，周艷兵*，潘瑜春，郜允兵，李 斌，胡茂桂

（1.北京市農業(yè)物聯網工程技術研究中心，北京100097；2.農業(yè)部農業(yè)信息技術重點實驗室，北京100097；3.北京農業(yè)信息技術研究中心，北京100097；4.國家農業(yè)信息化工程技術研究中心，北京100097；5.中國科學院地理科學與資源研究所資源與環(huán)境信息系統國家重點實驗室，北京100101）

隨著城鎮(zhèn)化和工業(yè)化的不斷發(fā)展，土壤重金屬污染問題日趨顯著，對人類健康和農產品安全構成嚴重威脅，亟需尋求有效的應對措施[1-2]。根據土壤環(huán)境質量將農田劃分為不同的等級，是進行差別化管理和治理的前提條件?！锻寥牢廴痉乐涡袆佑媱潯窂娬{，要實施農用地分類管理，保障農業(yè)生產環(huán)境安全[3]。土壤環(huán)境質量等級劃分一般通過將污染物含量與相關標準和研究規(guī)定的閾值對比來實現，如原環(huán)保局制定的《土壤環(huán)境質量標準》（GB 15618—1995）[4]和原環(huán)保部制定的《食用農產品產地環(huán)境質量評價標準》（HJ/T 332—2006）等[5]。目前，采樣依然是土壤重金屬含量調查最主要的技術手段[6-10]，基于稀疏樣本數據獲得目標區(qū)域農田的重金屬等級，需要基于樣本數據進行空間上的插值。

空間插值方法可以分為確定性方法和統計學方法兩類。確定性方法如反距離插值、樣條函數等，確定性方法只能給出未知點的估計值，但是不能給出估計值的不確定性。由于樣本數據不足及土壤重金屬含量機理不清，插值結果存在一定誤差，需要結合插值結果的不確定性對插值結果進行解讀。統計學方法能夠同時給出未知點的估計值和估計值的不確定性，適用于土壤重金屬含量插值。簡單克里金、普通克里金或趨勢克里金等方法常用于土壤重金屬含量插值[11]，并通過將插值結果與等級劃分閾值進行對比，進而劃分農田土壤重金屬含量等級。由于其插值結果存在一定的平滑效應，使得較小值和較大值向中間值靠攏，當閾值較大或者較小時，根據插值結果和閾值關系進行等級劃分容易產生錯誤[12-13]；另外，這幾種克里金插值需要基于估計誤差服從正態(tài)分布的前提假設，而實際土壤重金屬樣點數據一般具有較高的偏度，并不滿足該假設。土壤重金屬含量等級劃分具有不同于含量插值的特點，它只需準確估計未知點上的含量與等級劃分閾值之間的大小關系，而并不需要精確含量值。指示克里金是一種非參數地統計方法，它通過直接與閾值進行對比將含量數據轉化為0和1 的指示值，能夠直接推斷得到未知點重金屬含量小于閾值的概率，能夠減少含量值平滑后的等級劃分誤差，也不需要限定重金屬含量的分布類型，更適合用于農田土壤重金屬等級劃分[14-16]。

由于政府管理部門和研究機構對農田環(huán)境的長期重視，在實際工作中，很多區(qū)域積累了多期農田土壤重金屬采樣數據。而目前采用多期土壤重金屬樣點數據的研究主要集中于對多期土壤樣品重金屬含量的統計特征對比和分析[17-19]，或是采用相關時空地統計模型對重金屬含量分布進行預測[20]；針對土壤重金屬進行等級劃分的研究則主要基于某一期土壤重金屬含量數據計算地累積指數、潛在風險系數、內梅羅污染指數等來進行風險等級劃分[21-23]，利用多期數據對土壤重金屬的等級劃分目前研究較少。由于在時間維度上，不同時期的土壤重金屬含量數據間具有一定差異性，如果直接將時間尺度上所有時期數據整合成某一時期的數據直接進行插值，時間變異性將轉變成空間上的變異性，引入了時間變異的誤差，這將影響最終的等級劃分結果，因此需要充分利用多期時空采樣數據來提高農田土壤重金屬質量等級劃分的精度。在本研究中，采用基于時空指示克里金插值的方法從時空角度對土壤重金屬質量等級劃分展開研究，并分析土壤重金屬的等級分布特征及等級劃分的不確定性，以期為科學管理和決策提供一定數據基礎，促進土地資源有效利用，為政府的土壤環(huán)境等級劃分工作提供輔助支撐。

1 材料與方法

1.1 研究區(qū)

圖1 研究區(qū)樣點Figure 1 The sampling points of study area

北京市，位于華北平原北部，總面積16 410.54 km2，其中山地約占總面積的62%，平原約占總面積的38%。平均溫度為8～12 ℃，年平均降水量584 mm 左右。其區(qū)域圖如圖1 所示，地形西北高、東南低，母質類型多樣。北部丘陵地區(qū)成土母質以片麻巖、花崗巖、安山巖為主，西部以灰?guī)r和砂頁巖為主；平原地區(qū)成土母質質地疏松，以砂壤、輕壤和黃土性母質為主。受地形和母質類型的影響，遠郊區(qū)以生產小麥、玉米和水稻為主，而近郊區(qū)以蔬菜、果樹和農田防護林為主。隨著北京城郊農業(yè)和現代化農業(yè)發(fā)展，大量農藥、化肥和有機肥投入逐年增加。污灌面積也在急劇增加，目前北京地區(qū)污水灌溉面積近800 km2，占全市耕地面積的19%，占灌溉農田面積的24%，其中大多分布在通州、大興、朝陽、房山等北京城水流下游郊區(qū)縣。

1.2 采樣點

基于由北京市農林科學院相關團隊管理的農產品質量與農田環(huán)境基礎數據平臺，本研究從中分別提取了2005、2006、2007、2009、2011、2012、2013 年這7年的農田土壤重金屬采樣點數據，這些樣本均是基于相同的方法進行采集和測定。采集時使用GPS定位，基本覆蓋北京市農田區(qū)，其點位分布情況如圖1 所示。每個樣點單元均為邊長10 m 的正方形，并在每個單元內的土壤表層0～25 cm深度范圍內采集3～5個點。按四分法取1.0 kg 分析樣品作為代表該點的混合樣品，土壤樣品經王水水浴加熱消解，重金屬鎘的含量由石墨爐原子吸收法測定，在分析過程中引入重復樣品和質控樣品進行質量監(jiān)控。

1.3 時空指示克里金

時空指示克里金方法首先將原變量值轉換成為指示變量[24]，如公式（1）所示[25]：

式中：IZ（s，t；Zc）為基于閾值Zc的指示值；Zc為閾值；Z(s,t)為原變量值；s為某空間研究區(qū)域；t為某時間。

某時間t內具有全覆蓋數據的研究區(qū)域s，不超過閾值Zc的累積概率可以通過指示值進行計算：

式中：F（Zc）為不超過閾值的累計概率；Prob[Z（s，t）]≤Zc）為小于閾值的概率（或者頻率）；E[IZ（s，t；Zc）]為基于某時間段某空間范圍內指示值的期望值。研究區(qū)域基于閾值的累積概率可以轉化為指示值的期望值。

針對在未抽樣點處，隨機變量Z（s，t）低于閾值的累積概率同樣可以使用樣本點指示值的加權獲得，同時指示變量的無偏最優(yōu)估計結果即為Z（s，t）累積概率的無偏最優(yōu)估計結果。即：

式中：F[s，t；Zc|（n）]為隨機變量低于閾值的累計概率；λi為相應樣本點指示值的權重。

類似于普通克里金，以公式（3）中估計結果的無偏最優(yōu)作為限制條件，求解加權系數，計算未抽樣點上目標變量低于閾值的累積概率估計值及其估計誤差。

為了求解無偏最優(yōu)估計，基于平穩(wěn)假設，需要定義基于樣本指示值計算的半變異函數，基于樣本指示值計算的時空半變異函數[26]：

針對時空半變異函數相關學者們進行了大量研究[27-31]。如度量模型、乘積模型等。這些模型也被相關學者用于氣溫、水、輻射監(jiān)測、疾病等領域進行時空插值研究[32-37]。本研究是將時間差異性和空間差異性結合在一起進行研究，因此采用的是時空和度量模型（SumMetric Model）[38]，該模型能擬合出時空地理數據的時空變異結構，其時空半變異函數可以表示為：

式中：γ（hS，hT）為時空半變異函數；γ（hS）為空間尺度半變異函數；γ（hT）為時間尺度半變異函數；γST（hsT）為時間與空間相互關系的半變異函數。

1.4 基于時空指示克里金的等級劃分

在進行等級劃分中，根據一定的閾值將研究區(qū)域的研究對象劃分成相應的等級?；跁r空指示克里金的等級確定的基本步驟如下：

（1）指示變換：確定k 個等級閾值Zc1，Zc2，…，Zck，分別將目標變量分為C0，C1，C2，…，Ck類別，其中，C0=（0，Zc1]，C1=（Zc1，Zc2]，…，Ck=（Zc1，∞]。將每個閾值按照公式（1）對樣本值進行指示變換，得到k 個指示數據集。

（2）計算指示半變異函數：分別計算每個指示數據集所對應的空間半變異函數，或用中位數作為閾值的半變異函數替代每個等級閾值的半變異函數。

（3）指示克里金插值計算：依次針對閾值Zc1，Zc2，…，Zck，對每一個空間單元使用指示克里金進行空間插值，獲得每個空間單元小于第k 個閾值的概率P0，P1，P2，…，Pk及其估計誤差的標準差sp1，sp2，…，spk。

（4）空間單元等級確定：為每個閾值的概率設置判定的概率閾值Pc1，Pc2，…，Pck，結合上一步計算得到概率估計及其誤差對每個空間單元進行類別判定。針對每個空間單元，判定方法如下：

①針對第i 值Zci，計算估計概率的置信區(qū)間，[Pispi，Pi+spi]。

②判斷空間單元等級：如果i=1，則按照公式（7）對空間單元進行等級劃定[16]：

式中：ci為某等級；Ci為第i 級；Ci-1為第i-1 級；C-1表示未分等級；Pci為第ci 個概率閾值；Pi為第i 個概率閾值；spi為第pi個估計誤差的標準差。

如果i≠1，按照公式（8）進行等級劃定[16]

③判斷是否已經對所有閾值進行計算，即i 是否達到k，若達到，則結束，否則轉①繼續(xù)執(zhí)行。

（5）劃分等級邊界：合并同等級空間單元，形成等級邊界。

（6）等級劃分不確定性分析：通常越是在概率估計值接近概率閾值且估計誤差較大處，越容易劃錯，等級劃分的不確定性則越強，等級劃分不確定性可以從概率估計值及概率估計誤差兩方面來確定，并采用錯劃指數來反映其大小?；诟怕书撝颠M行等級劃分，當概率閾值估計誤差較小可忽略不計時，概率估計值越接近概率閾值，等級錯劃概率越大；而概率估計值具有一定的估計誤差，特別在一些變異較大且點位稀疏的位置，估計誤差通常較大，由此造成的等級劃分不確定性也較強。公式（9）計算得到基于概率估計值的錯劃指數，公式（10）計算得到基于估計誤差的錯劃指數，基于公式（11）獲得最終的綜合錯劃指數，其大小反映了等級劃分不確定性的強弱。

式中：Ic（x；Zc）為概率估計值的錯劃指數；Ie（x；Zc）為概率估計值誤差的錯劃指數；Is（x；Zc）為綜合錯劃指數；Pc為概率閾值(x；Zc)為概率估計值；δ（x；Zc）為概率估計值誤差。

1.5 等級劃分概率閾值的自適應確定方法

由上述劃分步驟可知，概率閾值Pc對最終等級劃分確定結果至關重要。時空指示克里金計算結果為小于閾值的概率，本文中一律采用小于閾值的概率。

在本研究中，通過設定不同的概率閾值來提取基于概率閾值判斷的等級劃分結果，將劃分結果與真實值的等級進行比對，提取不同概率閾值對應的錯劃指數，最終選擇較小錯劃指數所對應的概率閾值作為最終進行等級劃分的概率閾值，從而實現根據研究對象數據特征自適應確定相應的概率閾值。其具體步驟如下：

（1）首先將概率閾值分別設置為0.1、0.2、0.3、…、0.9；在某個概率閾值下，依次保留一個已知樣點，采用其他樣點基于時空指示克里金對該樣點進行等級劃分。

（2）將基于時空指示克里金的等級劃分結果與真實值的等級進行對比，統計由Juang 等[39]提出的第一類錯誤T1、第二類錯誤T2、綜合錯誤E。其中第一類錯誤T1為真實數據未超過閾值而統計推斷估計結果將其判斷為超過閾值的錯誤，第二類錯誤T2為真實數據超過閾值而統計推斷估計結果將其判斷為未超過閾值的錯誤。綜合錯誤E即為這兩類錯誤的總和。

式中：Zc表示等級閾值；Z（x）表示x 位置上的真實值；(x；Zc)表示x 位置上基于閾值Zc的估計值；n 表示研究區(qū)域總數據量；I（xi；zc）表示基于閾值的指示值；I*（xo；zc）表示克里金插值概率估計值；Pc為概率閾值。

（3）最后繪制等級劃分誤差圖，根據研究對象數據特征和研究需求自適應確定合適的概率閾值。如圖2 示例，在等級劃分錯誤比例圖中，隨著概率閾值的增大，第一類錯誤增加，第二類錯誤則會減少。如果需要第一類錯誤最小，概率閾值選擇0.1；如果需要第二類錯誤最小，概率閾值選擇0.9；若需要第一類錯誤和第二類錯誤最相近，則概率閾值選擇0.5；若需要錯劃比例最小，則概率閾值選擇0.4。

圖2 基于概率閾值的等級劃分錯誤比例示例圖Figure 2 The grade classification error ratio diagram based on probability threshold

根據上述時空指示克里金的等級劃分方法以及等級劃分概率閾值的自適應確定方法的實現步驟，研究采用Java及R語言，實現了基于時空不確定性分析的農田土壤環(huán)境質量等級劃分軟件。該軟件依托于gstat、rgeos、rgdal等R包，并調用vgmST函數和krigeST函數來構建相關的時空半變異函數和時空指示克里金插值計算，對土壤重金屬含量進行等級劃分和概率閾值的自適應確定。

1.6 性能對比分析

為了檢驗基于時空指示克里金的等級劃定精度，本研究選擇了基于時空普通克里金的等級劃分結果與基于時空指示克里金的等級劃分結果進行比對。具體的對比流程圖如圖3所示，其具體步驟如下：

（1）準備數據。基于研究中所用的北京市土壤重金屬含量數據，分別準備一套重金屬含量原始值數據和一套基于北京市背景值進行指示變換后的重金屬含量指示值數據。

（2）通過隨機抽樣，確定用于交叉檢驗的待預測點位。為確保待預測點位的代表性，以2013 年數據為基礎，隨機抽樣選取一半樣點作為最終的待預測點位。這套點位同時運用于下一步驟中兩種方法的插值。

（3）采用兩種方法分別進行插值。使用2013 年以前的所有年份采樣點以及2013 年未被選中而剩下的樣點作為插值樣點，分別基于各自的數據進行時空指示克里金插值和時空普通克里金插值。

圖3 等級劃分性能對比分析流程圖Figure 3 Flow chart of the grade classification comparison analysis

（4）統計預測點位的等級錯劃比例。以上一節(jié)確定的概率閾值為最終的概率閾值，基于時空指示克里金插值結果，按照公式（12）、公式（13）、公式（14）計算預測點位的等級錯劃比例；以北京市土壤重金屬的某值為最終的閾值，基于時空普通克里金插值結果，按照公式（12）、公式（13）、公式（14）計算預測點位的等級錯劃比例。

（5）基于時空指示克里金的等級劃定的錯劃比例結果和基于時空普通克里金的等級劃定的錯劃比例，分析基于時空指示克里金等級劃定的性能。

2 結果與討論

2.1 土壤重金屬元素的時空變異分析

土壤重金屬鎘含量的描述性統計結果如表1 所示，2005、2006、2007、2009、2011、2012 年和2013 年鎘含量是不斷變化的。其中鎘含量的均值在2011 年達到最大為0.232 mg?kg-1，在2007年最小0.156 mg?kg-1；變異系數最小值為2009 年的53.058%，最大值為2005年的176.476%；偏度最大值為2012年的11.760，最小值為2009 年的2.270；峰度最大值為2012 年的190.700，最小值為2009年的5.903。

以陳同斌等[40]研究獲得的北京市土壤中重金屬鎘背景值0.119 mg?kg-1為閾值，將所有年份的樣點數據按照公式（1）進行指示變化，小于背景值的賦值為1，大于背景值的賦值為0。對時空指示克里金的時空半變異函數進行了擬合，如圖4 所示。其中指示克里金的輸入值為鎘含量值的指示值。圖4a及圖4b分別為時空指示克里金最終擬合的時空半變異函數的散點圖和三維圖。

2.2 土壤重金屬等級劃分概率閾值確定

基于擬合的時空半變異函數進行時空指示克里金插值，得到了最終的土壤重金屬鎘的概率值。根據等級劃分中概率閾值的自適應確定方法，將預測點位的插值結果與其對應的真實值進行比較，分別統計對應的第一類錯誤、第二類錯誤和綜合錯劃比例，結果如圖5 所示?；诘谝活愬e誤與第二類錯誤的和在概率閾值為0.4的情況下達到最低，且后續(xù)趨近平穩(wěn)，因此本研究中確定最終等級劃分的概率閾值為0.4。

表1 土壤重金屬鎘含量的時空統計特征Table 1 Spatial-temporal statistics characteristics of soil metal cadmium content

圖4 IK的時空半變異函數散點圖及三維圖Figure 4 IK spatial-temporal semi-variogram plot and 3D graph

圖5 概率閾值錯誤比例統計Figure 5 Probability threshold error proportion statistic

2.3 等級劃分精度對比分析

按照1.6 節(jié)中的性能對比分析步驟，分別計算基于時空指示克里金插值結果和時空普通克里金插值結果的等級劃分錯劃比例，其中時空指示克里金插值結果以2.2節(jié)確定的概率閾值0.4為最終的概率閾值，而時空普通克里金插值結果以北京市農田土壤重金屬鎘的平均土壤背景值0.119 mg?kg-1為最終的閾值。兩種插值方法最終統計的錯誤比例如圖6 所示，其中基于時空指示克里金等級劃分的綜合錯誤數的錯誤比例為14.41%，而基于時空普通克里金等級劃分的綜合錯誤數的錯誤比例為18.92%，交叉檢驗結果表明基于時空指示克里金插值的等級劃分的錯誤數小于基于時空普通指示克里金插值的等級劃分所產生的錯誤數，由此可看出，在對北京市農田土壤重金屬鎘含量進行等級劃分時，相對于時空普通克里金，采用基于時空指示克里金的等級劃分方法能取得較高的等級劃分精度。

2.4 基于時空指示克里金的北京土壤重金屬等級劃分結果和不確定性分析

圖6 不同方法對應的錯誤個數Figure 6 The classification error proportion of IK and OK

基于時空指示克里金的等級劃分確定方法，最終得到了2013 年北京市農田土壤重金屬鎘含量等級分布圖以及鎘含量等級的錯劃指數分布圖，分別如圖7和圖8所示。圖7反映了研究區(qū)土壤重金屬鎘含量等級的分布情況，鎘含量等級分為高于背景值和低于背景值兩個級別，紅色區(qū)域為高于背景值等級，綠色區(qū)域為低于背景值區(qū)域。圖8 反映了鎘含量等級的錯劃指數分布情況，錯劃指數越大，即錯劃的可能性較大，說明鎘含量的等級劃分的不確定性越強；錯劃指數越小即錯劃的可能性較小，則說明鎘含量等級劃分的不確定性越低。

圖7 鎘含量的等級分布圖Figure 7 The grade distribution map of cadmium content

圖8 鎘含量等級的錯劃指數分布圖Figure 8 The error index distribution map of cadmium content grade

結合圖7鎘含量等級分布圖及圖8錯劃指數分布圖分析，在不確定性較低區(qū)域，大興南部地區(qū)、昌平部分地區(qū)、平谷中部地區(qū)、房山南部地區(qū)、朝陽區(qū)等地的鎘含量等級極大可能為高于背景值；延慶大部分地區(qū)、懷柔部分地區(qū)等地的鎘含量等級極大可能為低于背景值；其他區(qū)域還需進一步調查確認。其中鎘含量等級為大于背景值的區(qū)域與蔣紅群等[7]預測的鎘含量高風險區(qū)及中警風險區(qū)一致，且與實際情況吻合。根據以往研究[41-43]可知，農田區(qū)域鎘含量主要是來源于含鎘熱穩(wěn)定劑的地膜降解、污水灌溉、采礦活動、大氣沉降及其他人類活動。距離城鎮(zhèn)中心較近的農田區(qū)域，受人類活動影響相對更加強烈，重金屬鎘更容易累積，因此在大興南部地區(qū)、昌平部分地區(qū)、平谷中部地區(qū)、房山南部地區(qū)、朝陽區(qū)等距離城鎮(zhèn)中心較近的農田區(qū)域，其鎘含量等級極大可能為高于背景值；而距離城鎮(zhèn)中心較遠的農田區(qū)域本身成土母質含量低，受人類活動影響較小，重金屬鎘累積效應較弱，因此在延慶大部分地區(qū)、懷柔部分地區(qū)等遠郊區(qū)的農田區(qū)域，其鎘含量等級極大可能低于背景值。

部分農田區(qū)域鎘含量等級不確定性較強，其主要原因包括兩方面：一方面是受樣點分布影響，在點位稀疏區(qū)域其估計誤差較大，尤其是在變異較大位置；另一方面是受概率閾值影響，基于時空指示克里金插值獲取得到的概率估計值越接近概率閾值時，其等級錯劃概率就越大。這兩個因素在鎘含量等級不確定性的貢獻度方面可依據上述相應公式計算來進一步分析。針對不確定性較強的農田區(qū)域，可根據生產生活或監(jiān)測管理需求進行補充調查，為加密采樣提供了指導；或者結合相關輔助數據進一步確認，但是由于時空指示克里金無法利用輔助數據，后續(xù)研究將嘗試引入協變量對該方法進行改進。

本文基于鎘含量背景值利用基于時空指示克里金的等級劃分方法對北京市農田土壤重金屬鎘含量進行了等級劃分。根據原農業(yè)部印發(fā)的“土十條”實施意見，后續(xù)可將該方法進一步廣泛應用于耕地土壤環(huán)境質量類別劃分工作中，為農田土壤環(huán)境質量科學管理和決策提供支撐。

3 結論

（1）研究采用的基于時空指示克里金的等級劃分方法，能夠利用多期時空采樣數據對土壤環(huán)境質量進行等級劃分，可以自適應確定等級劃分的概率閾值，并能估計等級劃分的不確定性程度。該方法可為耕地土壤環(huán)境質量類別劃分工作提供輔助支撐。

（2）北京市2013 年農田土壤重金屬鎘含量的等級分布顯示：在昌平大部分地區(qū)、平谷中部地區(qū)、大興南部地區(qū)、房山南部等距離城鎮(zhèn)中心較近區(qū)域，其農用地鎘含量等級極大可能為高于背景值；而在延慶西部地區(qū)、懷柔北部地區(qū)等遠離城鎮(zhèn)中心區(qū)域，其農用地的鎘含量等級極大可能為低于背景值。

（3）北京市2013 年農田土壤重金屬鎘含量等級錯劃指數分布反映了北京市農田土壤重金屬鎘含量等級劃分的不確定性程度。其受樣點分布及概率閾值確定的影響。