基于序貫高斯條件模擬的土壤重金屬含量預(yù)測(cè)與不確定性評(píng)價(jià)
——以宜興市土壤Hg為例

徐辰星 濮勵(lì)杰，2? 朱 明，2 徐彩瑤 張 濛 許 艷

（1 南京大學(xué)地理與海洋科學(xué)學(xué)院，南京 210023）

（2 國(guó)土資源部海岸帶開(kāi)發(fā)與保護(hù)重點(diǎn)實(shí)驗(yàn)室，南京 210023）

（3 蘇州科技大學(xué)環(huán)境科學(xué)與工程學(xué)院，江蘇蘇州 215009）

土壤屬性制圖的主要目的是表達(dá)土壤異質(zhì)的空間信息，但由于空間預(yù)測(cè)模型的限制，由采樣點(diǎn)位置的觀測(cè)值來(lái)預(yù)測(cè)土壤性質(zhì)空間分布的過(guò)程中包含的不確定性，將通過(guò)一系列的傳遞并對(duì)最終結(jié)果產(chǎn)生深刻影響［1］。克里格估值只能提供簡(jiǎn)化變異的空間分布模式，在極值和空間連續(xù)性模式敏感的應(yīng)用中存在很大的局限性［2-4］，比如典型的土壤污染數(shù)據(jù)。隨機(jī)模擬方法能克服“平滑效應(yīng)”，它將數(shù)據(jù)作為一個(gè)整體，最大可能地接近真實(shí)的空間分布，因此近年來(lái)空間隨機(jī)模擬方法的應(yīng)用逐漸成為地統(tǒng)計(jì)學(xué)研究的主要趨勢(shì)之一［5-8］。目前建模土壤性質(zhì)常用的隨機(jī)模擬方法包括序貫高斯條件模擬、序貫高斯協(xié)同模擬、序貫高斯指示模擬和序貫高斯指示協(xié)同模擬［8-10］。

據(jù)中國(guó)2014年發(fā)布的《全國(guó)土壤污染狀況調(diào)查公報(bào)》顯示：全國(guó)不同區(qū)域的土壤樣品中總Hg含量為0.017～3.32 mg kg-1，表層0～20 cm的土壤Hg含量最高，點(diǎn)位超標(biāo)率為1.6%［11］。長(zhǎng)江三角洲重金屬污染嚴(yán)重，Hg是主要的重金屬污染元素。Hg進(jìn)入土壤后，95%以上被固定在土壤表層［12］，且受人為活動(dòng)影響較大。近年來(lái)蘇南城市化建設(shè)與經(jīng)濟(jì)發(fā)展如雨后春筍，土壤污染的報(bào)道較為常見(jiàn)，本文以宜興市土壤重金屬地質(zhì)調(diào)查數(shù)據(jù)為例，利用序貫高斯條件模擬方法描述土壤Hg的空間分布特征，并利用高斯序貫指示模擬進(jìn)行土壤Hg超過(guò)設(shè)定閾值的概率分布評(píng)價(jià)，探索該區(qū)域土壤重金屬含量空間預(yù)測(cè)的不確定性評(píng)價(jià)方法，對(duì)控制區(qū)域土壤重金屬風(fēng)險(xiǎn)具有重要的理論和現(xiàn)實(shí)意義。

1 材料與方法

1.1 研究區(qū)概況

宜興市地處江蘇省南部（31°07′N(xiāo)～31°37′N(xiāo)，119°31′E～120°03′E），總面積1 996 km2（含太湖水域面積242.3 km2）。宜興市四季分明，全年溫暖濕潤(rùn)，年平均氣溫15.7℃，1月平均氣溫2.8℃，7月平均氣溫28.3℃。年平均無(wú)霜期240多天。年平均雨日136.6 d，年平均降水量約1 177 mm，春夏雨水集中［13］。宜興市自改革開(kāi)放以來(lái)工業(yè)發(fā)展迅速，地形因素（地勢(shì)南高北低，南部為丘陵山區(qū)，北部為平原區(qū)，東部為太湖瀆區(qū)，西部為低洼圩區(qū)）也使得土壤重金屬分布的異質(zhì)性較強(qiáng)［14］。研究區(qū)的土壤以黃棕壤、水稻土、潮土為主。

1.2 樣品采集與分析

