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

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

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

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

由于政府管理部門和研究機構(gòu)對農(nóng)田環(huán)境的長期重視，在實際工作中，很多區(qū)域積累了多期農(nóng)田土壤重金屬采樣數(shù)據(jù)。而目前采用多期土壤重金屬樣點數(shù)據(jù)的研究主要集中于對多期土壤樣品重金屬含量的統(tǒng)計特征對比和分析[17-19]，或是采用相關(guān)時空地統(tǒng)計模型對重金屬含量分布進行預(yù)測[20]；針對土壤重金屬進行等級劃分的研究則主要基于某一期土壤重金屬含量數(shù)據(jù)計算地累積指數(shù)、潛在風(fēng)險系數(shù)、內(nèi)梅羅污染指數(shù)等來進行風(fēng)險等級劃分[21-23]，利用多期數(shù)據(jù)對土壤重金屬的等級劃分目前研究較少。由于在時間維度上，不同時期的土壤重金屬含量數(shù)據(jù)間具有一定差異性，如果直接將時間尺度上所有時期數(shù)據(jù)整合成某一時期的數(shù)據(jù)直接進行插值，時間變異性將轉(zhuǎn)變成空間上的變異性，引入了時間變異的誤差，這將影響最終的等級劃分結(jié)果，因此需要充分利用多期時空采樣數(shù)據(jù)來提高農(nóng)田土壤重金屬質(zhì)量等級劃分的精度。在本研究中，采用基于時空指示克里金插值的方法從時空角度對土壤重金屬質(zhì)量等級劃分展開研究，并分析土壤重金屬的等級分布特征及等級劃分的不確定性，以期為科學(xué)管理和決策提供一定數(shù)據(jù)基礎(chǔ)，促進土地資源有效利用，為政府的土壤環(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 所示，地形西北高、東南低，母質(zhì)類型多樣。北部丘陵地區(qū)成土母質(zhì)以片麻巖、花崗巖、安山巖為主，西部以灰?guī)r和砂頁巖為主；平原地區(qū)成土母質(zhì)質(zhì)地疏松，以砂壤、輕壤和黃土性母質(zhì)為主。受地形和母質(zhì)類型的影響，遠(yuǎn)郊區(qū)以生產(chǎn)小麥、玉米和水稻為主，而近郊區(qū)以蔬菜、果樹和農(nóng)田防護林為主。隨著北京城郊農(nóng)業(yè)和現(xiàn)代化農(nóng)業(yè)發(fā)展，大量農(nóng)藥、化肥和有機肥投入逐年增加。污灌面積也在急劇增加，目前北京地區(qū)污水灌溉面積近800 km2，占全市耕地面積的19%，占灌溉農(nóng)田面積的24%，其中大多分布在通州、大興、朝陽、房山等北京城水流下游郊區(qū)縣。

1.2 采樣點

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

1.3 時空指示克里金

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

①針對第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個估計誤差的標(biāo)準(zhǔn)差。

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

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

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

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

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

1.5 等級劃分概率閾值的自適應(yīng)確定方法

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

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

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

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

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

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

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

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

1.6 性能對比分析

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

（1）準(zhǔn)備數(shù)據(jù)?；谘芯恐兴玫谋本┦型寥乐亟饘俸繑?shù)據(jù)，分別準(zhǔn)備一套重金屬含量原始值數(shù)據(jù)和一套基于北京市背景值進行指示變換后的重金屬含量指示值數(shù)據(jù)。

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

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

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

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

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

2 結(jié)果與討論

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

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

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

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

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

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

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

圖5 概率閾值錯誤比例統(tǒng)計Figure 5 Probability threshold error proportion statistic

2.3 等級劃分精度對比分析

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

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

圖6 不同方法對應(yīng)的錯誤個數(shù)Figure 6 The classification error proportion of IK and OK

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

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

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

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

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

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

3 結(jié)論

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

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

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