基于序貫高斯條件模擬的土壤重金屬含量預測與不確定性評價
——以宜興市土壤Hg為例

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

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

（2 國土資源部海岸帶開發(fā)與保護重點實驗室，南京 210023）

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

1.2 樣品采集與分析

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

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

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

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

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

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

1.4 研究方法

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

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

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

2 結果與討論

2.1 土壤Hg含量空間變異結構特征

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

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

表2 Hg含量空間變異結構特征Table 2 Structural characteristics of spatial variation of Hg content

2.2 空間預測精度

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

簡單克里格得到的空間分布更為平滑，E-type型會隨著隨機模擬指定次數(shù)的增加而更趨向于克里格估值（圖2a、圖2b）。平滑效應造成估值的數(shù)值值域范圍較?。?1］，克里格估值的取值范圍為0.050～0.265 mg kg-1，E-type型模擬取值范圍為0.048～0.276 mg kg-1。模擬數(shù)據(jù)集的峰度稍大，兩種方法均會造成柵格數(shù)據(jù)向中值靠攏。在這種情形下，克里格插值的平滑效應尤為明顯，數(shù)值集聚在0.121～0.18 mg kg-1。抽取的5個單次實現(xiàn)的模擬值均呈現(xiàn)正偏態(tài)，對柵格上的數(shù)值統(tǒng)計：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)計一致。標準差在距離采樣點較遠的地方以及研究區(qū)邊界有較大的誤差方差，簡單克里格得到的標準差在采樣點半徑500 m內(nèi)較小，E-type型標準差（圖2c、圖2d）只在采樣點局部較小，主要是因為克里格方差僅依賴于采樣點的分布情況。

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

2.3 空間預測的不確定性評價

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

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

圖2 Hg空間預測圖和標準差圖Fig. 2 Spatial prediction map and standard deviation map of soil Hg content

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

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

圖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 結 論

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