數(shù)據(jù)來(lái)源于無(wú)錫市耕地質(zhì)量生態(tài)地球化學(xué)調(diào)查，野外采樣于2009年9月1日—15日完成。樣點(diǎn)布設(shè)采用網(wǎng)格法，網(wǎng)格寬度為1 km，實(shí)際野外采樣過(guò)程中，若樣點(diǎn)落在建筑物或交通道路上，則選擇采樣點(diǎn)周?chē)邆洳杉瘲l件的地點(diǎn)采自表層0～20 cm，按照梅花采樣法，采集半徑20 m內(nèi)的4～8個(gè)點(diǎn)混合樣，裝入樣品袋帶回實(shí)驗(yàn)室待測(cè)。各樣點(diǎn)均用GPS記錄其經(jīng)緯度，并詳細(xì)記錄周?chē)沫h(huán)境信息。共采集了1 157個(gè)土壤表層樣品（圖1），土樣風(fēng)干后用瑪瑙研缽研碎過(guò)20目尼龍篩，裝袋備用，土壤Hg含量依照GB/T 22105《土壤質(zhì)量總汞、總砷、總鉛的測(cè)定 原子熒光法》，采用電感耦合等離子體原子熒光光譜法（ICP-MS）測(cè)定。

圖1 研究區(qū)土壤采樣點(diǎn)分布圖Fig. 1 Sampling sites distribution of study area

1.3 數(shù)據(jù)處理

經(jīng)實(shí)驗(yàn)室化學(xué)分析后，有兩個(gè)樣點(diǎn)的值分別為1.34和0.825 mg kg-1，考慮到重金屬元素分布變異程度較高的固有屬性，參照拉依達(dá)準(zhǔn)則將誤差標(biāo)準(zhǔn)擴(kuò)充至平均值與10倍標(biāo)準(zhǔn)差的和，因此將這兩個(gè)數(shù)據(jù)視為實(shí)驗(yàn)誤差剔除，實(shí)際用于本研究的樣本數(shù)量為1 155個(gè)。原始數(shù)據(jù)不符合正態(tài)分布，偏度和峰度分別為1.795和9.495。研究區(qū)內(nèi)有25個(gè)樣品高于國(guó)家土壤環(huán)境質(zhì)量二級(jí)標(biāo)準(zhǔn)（GB15618.1995）0.3 mg kg-1，占樣本總數(shù)的2.16%，高于全國(guó)平均水平［11］。本文利用GS+9.0進(jìn)行半方差函數(shù)模型的擬合，運(yùn)用ArcGIS10.1的高斯地統(tǒng)計(jì)模擬模塊實(shí)現(xiàn)序貫高斯條件模擬。為了研究地統(tǒng)計(jì)模擬的特點(diǎn)和可靠性，將原始數(shù)據(jù)分為模擬數(shù)據(jù)集和驗(yàn)證數(shù)據(jù)集，以4 km×4 km網(wǎng)格對(duì)原采樣點(diǎn)進(jìn)行重采樣，得到由309個(gè)樣本組成的模擬數(shù)據(jù)集，其余846個(gè)樣本則作為驗(yàn)證數(shù)據(jù)集用于模擬精度的驗(yàn)證。

表1 土壤Hg含量描述性統(tǒng)計(jì)Table 1 Descriptive statistics for soil Hg content

本研究使用Box-Cox法對(duì)模擬數(shù)據(jù)集進(jìn)行轉(zhuǎn)換嘗試變化參數(shù)，然后檢驗(yàn)變換后的正態(tài)性是否足夠理想，最終設(shè)置參數(shù)為0.202，轉(zhuǎn)換后的偏度和峰度分別為0.0009和3.62，通過(guò)SPSS中的單樣本K-S檢驗(yàn)，說(shuō)明數(shù)據(jù)符合正態(tài)分布，可用于進(jìn)一步分析。對(duì)模擬數(shù)據(jù)集中的樣點(diǎn)做半方差云圖和趨勢(shì)分析，未發(fā)現(xiàn)明顯的各向異性和趨勢(shì)。

1.4 研究方法

