賈 琳，夏天翔*，張麗娜，賈曉洋，張 丹

基于不同空間插值的污染場地土壤銻修復(fù)范圍預(yù)測研究①

賈 琳1,2,3，夏天翔1,2,3*，張麗娜1,2,3，賈曉洋1,2,3，張 丹1,2,3

(1 北京市生態(tài)環(huán)境保護(hù)科學(xué)研究院，北京 100037；2 國家城市環(huán)境污染控制工程技術(shù)研究中心，北京 100037；3 污染場地風(fēng)險模擬與修復(fù)北京市重點實驗室，北京 100037)

以污染場地土壤銻超標(biāo)地層為研究對象，結(jié)合土壤每層銻含量的分布特點，分析和比較了克里格插值模型(OK、UK)、反距離插值(IDW)、泰森多邊形(TIN)和分區(qū)預(yù)測模型(OK+TIN)對不同深度土壤銻修復(fù)范圍及預(yù)測精度的差異。結(jié)果表明：各方法確定的0 ～ 1 m地層修復(fù)面積大小順序為IDW>OK=UK> TIN >OK+TIN，修復(fù)面積依次為74 154、57 427、37 338、32 707 m2；各方法確定的1 ～ 3.8 m地層修復(fù)面積大小順序為IDW>UK>TIN>OK+TIN>OK，修復(fù)面積依次為13 089、11 030、9 660、2 183、1 231 m2，克里格插值帶來的平滑作用導(dǎo)致超標(biāo)范圍的遺漏和未超標(biāo)區(qū)域納入修復(fù)范圍。在受成本和污染物空間分布制約的情況下，TIN方法從幾何原理的角度可以較為準(zhǔn)確地反映不同深度土壤銻的修復(fù)范圍；克里格插值(包括OK+TIN方法)和反距離插值受數(shù)據(jù)分布假設(shè)的約束，平滑效應(yīng)明顯，估計精度受限。

污染場地；銻；空間插值；修復(fù)范圍

近年來，隨著城市化進(jìn)程的加速，城市關(guān)停、搬遷工業(yè)企業(yè)后原址場地的污染問題凸顯，對生態(tài)環(huán)境、食品安全和人體健康構(gòu)成了嚴(yán)重威脅[1-3]?；谶@一現(xiàn)狀，國家和地方越來越重視污染場地的環(huán)境管理，相繼出臺場地相關(guān)法律法規(guī)及政策標(biāo)準(zhǔn)。按照《中華人民共和國土壤污染防治法》第三十五條的規(guī)定，土壤污染風(fēng)險管控和修復(fù)包括土壤污染狀況調(diào)查和土壤污染風(fēng)險評估、風(fēng)險管控、修復(fù)、風(fēng)險管控效果評估、修復(fù)效果評估、后期管理等活動。依據(jù)場地調(diào)查和風(fēng)險評估的結(jié)果，對于超過可接受風(fēng)險水平的點位需要建立其風(fēng)險管控范圍或修復(fù)范圍(下統(tǒng)稱為“修復(fù)范圍”)，為場地后續(xù)修復(fù)方案的建立及修復(fù)工程的實施提供依據(jù)。由于場地污染物及污染狀況的差異，從管理層面上國家尚未出臺修復(fù)范圍建立的技術(shù)規(guī)定，因此有必要開展這一方面的研究，以期為污染場地風(fēng)險管控及修復(fù)工作提供參考。

目前常用的土壤屬性制圖方法有地統(tǒng)計學(xué)插值模型、確定性插值模型、離散點的定性分析方法等[4-10]。地統(tǒng)計學(xué)插值模型是利用樣本點的統(tǒng)計規(guī)律，使樣本點之間的空間自相關(guān)定量化，從而在待預(yù)測點周圍構(gòu)建樣本點的空間結(jié)構(gòu)模型，比如克里格插值法[11]?？死锔袷菑淖兞肯嚓P(guān)性和變異性出發(fā)，在有限區(qū)域內(nèi)對區(qū)域化變量的取值進(jìn)行無偏、最優(yōu)估計的一種方法[12]。按照空間場是否存在漂移可將克里格插值分為普通克里格和泛克里格，其中普通克里格常稱作局部最優(yōu)線性無偏估計，當(dāng)數(shù)據(jù)服從二階平穩(wěn)假設(shè)時，可以采用普通克里格進(jìn)行插值；泛克里格是在區(qū)域化變量非平穩(wěn)條件下的一種最優(yōu)化無偏線性估值方法，它是從區(qū)域化變量的角度考慮數(shù)據(jù)的空間結(jié)構(gòu)性，內(nèi)插值更接近客觀實際[13]。確定性插值方法是基于信息點之間的相似程度或者整個曲面的光滑性來創(chuàng)建一個擬合曲面，如反距離加權(quán)平均插值法、趨勢面法、樣條函數(shù)法等。泰森多邊形是常用的空間分析方法之一[14-17]，它是由一組連接兩相鄰點線段的垂直平分線組成的連續(xù)多邊形組成。由于環(huán)境事件發(fā)生的空間集聚性，采用地統(tǒng)計方法進(jìn)行空間分析多數(shù)情況下會導(dǎo)致結(jié)果有偏，因此泰森多邊形也經(jīng)常用于反映區(qū)域環(huán)境因子的集聚效應(yīng)[8, 16-17]。

