郝易成，師華定，白中科，楊澤棟，秦治恒

（1.中國環(huán)境科學(xué)研究院土壤與固廢環(huán)境研究所，北京 100012；2.中國地質(zhì)大學(xué)土地科學(xué)技術(shù)學(xué)院，北京 100083）

土壤調(diào)查是獲取土壤環(huán)境狀況與污染分布最重要的手段，調(diào)查結(jié)果的精度直接關(guān)系到污染物風(fēng)險評價的準(zhǔn)確性和管理者對于土壤環(huán)境決策管理的合理性[1]。其中，點(diǎn)位布設(shè)、樣品采集、含量分析等各個環(huán)節(jié)都會導(dǎo)致調(diào)查結(jié)果的不確定性。土壤本身具有較強(qiáng)的空間異質(zhì)性，僅提高化驗(yàn)分析與樣品采集的精度難以提高土壤污染調(diào)查結(jié)果的精度。研究表明，針對空間異質(zhì)性大的要素，點(diǎn)位布設(shè)方案是影響調(diào)查結(jié)果精度的最重要因素之一[1-3]。Gustavsson等[3]在土壤污染調(diào)查誤差研究中發(fā)現(xiàn)，抽樣誤差是分析誤差的20倍以上。Jenkins等[4]在土壤三硝基甲苯污染的研究中表明，至少95%的變異度是由采樣位置導(dǎo)致。因此，優(yōu)化布點(diǎn)方案對提高土壤污染調(diào)查結(jié)果的精度有著重要作用。

目前，布點(diǎn)方法主要有均勻布點(diǎn)與非均勻布點(diǎn)兩種[5]。均勻布點(diǎn)是在缺乏背景信息的情況下，采取諸如隨機(jī)布點(diǎn)、網(wǎng)格布點(diǎn)等方法布點(diǎn)。而非均勻布點(diǎn)則是在初步調(diào)查的基礎(chǔ)上，根據(jù)已有背景信息，針對性布點(diǎn)?，F(xiàn)有的土壤污染調(diào)查大多缺乏先驗(yàn)背景信息，以均勻布點(diǎn)居多，多是對研究區(qū)內(nèi)污染物平均水平的估計[6]。而在土壤污染與防治的管理過程中，通常關(guān)注于污染位置、范圍與程度。因此，以平均水平估計為目的的均勻布點(diǎn)難以滿足管理決策需求[7-8]。而在實(shí)際的土壤污染調(diào)查中，一般分為普查和詳查等階段。普查主要目的是摸清土壤基本情況，識別土壤污染物與大致污染范圍，其布點(diǎn)范圍廣、分布均勻、點(diǎn)位較少。詳查則是在普查提供的背景信息的基礎(chǔ)上，有針對性地在潛在污染區(qū)增加樣點(diǎn)，更精確地確定污染類型、污染程度與污染范圍，其布點(diǎn)針對性強(qiáng)、區(qū)域集中、點(diǎn)位較多。土壤污染布點(diǎn)方案的誤差，主要來源于清潔區(qū)被高估和污染區(qū)被低估[9-10]。前者會增加不必要的修復(fù)投入和管理成本，后者則會使污染區(qū)不能得到有效的管理與治理。為了降低該類誤差，通常會增加樣本量，隨之而來也增加了成本。高效的布點(diǎn)方案就是要盡可能在降低誤差的同時減少成本投入，以最優(yōu)的資金效率完成土壤污染調(diào)查。布點(diǎn)方案優(yōu)化的目的就是尋找降低調(diào)查誤差的最優(yōu)布點(diǎn)區(qū)間與樣本量。土壤污染同時受氣候、地質(zhì)、水文、生態(tài)等背景條件，污染源匯關(guān)系以及污染物性質(zhì)等因素的綜合影響，其分布表現(xiàn)出不同的空間相關(guān)性和變異性，能否準(zhǔn)確判斷污染物分布是影響調(diào)查結(jié)果的關(guān)鍵部分[1]。