（1）半方差模型和序貫高斯條件模擬。本文采用半方差模型描述連續(xù)的土壤屬性的空間自相關(guān)結(jié)構(gòu)，揭示可能影響數(shù)據(jù)分布的因素，并用序貫高斯條件模擬（Sequential Gaussian Conditional Simulation, SGCS）定量刻畫(huà)土壤屬性的變異性和不確定性，通過(guò)隨機(jī)變量的系列結(jié)果來(lái)統(tǒng)計(jì)不確定性。

（2）精度評(píng)價(jià)。為了測(cè)試半方差模型的擬合效果以及驗(yàn)證高斯地統(tǒng)計(jì)模擬的精確度和可靠性，需利用一種客觀的檢驗(yàn)方法進(jìn)行評(píng)價(jià)。對(duì)序貫高斯條件模擬的精度，采用驗(yàn)證數(shù)據(jù)集中樣點(diǎn)的Hg含量的平均預(yù)測(cè)誤差（Mean Prediction Error, MPE）和均方根預(yù)測(cè)誤差（Root Mean Square Prediction Error, RMSPE）來(lái)評(píng)價(jià)［16］。

（3）不確定性評(píng)價(jià)。本文采用高斯序貫指示模擬探討土壤Hg的不確定性［10，17］。單一空間位置x′上Hg含量的不確定性可以用該點(diǎn)Hg含量y′大于某個(gè)給定閾值 yt的概率來(lái)定量表示，空間不確定性（spatial uncertainty）可以用于評(píng)價(jià)基于概率Pexceed描繪重金屬空間分布的可信度［20］。

2 結(jié)果與討論

2.1 土壤Hg含量空間變異結(jié)構(gòu)特征

研究區(qū)土壤表層Hg含量的平均值為0.127 mg kg-1，變異系數(shù)達(dá)52.4%。半方差模型參數(shù)見(jiàn)表2，塊金系數(shù)達(dá)0.504，屬于中等強(qiáng)度變異，說(shuō)明土壤Hg受到土壤母質(zhì)、地形等結(jié)構(gòu)因素和耕作、土壤改良等隨機(jī)因素的共同影響。本文依據(jù)決定系數(shù)的大小，選擇了擬合效果（R2=0.930）最佳的球狀模型，并用此模型的參數(shù)進(jìn)行克里格插值。

克里格插值產(chǎn)生的交叉驗(yàn)證報(bào)告，均方根預(yù)測(cè)誤差，標(biāo)準(zhǔn)平均值，平均標(biāo)準(zhǔn)誤差已達(dá)到最小，并且均方根預(yù)測(cè)誤差與平均標(biāo)準(zhǔn)誤差非常接近，符合一般情況下模型最優(yōu)的準(zhǔn)則。本研究將序貫高斯條件模擬方法的實(shí)現(xiàn)次數(shù)設(shè)置為100，E-type模擬的土壤Hg空間分布與簡(jiǎn)單克里格估值呈現(xiàn)了相似的趨勢(shì)（圖2），高值集中在研究區(qū)的東北部太湖瀆區(qū)和西部的洼地圩區(qū)的。

表2 Hg含量空間變異結(jié)構(gòu)特征Table 2 Structural characteristics of spatial variation of Hg content

2.2 空間預(yù)測(cè)精度

將序貫高斯條件模擬得到E-type型，與抽取的第1、25、50、75、100次實(shí)現(xiàn)，提取柵格上的數(shù)值到驗(yàn)證數(shù)據(jù)集中的846個(gè)點(diǎn)上，對(duì)比模擬值與真實(shí)測(cè)量值。隨機(jī)抽取的5個(gè)單次實(shí)現(xiàn)的平均預(yù)測(cè)誤差與均方根預(yù)測(cè)誤差均較簡(jiǎn)單克里格插值大（相差在0.01以?xún)?nèi)），E-type型與簡(jiǎn)單克里格的平均預(yù)測(cè)誤差和均方根誤差相差在0.000 1以?xún)?nèi)，可見(jiàn)兩種不同原理的方法在研究區(qū)尺度上運(yùn)用的精度可以視作無(wú)差別。將閾值設(shè)定為當(dāng)?shù)赝寥繦g的背景值0.2 mg kg-1和用于不確定性研究的0.15 mg kg-1的分類(lèi)錯(cuò)誤率范圍分別為11.70%～12.88%和27.19%～30.14%。雖然5個(gè)單次實(shí)現(xiàn)的分類(lèi)錯(cuò)誤率均高于E-type，但是簡(jiǎn)單克里格估值與E-type型的分類(lèi)錯(cuò)誤率不存在顯著差異，隨機(jī)模擬與傳統(tǒng)的地統(tǒng)計(jì)方法在精度上差異不顯著。