由于受歷史生產(chǎn)活動和污染源分布等因素的影響，工業(yè)污染場地土壤中的關(guān)注污染物往往呈現(xiàn)局部污染嚴(yán)重的特點，這些區(qū)域的污染物含量極高，具有很強的空間變異性[2, 18]。如何根據(jù)場地污染物空間分布選擇適用的空間插值方法是影響土壤修復(fù)范圍預(yù)測質(zhì)量的最主要因素之一?？死锔癫逯翟跀?shù)據(jù)服從二階平穩(wěn)假設(shè)時，預(yù)測精度更高；反距離插值在克里格插值無法取得可信賴的半變異函數(shù)時更占優(yōu)勢；泰森多邊形依據(jù)污染點位的空間位置進(jìn)行空間劃分，可以更好地反映區(qū)域的集聚效應(yīng)和空間結(jié)構(gòu)。研究表明，對于某種場地土壤污染物濃度呈偏態(tài)情況下，分區(qū)預(yù)測(克里格+TIN)的預(yù)測精度更高[19]。

銻是土壤中廣泛分布的一種有毒重金屬，其對人體及生物具有毒性及致癌性。中國的銻礦儲量和銻生產(chǎn)量為世界之最，也是場地上經(jīng)常出現(xiàn)的重金屬之一[20]，而環(huán)境中銻污染的報道卻較少。銻是子彈制作時常用的硬化劑，子彈大約含銻2% ～ 5%，研究表明，在打靶場廢棄的子彈被腐蝕后，分解釋放的銻會進(jìn)入土壤或者經(jīng)淋溶進(jìn)入土壤[21]。本研究選擇該場地土壤特征污染物銻作為研究對象，分析其含量的空間變異性，分別采用普通克里格插值(OK)、泛克里格插值(UK)、反距離插值(IDW)、泰森多邊形(TIN)、分區(qū)預(yù)測(OK+TIN)建立場地土壤銻修復(fù)范圍，從預(yù)測精度和修復(fù)范圍的角度評價幾種預(yù)測模型的適用性，以期為土壤銻污染狀況及修復(fù)范圍劃定方法的選擇提供參考。

1 材料與方法

1.1 研究區(qū)域概況

研究場地為我國某大型軍工場地，場地占地面積約27 hm2，按照當(dāng)?shù)匾?guī)劃的要求，企業(yè)目前已停產(chǎn)搬遷。企業(yè)有近70 a的生產(chǎn)歷史，生產(chǎn)初期工藝設(shè)備相對落后、經(jīng)營管理較為粗放、環(huán)保設(shè)施缺少或不完善，在多年生產(chǎn)經(jīng)營過程中，導(dǎo)致一些對人體有害物質(zhì)的排放。根據(jù)場地生產(chǎn)工藝及現(xiàn)場調(diào)查的需要，本次調(diào)查將場地劃分為A(污水排放區(qū))、B(底火、火工區(qū))、C(曳光彈生產(chǎn)區(qū))、D(射釘彈生產(chǎn)區(qū))4個區(qū)域(圖1、表1)。該企業(yè)主要生產(chǎn)原料包括雷汞、鉛、鋁粉、銅、鉀銻和石油烴等。場地地勢在總體上是東北、東南兩側(cè)向河谷逐漸降低，河谷部位則呈由西向東逐漸降低的趨勢。

區(qū)內(nèi)地貌依成因類型可劃分為堆積地形及河漫灘，河漫灘可分為高漫灘和低漫灘。場地土壤分為4層：第一層(①層)：人工堆積土，埋深為0 ～ 3.5 m，主要為低液限黏土夾爐灰、碎磚及建筑垃圾，將土壤巖性概化為低液限黏土進(jìn)行風(fēng)險評估；第二層(② ～ ③層)：老土層，埋深為3.5 ～ 5.5 m，主要為低液限黏土，將巖性統(tǒng)一概化為黏土；第三層(④ ～ ⑦層)：含水層，埋深為5.5 ～ 22.0 m，主要為級配不良中砂、粗砂組成，將巖性統(tǒng)一概化為砂土；第四層(⑧層)：含水層底板，埋深在22 m以下，為低液限黏土，將巖性統(tǒng)一概化為黏土。

該地塊土壤銻污染主要集中在埋深3.8 m以上的人工堆積層，隨著深度的增加土壤污染程度減輕，同時考慮到土壤巖性對污染物遷移風(fēng)險的影響，結(jié)合地質(zhì)分層(圖2)，本場地土壤銻風(fēng)險評估及修復(fù)范圍劃定將地層劃為0～1.0、1.0～3.8、>3.8 m。

(A區(qū)：污水排放區(qū)；B區(qū)：底火、火工區(qū)；C區(qū)：曳光彈生產(chǎn)區(qū)；D區(qū)：射釘彈生產(chǎn)區(qū))

表1 場地各分區(qū)主要生產(chǎn)設(shè)施及功能

圖2 場地土壤分層概化