近年來，應(yīng)用地統(tǒng)計學(xué)的方法進(jìn)行土壤污染調(diào)查布點(diǎn)與加密布點(diǎn)已成為研究熱點(diǎn)[1，11-17]。該方法基于土壤污染的空間自相關(guān)性，運(yùn)用空間插值或條件模擬的方法，預(yù)測土壤污染空間分布，優(yōu)化點(diǎn)位布局，提高布點(diǎn)效率[9，11，15，18-20]。地統(tǒng)計學(xué)中各類方法應(yīng)用于加密布點(diǎn)的研究已有很多，但大多集中于不同方法在加密布點(diǎn)應(yīng)用中的可行性研究，而對基于某一方法的不同布點(diǎn)方案的效率（即在固定投入如點(diǎn)位個數(shù)情況下，采用不同布點(diǎn)方案帶來的調(diào)查結(jié)果準(zhǔn)確性的高低）變化的研究依然十分匱乏。

指示克里金（Indicator kriging）是由 Journel（1982）提出的，最初被用于地質(zhì)統(tǒng)計學(xué)處理偏態(tài)分布和極值的技術(shù)，常用于礦床礦石品位評估。由于土壤污染數(shù)據(jù)常具有極端值和高偏離度的特性，指示克里金也被廣泛應(yīng)用于土壤污染的研究[21-24]。該方法在普查數(shù)據(jù)的基礎(chǔ)上，通過指示克里金模擬來預(yù)測污染物含量超過閾值的概率。指示克里金是一種非參數(shù)分布估計方法，該方法基于簡單的二元變換，其數(shù)據(jù)的0～1指示變換使得預(yù)測變量對離群值具有的魯棒性，不存在均勻分布的假設(shè)[25]。指示克里金預(yù)測值的大小代表大于或小于特定閾值的概率，即從指示數(shù)據(jù)得出的非抽樣地點(diǎn)的期望值等于變量的累計分布[26]。

本研究即聚焦于不同布點(diǎn)方案（加密區(qū)域、加密程度）的效率問題和方法優(yōu)化，以土壤污染詳查工作的需求為導(dǎo)向，采用指示克里金的方法確定加密布點(diǎn)區(qū)域，對該方法下不同概率區(qū)間布點(diǎn)的效率進(jìn)行比較分析，并對在最優(yōu)區(qū)間內(nèi)的最優(yōu)樣本量進(jìn)行了研究，以提高土壤污染布點(diǎn)方案對污染物分布的估計精度，為土壤污染調(diào)查提供科學(xué)方法的支持。

1 材料與方法

1.1 研究區(qū)概況

本研究以高原區(qū)某山區(qū)為例。考慮到樣點(diǎn)代表性，選取海拔4200 m以下人類活動較多的區(qū)域?yàn)榘咐齾^(qū)。案例區(qū)南北跨度25 km，東西跨度40 km，總面積484.01 km2。該地區(qū)為灌叢草甸區(qū)，以牧草地為主要地類。空氣稀薄、干燥，雨季集中在6—9月，年降水量約450 mm，年均氣溫7.4℃，屬高原溫帶半干旱季風(fēng)氣候區(qū)。土壤類型主要為高山草甸土，兼有山地灌叢草原土，土體礫石含量高，土層較薄，肥力較低。

研究區(qū)內(nèi)有鉛礦開采企業(yè)，礦區(qū)尾礦中含有大量鉛，雨季地表徑流沖刷嚴(yán)重，污染物隨雨水、河水?dāng)U散，導(dǎo)致周邊土壤存在較大的鉛污染可能性。

1.2 樣品采集與測試

本案例共采集樣點(diǎn)588個，參考500 m×500 m的網(wǎng)格布點(diǎn)，每個樣點(diǎn)均采用五點(diǎn)法采樣，采樣深度20 cm。為防止采樣工具對樣品造成污染，用不銹鋼鏟挖好的土方需用木鏟對表面剔除，并用木鏟取樣。

土壤樣品鉛含量測試：經(jīng)消解預(yù)處理后，采用電感耦合等離子體質(zhì)譜儀進(jìn)行檢測，根據(jù)元素的質(zhì)譜圖或特征離子進(jìn)行定性，內(nèi)標(biāo)法定量。

1.3 方法

1.3.1 基于指示克里金的加密布點(diǎn)區(qū)域確定方法

在決策管理中，需要治理與否僅取決于污染物含量是否超過某一閾值，基于指示克里金的加密布點(diǎn)方法正是基于該需求提出。