簡(jiǎn)單克里格得到的空間分布更為平滑，E-type型會(huì)隨著隨機(jī)模擬指定次數(shù)的增加而更趨向于克里格估值（圖2a、圖2b）。平滑效應(yīng)造成估值的數(shù)值值域范圍較小［21］，克里格估值的取值范圍為0.050～0.265 mg kg-1，E-type型模擬取值范圍為0.048～0.276 mg kg-1。模擬數(shù)據(jù)集的峰度稍大，兩種方法均會(huì)造成柵格數(shù)據(jù)向中值靠攏。在這種情形下，克里格插值的平滑效應(yīng)尤為明顯，數(shù)值集聚在0.121～0.18 mg kg-1。抽取的5個(gè)單次實(shí)現(xiàn)的模擬值均呈現(xiàn)正偏態(tài)，對(duì)柵格上的數(shù)值統(tǒng)計(jì)：0.091～0.12 mg kg-1的頻數(shù)最大，0.061～0.09 mg kg-1和0.121～0.15 mg kg-1的頻數(shù)次之，與輸入數(shù)據(jù)的描述性統(tǒng)計(jì)一致。標(biāo)準(zhǔn)差在距離采樣點(diǎn)較遠(yuǎn)的地方以及研究區(qū)邊界有較大的誤差方差，簡(jiǎn)單克里格得到的標(biāo)準(zhǔn)差在采樣點(diǎn)半徑500 m內(nèi)較小，E-type型標(biāo)準(zhǔn)差（圖2c、圖2d）只在采樣點(diǎn)局部較小，主要是因?yàn)榭死锔穹讲顑H依賴(lài)于采樣點(diǎn)的分布情況。

序貫高斯條件模擬借助Monto-Carlo方法，依據(jù)原始數(shù)據(jù)遵循的概率密度函數(shù)求解更具優(yōu)勢(shì)，模擬值的方差在任何時(shí)候均等于原始數(shù)據(jù)所遵循的高斯分布的方差，避免了平滑效應(yīng)。此外，序貫高斯模擬條件將局部變異性重新添加到了其生成的表面中，克服了克里格方差與實(shí)際樣本值無(wú)關(guān)的缺點(diǎn)［21］。克里格估值在每個(gè)位置得到的是最優(yōu)估值，不存在任何可能的變化，誤差表征的不確定性表現(xiàn)在不同的預(yù)測(cè)點(diǎn)位置之間。序貫高斯條件模擬的隨機(jī)性服從條件，每次模擬結(jié)果再現(xiàn)輸入數(shù)據(jù)的統(tǒng)計(jì)特征，即均值、方差、協(xié)方差、變異函數(shù)均與原始數(shù)據(jù)的均值、方差、協(xié)方差、變異函數(shù)保持一致。若干次模擬實(shí)現(xiàn)就在每個(gè)位置上產(chǎn)生了多個(gè)模擬結(jié)果，具備了以概率論為基礎(chǔ)的不確定性評(píng)價(jià)，可以對(duì)該位置空間變量可能的取值區(qū)間做數(shù)理統(tǒng)計(jì)。

2.3 空間預(yù)測(cè)的不確定性評(píng)價(jià)

本文以評(píng)價(jià)標(biāo)準(zhǔn)0.15 mg kg-1作為土壤Hg含量的閾值為例進(jìn)行序貫指示模擬，模擬次數(shù)設(shè)為1 000次，利用局部不確定性評(píng)價(jià)［17］計(jì)算各個(gè)點(diǎn)位超過(guò)該臨界閾值的概率，對(duì)這些區(qū)域做超閾值的聯(lián)合分布概率，可以評(píng)價(jià)“污染區(qū)”的風(fēng)險(xiǎn)程度。