1.2 樣品采集與分析

根據(jù)我國目前污染場地的管理要求，本次調(diào)查按照初步調(diào)查和詳細(xì)調(diào)查兩個步驟進(jìn)行。本次初步調(diào)查主要采用判斷布點，在可能造成潛在污染的熱點區(qū)域進(jìn)行采樣點布設(shè)。結(jié)合初步調(diào)查結(jié)果，進(jìn)行詳細(xì)調(diào)查的布點。整個調(diào)查區(qū)域通過開展初步和詳細(xì)調(diào)查共布置土壤采樣點159個，最大鉆探深度為31.5 m。

為了判斷不同深度土壤污染情況，按照地層分布進(jìn)行了不同深度的采樣，取樣深度范圍為0.2 ～ 31 m。本研究實際的鉆孔數(shù)目為159個，現(xiàn)場鉆探的同時采用RTK對各采樣點地面高程和地理位置進(jìn)行測量，對采樣點的坐標(biāo)和高程進(jìn)行記錄。場地土壤鉆探采用傳統(tǒng)的鋼索沖擊鉆法，獲取的土壤剖面樣品全部放置于木制巖心箱，表層土壤采樣時，先清除地表石塊、植被等，再進(jìn)行樣品采集。本次土壤調(diào)查共采集樣品數(shù)量為1 264個(含100個平行樣)，采集的樣品裝于200 ml的廣口玻璃瓶內(nèi)進(jìn)行保存。所有樣品采集完成后，置于硬質(zhì)塑料樣品轉(zhuǎn)運箱內(nèi)，加入藍(lán)冰以確保運輸過程中維持低溫。樣品采集完畢后送至實驗室進(jìn)行檢測。

土壤銻的檢測方法采用USEPA 6020B-2014電感耦合等離子體質(zhì)譜法進(jìn)行測定。樣品分析取0.1 ～ 0.2 g風(fēng)干土，過100目篩后，加入5 ml硝酸、2 ml高氯酸、2 ml氫氟酸，加蓋并放入160 ℃ 石墨爐消解，開蓋于約160 ℃ 加熱至透明液滴，加入0.5 ml硼酸溶液冷卻定容至50 ml待測。在樣品分析時進(jìn)行加標(biāo)回收率和平行樣的測定，標(biāo)樣(編號GSB04-1767-2004，批號174001)由國家有色金屬及電子材料分析測試中心提供，加標(biāo)樣品回收率為95% ～ 97%，現(xiàn)場平行樣的相對偏差為0.0% ～ 20.69%。

1.3 風(fēng)險控制值

根據(jù)健康風(fēng)險評估結(jié)果，確定銻為該場地土壤的目標(biāo)污染物，需基于健康風(fēng)險評估結(jié)果采用劑量–效應(yīng)模型反算其土壤風(fēng)險控制值?；谠O(shè)定為1的可接受非致癌風(fēng)險水平，通過風(fēng)險評估模型反推獲得其風(fēng)險控制值為20 mg/kg。我國現(xiàn)行《土壤環(huán)境質(zhì)量建設(shè)用地土壤污染風(fēng)險管控標(biāo)準(zhǔn)(試行)》(GB 36600—2018)[22]規(guī)定的第一類用地銻篩選值也為20 mg/kg，因此確定場地土壤銻的風(fēng)險控制值為20 mg/kg。

1.4 研究方法

1.4.1 土壤銻修復(fù)范圍預(yù)測方法 以土壤銻超標(biāo)地層為研究對象，結(jié)合土壤每層銻含量的分布特點，本研究選用以下5種方法分析土壤超標(biāo)銻修復(fù)范圍及精度：①普通克里格插值(OK)：基于普通克里格插值的“二階平穩(wěn)假設(shè)”，將原始數(shù)據(jù)進(jìn)行對數(shù)轉(zhuǎn)換后進(jìn)行插值；②泛克里格插值(UK)：通過數(shù)據(jù)分析，去除數(shù)據(jù)來自于背景值的主導(dǎo)趨勢，采用泛克里格插值對剔除趨勢后數(shù)據(jù)進(jìn)行插值；③反距離插值：采用原始數(shù)據(jù)進(jìn)行反距離插值；④泰森多邊形方法：作為污染場地常用修復(fù)范圍的確定方法，基于采樣點建立泰森多邊形，以超過修復(fù)目標(biāo)值的泰森多邊形作為修復(fù)范圍；⑤分區(qū)預(yù)測方法：基于正態(tài)分布的“3σ”原則，以銻含量超過平均值加4倍標(biāo)準(zhǔn)差的數(shù)值作為異常值，并以原始數(shù)據(jù)中值代替其高值部分，進(jìn)行拆除高值后數(shù)據(jù)的普通克里格插值(OK)，最后將兩部分插值數(shù)據(jù)疊加，將超過修復(fù)目標(biāo)值的區(qū)域劃定為修復(fù)范圍。

1.4.2 插值精度驗證 土壤污染物空間插值精度的驗證受驗證方法和驗證指標(biāo)的影響。本研究通過調(diào)整半方差模型的參數(shù)進(jìn)行交叉驗證，確定特定插值方法下的適用參數(shù)。交叉驗證經(jīng)常用來比較不同方法或同一方法不同參數(shù)選擇的插值結(jié)果，如可以通過半方差模型中選擇最佳搜索半徑和步長距離或鄰域點個數(shù)使克里格方差最小。