指示克里金[25]首先將污染物濃度值轉(zhuǎn)換為基于某閾值的指示值，如式（1）。假設(shè)在研究區(qū)A內(nèi)，設(shè)置污染物含量閾值為ZC，則在A內(nèi)每一個樣點(diǎn)X∈A上定義一個Zc的如下階梯函數(shù)：

然后預(yù)測單點(diǎn)位置或區(qū)域的累積概率，如式（2）和式（3）。

將λi帶入公式（2）或公式（3）即可得到不超過Zc的累積概率。該方法的誤差可以通過式（4）計算[27]。

累積概率的取值范圍為0～1。某點(diǎn)位概率越高，其污染的可能性也就越大。污染與否的誤判主要來源于污染區(qū)與清潔區(qū)的過渡區(qū)域[7，28]，且為了準(zhǔn)確估計污染區(qū)范圍，該過渡區(qū)域應(yīng)更靠近污染區(qū)。過渡區(qū)域在概率地圖上即表現(xiàn)為某概率區(qū)間。只有在該概率區(qū)間內(nèi)加密布點(diǎn)才能提高對污染范圍的估計精度。假設(shè)過渡區(qū)域的概率區(qū)間為0.50～0.95，即已經(jīng)認(rèn)為概率大于0.95的范圍為污染區(qū)，概率小于0.50的范圍為清潔區(qū)。概率在0.50～0.95的范圍即為加密布點(diǎn)區(qū)，需要進(jìn)一步調(diào)查。

1.3.2 數(shù)據(jù)處理方法

運(yùn)用指示克里金方法進(jìn)行概率制圖。對污染范圍的預(yù)測采用序貫高斯模擬方法，該方法是貝葉斯理論的應(yīng)用，消除了普通插值的平滑效應(yīng)。運(yùn)用序貫高斯進(jìn)行面積預(yù)測時，單次模擬存在一定的波動性，需要大量模擬后取均值實(shí)現(xiàn)[11]。借鑒已有相關(guān)研究經(jīng)驗(yàn)[11]，以及不同模擬次數(shù)的實(shí)驗(yàn)（500次以上模擬結(jié)果基本不會發(fā)生變化），為保證研究的準(zhǔn)確性，本研究模擬次數(shù)為1000次。指示克里金概率制圖與序貫高斯模擬預(yù)測均在ArcGIS 10.2中實(shí)現(xiàn)。

2 案例驗(yàn)證及效率分析

以該區(qū)域鉛污染為例，進(jìn)行土壤污染調(diào)查加密布點(diǎn)區(qū)域優(yōu)化及效率研究。該案例為模擬研究，以500 m×500 m網(wǎng)格輔助參考，一次性采集588個點(diǎn)位，假設(shè)以這588個點(diǎn)位近似代表案例區(qū)土壤真實(shí)情況。運(yùn)用類似抽樣調(diào)查的方法，在研究區(qū)上疊加1 km×1 km的網(wǎng)格，落在網(wǎng)格內(nèi)的點(diǎn)作為普查點(diǎn)位，當(dāng)網(wǎng)格內(nèi)有多個點(diǎn)位時，隨機(jī)選取一個作為普查點(diǎn)位，由此獲得普查點(diǎn)位共176個。參照《土壤環(huán)境質(zhì)量農(nóng)用地土壤污染風(fēng)險管控標(biāo)準(zhǔn)（試行）》（GB 15618—2018），運(yùn)用序貫高斯模擬分別對總體點(diǎn)位、普查點(diǎn)位、加密后點(diǎn)位進(jìn)行污染范圍預(yù)測，以通過總體點(diǎn)位預(yù)測的污染面積為基準(zhǔn)，分別將通過普查點(diǎn)位和加密后點(diǎn)位預(yù)測的污染面積與總體點(diǎn)位預(yù)測的污染面積相比較，以此作為普查點(diǎn)位相對總體點(diǎn)位的精度與加密后點(diǎn)位相對總體點(diǎn)位的精度，并將二者進(jìn)行比較分析，研究不同概率區(qū)間內(nèi)布點(diǎn)效率與最優(yōu)概率區(qū)間內(nèi)樣本量的效率。