本文以1 000次高斯序貫指示模擬為例探討空間多個(gè)位置的不確定性，最大的標(biāo)準(zhǔn)差為0.053 8，落在了未采土壤樣本的水域，其次較大的標(biāo)準(zhǔn)差0.045落在臨界概率為55%～65%的區(qū)域內(nèi)。標(biāo)準(zhǔn)差與臨界概率兩者比值的最大值可用作表征相對(duì)最大誤差，最大比值為8.18%表明結(jié)果可靠。將指示柵格的大小設(shè)定為100 m2，用超閾值概率［17］表征不同置信度下的土壤Hg超閾值區(qū)域。按95%、90%、85%、75%的閾值概率分別有206個(gè)、307個(gè)、511個(gè)和1 838個(gè)柵格，再對(duì)這些區(qū)域進(jìn)行超閾值的聯(lián)合評(píng)價(jià)，由表5所示，單點(diǎn)閾值概率達(dá)到95%、90%、85%區(qū)域的聯(lián)合概率分別為24.76%、16.61%和10.00%，多點(diǎn)聯(lián)合概率較單點(diǎn)臨界概率的評(píng)價(jià)標(biāo)準(zhǔn)更加嚴(yán)格。當(dāng)單點(diǎn)的臨界概率達(dá)到75%時(shí)，采用聯(lián)合概率進(jìn)行計(jì)算只有2.77%，說(shuō)明盡管臨界概率比較高，但置信度不足以用作劃定依據(jù)。

圖2 Hg空間預(yù)測(cè)圖和標(biāo)準(zhǔn)差圖Fig. 2 Spatial prediction map and standard deviation map of soil Hg content

圖3 序貫高斯條件模擬的單次實(shí)現(xiàn)Fig. 3 Values distributions of SGCS

表層（0～20 cm）土壤Hg含量空間表達(dá)波動(dòng)性大的區(qū)域，主要集中在該市距離太湖約3 km的新莊鎮(zhèn)及和橋鎮(zhèn)，同時(shí)該區(qū)域也是Hg含量較高的一帶。結(jié)合土地利用類(lèi)型和地圖，新莊鎮(zhèn)有五十多家燈具企業(yè)，江蘇省首座含汞燈管處理裝置在和橋鎮(zhèn)投入運(yùn)行，用于粉碎和酸洗回收金屬汞和硫酸汞。固體廢棄物如燈管的露天堆放及管理不善，易造成包膜破碎和污染物質(zhì)向土壤滲透。污染物可通過(guò)大氣沉降、徑流等途徑［22］，造成汞化合物對(duì)農(nóng)田、水源的嚴(yán)重污染。經(jīng)核實(shí)地圖變更信息，宜興市曾有六十余家涉鉛、鎘、汞、鉻等重金屬的化工企業(yè)，上百家涉汞燈具企業(yè)，數(shù)十家鉛蓄電池企業(yè)。以及紫砂壺、玻璃、陶瓷、搪瓷、橡膠使用到的化工無(wú)機(jī)染料，農(nóng)業(yè)有機(jī)肥的大量累積施用，均有可能造成土壤Hg含量的累積。土壤母質(zhì)是土壤汞的最基本來(lái)源，但外界干擾是引起土壤Hg含量和空間波動(dòng)性升高的重要因素［12］。結(jié)合高程DEM，土地利用方式和坡度是影響汞流失的重要原因，南部丘陵地區(qū)多個(gè)山莊和風(fēng)景區(qū)的土壤Hg含量低，且空間表達(dá)的波動(dòng)性很小。

圖4 臨界概率Pc條件下獲取的土壤Hg含量＞0.15 mg kg-1的區(qū)域分布Fig. 4 Spatial distribution of Hg content＞0.15 mg kg-1 under the conditional of critical probability Pc

3 結(jié) 論

兩種方法在本文的應(yīng)用中精度無(wú)異；克里格估值具有明顯的“平滑效應(yīng)”，體現(xiàn)在空間分布的連續(xù)、估值值域范圍的收束與數(shù)值分布向擬合面的集中，而序貫高斯條件模擬的數(shù)值服從概率密度函數(shù)，尊重了原始數(shù)據(jù)的分布，較好地保留了數(shù)據(jù)的變異特征；高斯序貫指示模擬可以評(píng)價(jià)空間任意位置上Hg含量大于該閾值的概率，衍生的不確定性研究結(jié)果有助于定量評(píng)估環(huán)境風(fēng)險(xiǎn)或健康風(fēng)險(xiǎn)、檢測(cè)各種資源環(huán)境應(yīng)用模型的穩(wěn)健性以及誤差傳播等問(wèn)題。