插值精度驗證指標(biāo)中平均誤差(ME)、均方根誤差(RMSE)為常用的衡量已知數(shù)據(jù)和預(yù)測數(shù)據(jù)差別的指標(biāo)，根據(jù)每個驗證點位置的觀測值和預(yù)測值計算，公式如下：

1.5 數(shù)據(jù)分析與處理

本研究土壤銻含量的描述性統(tǒng)計分析采用軟件SPSS 19.0，空間插值方法及精度驗證采用ESRI ArcGIS10.5進(jìn)行處理。

2 結(jié)果與分析

2.1 場地土壤銻含量的統(tǒng)計與分析

從表2中可以看出，場地0 ～ 1 m土壤銻含量范圍為0.40 ～ 507.39 mg/kg，平均值為23.64 mg/kg，中位數(shù)為5.50 mg/kg，銻平均含量超過我國《土壤環(huán)境質(zhì)量建設(shè)用地土壤污染風(fēng)險管控標(biāo)準(zhǔn)(試行)》(GB 36600—2018)第一類用地篩選值標(biāo)準(zhǔn)(20 mg/kg)，204個土壤樣品中有37個土壤銻超標(biāo)，超標(biāo)率為18.14%；1.0 ～ 3.8 m土壤銻含量范圍為0.10 ～ 260.83 mg/kg，平均值為3.78 mg/kg，418個樣品中有10個土壤銻超標(biāo)，超標(biāo)率為2.39%。與0 ～ 1.0 m土層相比，1.0 ～ 3.8 m層土壤銻含量最大值和平均值均呈現(xiàn)下降趨勢。

0 ～ 1.0 m和1.0 ～ 3.8 m土壤銻含量均呈偏態(tài)分布(表2，圖3)，其偏度分別為5.15、14.46，峰度分別為30.04、246.38，極差分別為506.99、260.73 mg/kg，說明土壤銻含量具有明顯的集聚性。基于OK+TIN方法識別的0 ～ 1.0 m土壤銻含量的異常值有4個，1.0 ～ 3.8 m土壤銻含量的異常值有2個，分別用其中位數(shù)替代。0 ～ 1.0 m、1.0 ～ 3.8 m替換前后的土壤銻含量對數(shù)正態(tài)累積頻率如圖4所示。從圖中可以看出，中值替換對每層的數(shù)據(jù)分布基本沒有影響，0 ～1.0 m土壤銻含量基本均位于95% 置信區(qū)間內(nèi)，數(shù)據(jù)呈對數(shù)正態(tài)分布；而1.0 ～ 3.8 m土壤銻含量均不服從對數(shù)正態(tài)分布，原因在于1.0 ～ 3.8 m土壤銻濃度整體較低，其超標(biāo)率為2.4%，而最大超標(biāo)倍數(shù)為12，說明1.0 ～ 3.8 m土壤銻含量的極端效應(yīng)十分明顯。

表2 場地土壤銻含量統(tǒng)計特征

注：①《土壤環(huán)境質(zhì)量建設(shè)用地土壤污染風(fēng)險管控標(biāo)準(zhǔn)(試行)》(GB36600—2018)第一類用地篩選值。

圖3 場地土壤銻含量直方圖

圖4 場地土壤銻含量對數(shù)正態(tài)累積頻率圖

2.2 場地土壤銻含量空間分布

將土壤各采樣點銻濃度與原廠區(qū)平面布置圖進(jìn)行空間疊加后可以發(fā)現(xiàn)(圖5)，0 ～ 1.0、1.0 ～ 3.8 m土壤銻含量具有明顯的空間差異性，以0 ～ 1.0 m土壤尤為明顯。0 ～ 1.0 m土壤銻含量超過篩選值特別是超標(biāo)嚴(yán)重的采樣點主要分布在以銻為原料的生產(chǎn)與加工區(qū)域，即場區(qū)西南部的曳光彈生產(chǎn)區(qū)、射釘彈生產(chǎn)區(qū)和東中部的底火火工區(qū)。從圖5中可以看出，土壤銻超標(biāo)嚴(yán)重的采樣點分布在B、C兩個分區(qū)，銻最高濃度507.39 mg/kg出現(xiàn)在曳光彈生產(chǎn)區(qū)域0 ～ 1.0 m地層，1.0 ～ 3.8 m土壤銻含量最大值出現(xiàn)在C區(qū)底火加工區(qū)。原因在于這些區(qū)域都屬于場地的子彈生產(chǎn)區(qū)，其生產(chǎn)過程中采用的起爆藥含有硫化銻成分，生產(chǎn)過程中炸藥的遺撒和污水排放都會導(dǎo)致銻在表層土壤累積。

圖5 場地土壤銻含量空間分布

3 土壤銻修復(fù)范圍預(yù)測

3.1 模型擬合