式（5）為精度計算方法。P1、P2分別表示普查精度和加密布點(diǎn)后精度；A0、A1、A2分別表示總體點(diǎn)位、普查點(diǎn)位、加密布點(diǎn)后點(diǎn)位預(yù)測的污染面積。

式（6）為布點(diǎn)效率計算方法。E表示某一概率區(qū)間內(nèi)加密布點(diǎn)的平均效率；n表示該概率區(qū)間內(nèi)加密布點(diǎn)的數(shù)量。

2.1 加密布點(diǎn)區(qū)域確定

基于普查點(diǎn)位數(shù)據(jù)，分析其空間分布規(guī)律，運(yùn)用指示克里金方法預(yù)測案例區(qū)鉛污染概率分布，如圖1所示。

圖1 基于指示克里金的土壤Pb污染概率Figure 1 Probability of soil Pb pollution

已有研究通過概率制圖確定土壤污染調(diào)查加密布點(diǎn)區(qū)域[1]，但目前還沒有對不同概率區(qū)間加密點(diǎn)位效率的研究。本研究以分析比較不同概率區(qū)間加密點(diǎn)位效率為目的，劃分了8個概率區(qū)間：（0，0.05）、（0.05，0.20）、（0.20，0.35）、（0.35，0.50）、（0.50，0.65）、（0.65，0.80）、（0.80，0.95）、（0.95，1.00）。其中，（0，0.05）與（0.95，1.00）判斷為清潔或污染的概率極大，故單獨(dú)設(shè)立區(qū)間研究，其余區(qū)間均以0.15為步長設(shè)立。分別在8個概率區(qū)間內(nèi)均勻隨機(jī)選擇普查點(diǎn)位以外的點(diǎn)模擬加密布點(diǎn)，考慮到數(shù)據(jù)的有限性以及科學(xué)性，在概率介于（0.05，0.95）之間的每個區(qū)間設(shè)置加密點(diǎn)位14個，在（0，0.05）以及（0.95，1.00）兩個概率區(qū)間設(shè)置4個加密點(diǎn)位，共計92個（圖2a），加密點(diǎn)位與普查點(diǎn)位的分布見圖2b。

圖2 土壤Pb污染調(diào)查點(diǎn)位分布Figure 2 Soil Pb pollution investigation point distribution

分別對8個概率區(qū)間加密后樣點(diǎn)進(jìn)行序貫高斯模擬土壤污染范圍，與普查點(diǎn)位模擬的土壤污染范圍精度進(jìn)行比較，以精度的增減來表示效率的高低。結(jié)果表明：在非均勻布點(diǎn)的情況下，加密布點(diǎn)不一定能提高土壤污染范圍的估計精度，在不合理的區(qū)域布點(diǎn)反而會降低土壤污染范圍的估計精度。8個概率區(qū)間內(nèi)平均每加密一個點(diǎn)帶來的精度變化均不相同，如表1所示。

表1 加密區(qū)域效率Table 1 Efficiency of encryption area

2.2 加密樣本量比較

在布點(diǎn)方案的加密工作中，選擇合適的樣本量既能滿足土壤污染范圍的估計精度又能減少不必要的投入。在上節(jié)研究的基礎(chǔ)上，選擇概率范圍在0.50～0.95的區(qū)域進(jìn)行階梯式加密布點(diǎn)。共設(shè)置了7個階梯，分別加密5.7%（10個）、11.4%（20個）、17.0%（30個）、22.7%（40個）、28.4%（50個）、34.1%（60個）、39.8%（70個）點(diǎn)位。對7個加密后樣點(diǎn)分別進(jìn)行序貫高斯模擬，預(yù)測其污染面積。將7個加密樣點(diǎn)的精度與176個普查點(diǎn)位精度相比較，以精度的增減來表示效率的高低，結(jié)果見表2。結(jié)果表明，在非均勻布點(diǎn)的情況下，加密布點(diǎn)不一定能提高土壤污染范圍的估計精度。過密的點(diǎn)位不僅會造成成本的提高，還會降低土壤污染范圍的估計精度。在本案例中，加密點(diǎn)數(shù)在11.4%（20個）以內(nèi)會提高土壤污染范圍的估計精度，其中以加密5.7%（10個）點(diǎn)位效率最高。當(dāng)加密點(diǎn)數(shù)超過11.4%（20個）時，土壤污染范圍的估計精度反而降低，且呈現(xiàn)出隨點(diǎn)數(shù)增加而效率降低的趨勢。