同一插值方法中，不同參數(shù)的設(shè)定對插值精度有顯著影響。影響IDW插值的因素有指數(shù)值的選擇、估計的鄰域樣點數(shù)目及搜索半徑；影響克里格插值的因素包括數(shù)據(jù)的變異性和空間結(jié)構(gòu)、半方差函數(shù)模型和搜索半徑，這些方法大部分參數(shù)的確定沒有太多可遵循的規(guī)律，需通過交叉驗證等方法來確定。

3.1.1 IDW擬合參數(shù) 本研究通過選取不同的指數(shù)值進(jìn)行IDW插值，鄰域樣點數(shù)據(jù)和搜索半徑根據(jù)場地生產(chǎn)工藝分布確定合理的數(shù)值，對插值結(jié)果進(jìn)行交叉驗證并計算均方根誤差(RMSE)。由圖6可以看出，RMSE隨著指數(shù)值增大呈現(xiàn)先減小后增大的趨勢，0 ～ 1.0 m土壤銻含量RMSE值在=1.88處達(dá)到最小值，1.0 ～ 3.8 m土壤銻含量則在=9.59處達(dá)到最小值，此時其IDW插值效果最優(yōu)，因此本研究分別選取=1.88及=9.59進(jìn)行IDW插值。

3.1.2 土壤銻含量半變異函數(shù) 普通克里格和泛克里格插值擬合的關(guān)鍵是變量的變異函數(shù)分析，因此對0 ～ 1.0 m和1.0 ～ 3.8 m土壤銻含量整體數(shù)據(jù)經(jīng)對數(shù)轉(zhuǎn)換后的數(shù)據(jù)、去掉主導(dǎo)趨勢的對數(shù)轉(zhuǎn)換數(shù)據(jù)和拆分后待插值數(shù)據(jù)分別進(jìn)行半變異函數(shù)分析。塊金值和基臺值是半變異函數(shù)的兩個重要的參數(shù)，塊金值反映的是最小抽樣尺度以下變量的變異性及測量誤差，塊金值與基臺值的比值表示隨機部分引起的空間異質(zhì)性占系統(tǒng)總變異的比例。如果該比值高，說明樣本間的變異更多的是由隨機因素引起的，經(jīng)交叉驗證確定的最優(yōu)模型及參數(shù)如表3。從表中可以看出，穩(wěn)態(tài)模型可較好地模擬3種情形下不同深度土壤銻含量的空間變異。0 ～ 1.0 m土壤銻數(shù)據(jù)去掉主導(dǎo)趨勢后半變異模型參數(shù)與整體數(shù)據(jù)保持一致，說明0 ～ 1.0 m土壤銻含量主要受場地活動的影響，土壤背景含量影響有限；拆分后數(shù)據(jù)模型基臺值降低，塊金效應(yīng)增加，說明場地土壤銻含量空間變異減少，空間相關(guān)性增加。1.0 ～ 3.8 m土壤銻數(shù)據(jù)去掉主導(dǎo)趨勢后塊金效應(yīng)增加明顯，說明1.0 ～ 3.8 m土壤銻主要受場地土壤隨機因素的影響，來自于場地生產(chǎn)活動的影響減少。

圖6 指數(shù)p值對IDW插值RMSE的影響

表3 半變異函數(shù)模型參數(shù)

3.2 不同方法對土壤銻修復(fù)范圍的預(yù)測精度分析

表4列出了各個插值方法預(yù)測統(tǒng)計結(jié)果與預(yù)測誤差。從表中可以看出，0 ～ 1.0 m土壤深度OK、UK、IDW、OK+TIN預(yù)測銻含量最大值為160.70、160.69、238.48、41.72 mg/kg，最小值為2.47、2.47、1.15、2.84 mg/kg，平均值為23.36、27.84、28.30、18.23 mg/kg。克里格插值對于高值和低值的平滑作用特別明顯，由于去除了極端高值的影響，數(shù)據(jù)明顯呈現(xiàn)左側(cè)拖尾效應(yīng)，從而導(dǎo)致OK+TIN的平滑作用更為明顯，估計的平均值也明顯低于其他模型。從預(yù)測的平均誤差和均方根誤差來看，OK+TIN的分區(qū)預(yù)測方法要優(yōu)于其他3種插值方法，說明OK+TIN對于污染物空間分布趨勢有更好的精度，但由于正態(tài)分布轉(zhuǎn)換過程中引起的平滑效應(yīng)，導(dǎo)致其高值的估計誤差太大，因此修復(fù)范圍確定的誤差較大。IDW方法對樣點分布位置和樣點數(shù)量要求較高，樣點越多且分布越均勻其預(yù)測精度要高，而場地由于受建筑物、污染識別環(huán)節(jié)及成本的影響，很難做到高密度和均勻網(wǎng)格布點采樣，因此插值方法對于污染范圍的預(yù)測精度較低。