表2 加密樣本量效率Table 2 Efficiency of encrypted sample size

3 討論

對比588個總體點(diǎn)位，除了0～0.05與0.95～1.00兩個概率區(qū)間因面積較小加密4個點(diǎn)位外，其余各區(qū)間均加密14個點(diǎn)位，即0～0.05與0.95～1.00兩個區(qū)間加密后點(diǎn)位數(shù)為180個，其余區(qū)間均為190個。根據(jù)結(jié)果來看，概率在0.50～0.95之間的區(qū)域加密布點(diǎn)時，加密布點(diǎn)對土壤污染范圍估計精度有所提升。其中概率在0.65～0.95之間的區(qū)域加密布點(diǎn)效率顯著高于0.50～0.65區(qū)域，點(diǎn)位平均效率達(dá)到3%以上，總體精度提高了40%以上。概率在0.05～0.50之間的區(qū)域加密布點(diǎn)后效率反而下降，且在這段概率區(qū)間內(nèi)，效率隨概率增大而降低。概率在0～0.05以及0.95～1.00之間的極值區(qū)加密布點(diǎn)也會造成土壤污染范圍估計精度的降低。

這是因?yàn)?，插值的方法對?shù)據(jù)密度的敏感性較高，在數(shù)據(jù)更密的區(qū)域會夸大該類型區(qū)域的預(yù)測范圍。從整個研究區(qū)來看，概率0～0.05以及0.95～1.00之間的區(qū)域已經(jīng)有極大概率被判斷為清潔區(qū)和污染區(qū)，在清潔區(qū)加密布點(diǎn)會增大清潔區(qū)的預(yù)測范圍，同理，在污染區(qū)加密布點(diǎn)會增加污染區(qū)的預(yù)測范圍。即清潔區(qū)與污染區(qū)內(nèi)固有的普查點(diǎn)位已足夠代表該區(qū)域的污染情況，進(jìn)一步加密布點(diǎn)反而會造成精度的降低與成本的浪費(fèi)。同樣，在靠近污染區(qū)的過渡區(qū)域（0.50～0.95），其污染可能性較大，但固有的普查點(diǎn)位難以達(dá)到對污染范圍較高精度的估計，在該區(qū)域加密布點(diǎn)會增加污染范圍的估計精度。而在靠近清潔區(qū)的過渡區(qū)域（0.05～0.50），其污染可能性較低，在該區(qū)域布點(diǎn)會增加清潔區(qū)的估計范圍，從而降低污染區(qū)范圍的估計精度。

從指示克里金的機(jī)制來看，概率越接近0.5的區(qū)域點(diǎn)位越稀疏，其不確定性越高。本研究為了比較不同概率區(qū)間布點(diǎn)效率，對各區(qū)間采用均勻布點(diǎn)的方法，因此對于不確定性較高的區(qū)域加密布點(diǎn)相對較少。因此也造成了概率區(qū)間在0.05～0.50之間區(qū)域，效率隨概率增加而降低的現(xiàn)象。

以往的調(diào)查經(jīng)驗(yàn)認(rèn)為，在調(diào)查區(qū)域布點(diǎn)越密，調(diào)查結(jié)果越能接近真實(shí)情況。造成這個結(jié)論的原因是目前大多數(shù)土壤調(diào)查均屬于均勻布點(diǎn)，即在調(diào)查區(qū)域采用網(wǎng)格布點(diǎn)、隨機(jī)布點(diǎn)等方法進(jìn)行布點(diǎn)方案設(shè)計。布點(diǎn)越密結(jié)果越精確正是在均勻布點(diǎn)的前提假設(shè)下才可成立。