從表4中可以看出1.0 ～ 3.8 m土壤深度 OK、UK、IDW、OK+TIN預(yù)測銻含量最大值為26.72、29.16、57.42、18.89 mg/kg，預(yù)測的最小值1.68、1.69、0.26、1.50 mg/kg，平均值為6.38、7.70、6.97、5.04 mg/kg。4種插值方法對于高值的平滑作用均超過10倍，這可能與土壤銻含量不服從對數(shù)正態(tài)分布有關(guān)。OK、UK、IDW 3種插值方法的均方根誤差差別不大。盡管OK+TIN具有相對低的均方根誤差，但是由于本地層土壤大部分采樣點銻含量較低，加上OK+TIN對于高值的平滑作用，導(dǎo)致插值后土壤銻含量均低于風(fēng)險控制值20 mg/kg，這進(jìn)一步表明OK+TIN插值方法跟OK方法一樣受限于數(shù)據(jù)分布的要求。

3.3 不同方法對于土壤銻修復(fù)范圍預(yù)測結(jié)果比較

對原始數(shù)據(jù)進(jìn)行了反距離插值，同時對經(jīng)對數(shù)轉(zhuǎn)換的數(shù)據(jù)進(jìn)行普通克里格和泛克里格插值，將不同方法插值后超過修復(fù)目標(biāo)值20 mg/kg的區(qū)域劃定為修復(fù)范圍，如圖7所示。從圖中可以看出，針對0 ～ 1.0 m土壤銻修復(fù)范圍，普通克里格和泛克里格方法確定的修復(fù)范圍一致，進(jìn)一步驗證了土壤銻主要來自于場地的生產(chǎn)活動。OK、UK、IDW對于低值的平滑效應(yīng)也特別明顯，未超標(biāo)點位被納入修復(fù)范圍，從而帶來修復(fù)規(guī)模的增加。與其他方法相比，OK+TIN方法可以更好地刻畫污染嚴(yán)重區(qū)域的修復(fù)范圍，但對于周圍大部分清潔點存在的情況下，超標(biāo)點極易被平滑為未超標(biāo)點而導(dǎo)致修復(fù)范圍的縮小。以土壤鉆孔建立泰森多邊形，并將超過20 mg/kg的點位所在的泰森多邊形劃為修復(fù)范圍，表5統(tǒng)計了各種方法確定的修復(fù)面積。從表中可以看出，各方法確定的修復(fù)面積大小順序為IDW>OK=UK>TIN>OK+TIN，修復(fù)面積依次為74 154、57 427、37 338、32 707 m2，表明由于受采樣方案的限制，IDW方法平滑效益最明顯，將會帶來修復(fù)成本的增加；OK+TIN與TIN 確定的修復(fù)范圍接近，但OK+TIN仍存在一定的平滑效應(yīng)。從1.0 ～ 3.8 m各插值方法確定的場地土壤銻的修復(fù)范圍如圖8，從圖中可以看出，OK、UK、IDW方法對于高值的平滑作用較為明顯，以O(shè)K、OK+TIN尤為突出，OK+TIN僅保留了TIN的部分，OK插值導(dǎo)致了8個超標(biāo)點均平滑掉，從而導(dǎo)致修復(fù)范圍的減少，OK確定的修復(fù)面積也最小(表5)；IDW方法相比較其他兩種方法較為合理，但其對于未超標(biāo)點的平滑作用帶來修復(fù)面積的增加。

表4 不同插值方法預(yù)測統(tǒng)計結(jié)果與預(yù)測誤差(mg/kg)

注：平均誤差為數(shù)據(jù)的絕對值。

從不同深度土壤銻的模型預(yù)測精度和修復(fù)范圍可以看出，克里格插值由于其對于高值和低值具有較強的平滑效應(yīng)，預(yù)測結(jié)果精度較低；反距離插值由于受場地數(shù)據(jù)布點方法的影響，插值結(jié)果存在一定的局限性；OK+TIN在一定程度上彌補了克里格插值的平滑效應(yīng)，但其對于高值的平滑作用仍存在，仍可能存在遺漏高風(fēng)險區(qū)域范圍的情況。從插值效果來看，考慮空間相關(guān)性的克里格插值對于數(shù)據(jù)質(zhì)量要求較高，在數(shù)據(jù)不符合假設(shè)條件的情況下，容易帶來修復(fù)面積的增加或者遺漏，這表明了地統(tǒng)計學(xué)插值方法對于樣本數(shù)量和質(zhì)量的依賴性很難滿足污染物空間分布預(yù)測對于精度的要求[23-24]。泰森多邊形為反距離插值的一種極端情況，通過以土壤鉆孔代表特定修復(fù)地塊的方法，從幾何原理上給出了較為合理的修復(fù)范圍，其優(yōu)勢在于確定樣點較少的深層土壤修復(fù)范圍時，可以以鉆孔代表的地塊為修復(fù)范圍，從而保持了不同深度超標(biāo)點位修復(fù)范圍的一致性。

4 結(jié)論與建議

受場地生產(chǎn)活動的影響，場地土壤銻含量呈明顯的空間變異性，在污染源區(qū)域銻含量呈現(xiàn)明顯的集聚性。0 ～ 1.0 m深度土壤銻含量主要受場地活動影響，基本接近對數(shù)正態(tài)分布；1.0 ～ 3.8 m土壤銻含量僅個別點位受場地活動影響，呈偏態(tài)分布。

采用普通克里格插值(OK)、泛克里格插值(UK)、分區(qū)預(yù)測(OK+TIN)、反距離插值(IDW)、泰森多邊形(TIN)5種模型確定場地土壤需修復(fù)的范圍。OK、UK方法可以結(jié)合數(shù)據(jù)的空間分布確定預(yù)測精度和預(yù)測誤差，但其對高值和低值的平滑作用特別明顯，容易導(dǎo)致超標(biāo)點位被平滑掉或未超標(biāo)點納入修復(fù)范圍，從而導(dǎo)致修復(fù)范圍的減少或增加。OK+TIN本質(zhì)上是克里格插值和TIN的結(jié)合，其通過調(diào)整數(shù)據(jù)分布進(jìn)行克里格插值，因此估計精度較高，但其在大部分點位都是清潔點(1.0 ～ 3.8 m)的情況下很難克服OK的局限，容易導(dǎo)致修復(fù)范圍的減少甚至遺漏。IDW的估計精度受限于樣點分布位置和樣本數(shù)量，在數(shù)據(jù)不均勻的情況下容易將未超標(biāo)點納入修復(fù)范圍。在受成本和污染物空間分布制約的情況下，建議采用TIN方法確定污染物超標(biāo)范圍，其從幾何原理的角度可以較為準(zhǔn)確地劃定不同深度土壤的修復(fù)范圍，其突出優(yōu)勢在于可以實現(xiàn)不同深度超標(biāo)點位修復(fù)范圍的一致性且可操作性強。盡管TIN方法無法提供估計精度，但從業(yè)者可以結(jié)合地塊的污染識別情況和增加采樣點的方法進(jìn)一步優(yōu)化修復(fù)邊界，進(jìn)而較為準(zhǔn)確地確定污染物的修復(fù)范圍。

圖7 不同插值方法對0 ～ 1.0 m土壤銻修復(fù)范圍的預(yù)測

圖8 不同插值方法對1.0 ～ 3.8 m土壤銻修復(fù)范圍的預(yù)測

表5 不同方法確定的土壤銻修復(fù)范圍(m2)

污染場地調(diào)查中，準(zhǔn)確地反映場地污染修復(fù)范圍十分重要。由于場地的復(fù)雜性和成本的約束，調(diào)查數(shù)據(jù)往往很難滿足克里格插值或其變身(OK+TIN)及IDW的假設(shè)條件，因此其確定的修復(fù)范圍往往與實際情況不相符。采用幾何(TIN)和布點優(yōu)化相結(jié)合的方法不僅可以準(zhǔn)確地劃定場地修復(fù)范圍，而且對污染地塊的修復(fù)有著重要的意義。

[1] 駱永明. 中國污染場地修復(fù)的研究進(jìn)展、問題與展望[J]. 環(huán)境監(jiān)測管理與技術(shù), 2011, 23(3): 1–6.

[2] 駱永明, 滕應(yīng). 我國土壤污染的區(qū)域差異與分區(qū)治理修復(fù)策略[J]. 中國科學(xué)院院刊, 2018, 33(2): 145–152.

[3] 廖曉勇, 崇忠義, 閻秀蘭, 等. 城市工業(yè)污染場地: 中國環(huán)境修復(fù)領(lǐng)域的新課題[J]. 環(huán)境科學(xué), 2011, 32(3): 784–794.

[4] 史文嬌, 岳天祥, 石曉麗, 等. 土壤連續(xù)屬性空間插值方法及其精度的研究進(jìn)展[J]. 自然資源學(xué)報, 2012, 27(1): 163–175.

[5] 馬宏宏, 余濤, 楊忠芳, 等. 典型區(qū)土壤重金屬空間插值方法與污染評價[J]. 環(huán)境科學(xué), 2018, 39(10): 4684–4693.

[6] 謝龍濤, 潘劍君, 白浩然, 等. 基于GIS的農(nóng)田土壤重金屬空間分布及污染評價——以南京市江寧區(qū)某鄉(xiāng)鎮(zhèn)為例[J]. 土壤學(xué)報, 2020, 57(2): 316–325.

[7] 謝云峰, 杜平, 曹云者, 等. 基于地統(tǒng)計條件模擬的土壤重金屬污染范圍預(yù)測方法研究[J]. 環(huán)境污染與防治, 2015, 37(1): 1–6.

[8] 謝紹鋒, 肖化順, 儲蓉, 等. 基于泰森多邊形的廣州市林火空間分布規(guī)律研究[J]. 西北林學(xué)院學(xué)報, 2018, 33(3): 178–185.

[9] 謝云峰, 陳同斌, 雷梅, 等. 空間插值模型對土壤Cd污染評價結(jié)果的影響[J]. 環(huán)境科學(xué)學(xué)報, 2010, 30(4): 847–854.

[10] 劉庚, 郭觀林, 南鋒, 等. 某大型焦化企業(yè)污染場地中多環(huán)芳烴空間分布的分異性特征[J]. 環(huán)境科學(xué), 2012, 33(12): 4256–4262.

[11] Bilonick R A. An introduction to applied geostatistics[J]. Technometrics, 2012, 33(4): 483–485.

[12] Breikreuz H. Handbook of applied advanced geostatistical ore reserve estimation[J]. Earth-Science Reviews, 1992, 32(3): 222–223.