在非均勻布點(diǎn)的情況下，在確定的清潔區(qū)或污染區(qū)加密點(diǎn)位，利用插值方法估計面積，勢必會造成對整個研究區(qū)預(yù)測的清潔面積或污染面積的增加，反而降低了土壤污染范圍的估計精度。因此不應(yīng)在大概率確定污染（清潔）的區(qū)域加密布點(diǎn)，而應(yīng)該在疑似污染但不能確定的區(qū)域加密布點(diǎn)。并且，在疑似污染區(qū)加密布點(diǎn)的結(jié)果精度隨加密點(diǎn)數(shù)增加呈先增加后降低的趨勢。這是由于疑似污染區(qū)隨著加密點(diǎn)位的增多，判斷為污染區(qū)或清潔區(qū)的概率逐漸增大，最終變成確定的污染區(qū)或清潔區(qū)，若繼續(xù)加密布點(diǎn)就會夸大該區(qū)域的面積，從而降低精度。這與往常采用均勻布點(diǎn)方法的情況有所不同。

簡單來說，假設(shè)研究區(qū)為A，加密區(qū)為B，即有B?A，若B為大概率判斷為清潔區(qū)或污染區(qū)（已知區(qū)域），用加密布點(diǎn)的數(shù)據(jù)來估計A的情況勢必會造成誤判（若只用加密布點(diǎn)的數(shù)據(jù)僅對B進(jìn)行估計是可以有效提升精度的）；若B為不確定區(qū)域（未知區(qū)域），用加密布點(diǎn)的數(shù)據(jù)對A的情況進(jìn)行估計是有可能提高估計精度的（是否可以提高還需看加密的點(diǎn)數(shù)是否合適）。

4 結(jié)論與局限

4.1 研究結(jié)論

本研究優(yōu)化了基于指示克里金的確定土壤污染調(diào)查加密布點(diǎn)區(qū)域的方法，即利用已有調(diào)查數(shù)據(jù)，采用類似抽樣調(diào)查的方法先進(jìn)行模擬加密布點(diǎn)，確定加密布點(diǎn)最優(yōu)布點(diǎn)區(qū)域以及加密程度，在此基礎(chǔ)上再進(jìn)行正式加密。該方法的意義不僅在于節(jié)省了調(diào)查成本，更在于提高了結(jié)果精度，避免了無目的加密造成的預(yù)測偏差。

在該方法的基礎(chǔ)上研究了不同概率區(qū)間內(nèi)的加密布點(diǎn)效率與加密樣本量的效率。結(jié)果表明：在非均勻布點(diǎn)的前提下，加密布點(diǎn)不一定能提高土壤污染范圍的估計精度。在確定的清潔區(qū)或污染區(qū)加密布點(diǎn)以及在疑似污染區(qū)加密布點(diǎn)過密，都會降低整個研究區(qū)土壤污染范圍的估計精度，在疑似污染區(qū)加密適當(dāng)?shù)狞c(diǎn)位數(shù)才有助于估計精度的提高。案例研究顯示：

（1）概率在0.50～0.95之間的區(qū)域，加密布點(diǎn)對土壤污染范圍估計精度帶來了提升；概率在0.05～0.50之間的區(qū)域加密布點(diǎn)效率反而下降；概率在0～0.05以及0.95～1.00之間的極值區(qū)加密布點(diǎn)也會造成土壤污染范圍估計精度的降低。

（2）加密點(diǎn)數(shù)在11.4%（20個）以內(nèi)會提高土壤污染范圍的估計精度，而當(dāng)加密點(diǎn)數(shù)超過11.4%（20個）時，土壤污染范圍的估計精度反而降低，且呈現(xiàn)出隨點(diǎn)數(shù)增加而效率降低的趨勢。

4.2 局限性

（1）由于土壤污染空間分布較為復(fù)雜，不同尺度下土壤污染分布規(guī)律可能會有所不同，由本案例得出的結(jié)論是否適用于更大尺度的區(qū)域，還需要做進(jìn)一步的研究。

（2）本案例區(qū)屬于高原區(qū)土壤鉛污染較重區(qū)域，由本案例得出的結(jié)論是否適用于不同地形區(qū)域、不同污染物、不同污染程度區(qū)域，還有待進(jìn)一步的研究。

（3）本研究所提出的優(yōu)化基于指示克里金的確定土壤污染調(diào)查加密布點(diǎn)區(qū)域的方法僅適用于達(dá)到一定點(diǎn)位數(shù)量的大中型規(guī)模調(diào)查，小規(guī)模調(diào)查點(diǎn)位較少不足以支撐模擬研究。