[13] Bilonick R A. An introduction to applied geostatistics[J]. Technometrics, 2012, 33(4): 483–485.

[14] She B, Zhu X Y, Ye X Y, et al. Weighted network Voronoi Diagrams for local spatial analysis[J]. Computers, Environment and Urban Systems, 2015, 52: 70–80.

[15] Xiong C K, Chen D F, Lu D, et al. Path planning of multiple autonomous marine vehicles for adaptive sampling using Voronoi-based ant colony optimization[J]. Robotics and Autonomous Systems, 2018, 115: 90–103.

[16] Gahegan M, Lee I. Data structures and algorithms to support interactive spatial analysis using dynamic Voronoi diagrams[J]. Computers, Environment and Urban Systems, 2000, 24(6): 509–537.

[17] Barclay M J, Galton A. Comparison of region approximation techniques based on Delaunay triangulations and Voronoi diagrams[J]. Computers, Environment and Urban Systems, 2008, 32(4): 261–267.

[18] 譚萬能, 李志安, 鄒碧, 等. 地統(tǒng)計學(xué)方法在土壤學(xué)中的應(yīng)用[J]. 熱帶地理, 2005, 25(4): 307–311.

[19] 劉庚, ?？〗? 張朝, 等. 某鉛酸蓄電池污染場地表層土壤重金屬Pb空間分布預(yù)測研究[J]. 環(huán)境科學(xué), 2014, 35(12): 4712–4719.

[20] 駱永明, 滕應(yīng). 中國土壤污染與修復(fù)科技研究進(jìn)展和展望[J]. 土壤學(xué)報, 2020, 57(5): 1137–1142.

[21] Annette J C, Hermann M, Paul W, et al. Solubility of antimony and other elements in samples taken from shooting ranges[J]. Journal of Environmental Quality, 2005, 34(1): 248–54.

[22] 中華人民共和國生態(tài)環(huán)境部. 土壤環(huán)境質(zhì)量建設(shè)用地土壤污染風(fēng)險管控標(biāo)準(zhǔn)(試行)(GB36600—2018)[S]. 北京: 標(biāo)準(zhǔn)出版社, 2018.

[23] 王幼奇, 張興, 趙云鵬, 等. 基于GIS和地理加權(quán)回歸的砂田土壤陽離子交換量空間預(yù)測[J]. 土壤, 2020, 52(2): 421–426.

[24] 陳思明, 王寧, 秦艷芳, 等. 基于特征變量與支持向量機回歸克里格(SVRK)法的濕地土壤有機質(zhì)空間變異特征分析[J]. 土壤, 2020, 52(6): 1298–1305.

Remediation Area Prediction of Heavy Metal Sb in Soil Based On Spatial Interpolation

JIA Lin1,2,3, XIA Tianxiang1,2,3*, ZHANG Lina1,2,3, JIA Xiaoyang1,2,3, ZHANG Dan1,2,3

(1 Beijing Municipal Research Institute of Ecological Environmental Protection, Beijing 100037, China; 2 National Engineering Research Centre of Urban Environmental Pollution Control, Beijing 100037, China; 3 Beijing Key Laboratory for Risk Modeling and Remediation of Contaminated Sites, Beijing 100037, China)

Based on the distribution characteristics of Sb content in each layer of soil, the prediction accuracy and remediation of Sb contamination in each soil layer were analyzed and compared by Kriging interpolation model (OK, UK), inverse distance interpolation (IDW), Voronoi diagram (TIN) and zone prediction model (OK + TIN). The results showed that the order of remediation area determined by each method in 0-1 m soil layer was IDW>OK =UK>TIN>OK+TIN and the remediation area was 74 154, 57 427, 37 338 and 32 707 m2respectively; the order of remediation area determined by each method in 1-3.8 m soil layer was IDW>UK>TIN>OK+TIN>OK, the remediation area was 13 089, 11 030, 9 660, 2 183 and 1 231 m2respectively. The smoothing effect brought by Kriging interpolation resulted in the omission of area exceeding 20 mg/kg and inclusion of area ≤20 mg/kg in the prediction of remediation area. Under the constraints of cost and spatial distribution of pollutants, TIN method could accurately reflect the remediation range of soil Sb within different depths from the perspective of geometric principles; Kriging interpolation including OK+TIN method and IDW interpolation were constrained by the data distribution assumption, with obvious smoothing effect and limited estimation accuracy. The above results have certain guiding significance for the determination of remediation scope of contaminated sites.

Contaminated site; Antimony (Sb); Spatial interpolation; Remediation scope

X131.3；X53

A

10.13758/j.cnki.tr.2022.04.021

賈琳, 夏天翔, 張麗娜, 等. 基于不同空間插值的污染場地土壤銻修復(fù)范圍預(yù)測研究. 土壤, 2022, 54(4): 817–826.

國家重點研發(fā)計劃項目(2018YFC1801405，2018YFC1801401, 2020YFC180750)）資助。

(xiatianxiang@cee.cn)

賈琳(1980—)，女，山東萊陽人，博士，副教授，主要從事污染場地調(diào)查與修復(fù)研究。E-mail: jl_0706@163.com