基于序贯高斯条件模拟的土壤重金属含量预测与不确定性评价——以宜兴市土壤Hg为例

徐辰星 濮励杰，2† 朱 明，2 徐彩瑶 张 濛 许 艳

（1 南京大学地理与海洋科学学院，南京 210023）

（2 国土资源部海岸带开发与保护重点实验室，南京 210023）

（3 苏州科技大学环境科学与工程学院，江苏苏州 215009）

土壤属性制图的主要目的是表达土壤异质的空间信息，但由于空间预测模型的限制，由采样点位置的观测值来预测土壤性质空间分布的过程中包含的不确定性，将通过一系列的传递并对最终结果产生深刻影响［1］。克里格估值只能提供简化变异的空间分布模式，在极值和空间连续性模式敏感的应用中存在很大的局限性［2-4］，比如典型的土壤污染数据。随机模拟方法能克服“平滑效应”，它将数据作为一个整体，最大可能地接近真实的空间分布，因此近年来空间随机模拟方法的应用逐渐成为地统计学研究的主要趋势之一［5-8］。目前建模土壤性质常用的随机模拟方法包括序贯高斯条件模拟、序贯高斯协同模拟、序贯高斯指示模拟和序贯高斯指示协同模拟［8-10］。

据中国2014年发布的《全国土壤污染状况调查公报》显示：全国不同区域的土壤样品中总Hg含量为0.017～3.32 mg kg-1，表层0～20 cm的土壤Hg含量最高，点位超标率为1.6%［11］。长江三角洲重金属污染严重，Hg是主要的重金属污染元素。Hg进入土壤后，95%以上被固定在土壤表层［12］，且受人为活动影响较大。近年来苏南城市化建设与经济发展如雨后春笋，土壤污染的报道较为常见，本文以宜兴市土壤重金属地质调查数据为例，利用序贯高斯条件模拟方法描述土壤Hg的空间分布特征，并利用高斯序贯指示模拟进行土壤Hg超过设定阈值的概率分布评价，探索该区域土壤重金属含量空间预测的不确定性评价方法，对控制区域土壤重金属风险具有重要的理论和现实意义。

1 材料与方法

1.1 研究区概况

宜兴市地处江苏省南部（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］。宜兴市自改革开放以来工业发展迅速，地形因素（地势南高北低，南部为丘陵山区，北部为平原区，东部为太湖渎区，西部为低洼圩区）也使得土壤重金属分布的异质性较强［14］。研究区的土壤以黄棕壤、水稻土、潮土为主。

1.2 样品采集与分析

数据来源于无锡市耕地质量生态地球化学调查，野外采样于2009年9月1日—15日完成。样点布设采用网格法，网格宽度为1 km，实际野外采样过程中，若样点落在建筑物或交通道路上，则选择采样点周围具备采集条件的地点采自表层0～20 cm，按照梅花采样法，采集半径20 m内的4～8个点混合样，装入样品袋带回实验室待测。各样点均用GPS记录其经纬度，并详细记录周围的环境信息。共采集了1 157个土壤表层样品（图1），土样风干后用玛瑙研钵研碎过20目尼龙筛，装袋备用，土壤Hg含量依照GB/T 22105《土壤质量总汞、总砷、总铅的测定 原子荧光法》，采用电感耦合等离子体原子荧光光谱法（ICP-MS）测定。

图1 研究区土壤采样点分布图Fig. 1 Sampling sites distribution of study area

1.3 数据处理

经实验室化学分析后，有两个样点的值分别为1.34和0.825 mg kg-1，考虑到重金属元素分布变异程度较高的固有属性，参照拉依达准则将误差标准扩充至平均值与10倍标准差的和，因此将这两个数据视为实验误差剔除，实际用于本研究的样本数量为1 155个。原始数据不符合正态分布，偏度和峰度分别为1.795和9.495。研究区内有25个样品高于国家土壤环境质量二级标准（GB15618.1995）0.3 mg kg-1，占样本总数的2.16%，高于全国平均水平［11］。本文利用GS+9.0进行半方差函数模型的拟合，运用ArcGIS10.1的高斯地统计模拟模块实现序贯高斯条件模拟。为了研究地统计模拟的特点和可靠性，将原始数据分为模拟数据集和验证数据集，以4 km×4 km网格对原采样点进行重采样，得到由309个样本组成的模拟数据集，其余846个样本则作为验证数据集用于模拟精度的验证。

表1 土壤Hg含量描述性统计Table 1 Descriptive statistics for soil Hg content

本研究使用Box-Cox法对模拟数据集进行转换尝试变化参数，然后检验变换后的正态性是否足够理想，最终设置参数为0.202，转换后的偏度和峰度分别为0.0009和3.62，通过SPSS中的单样本K-S检验，说明数据符合正态分布，可用于进一步分析。对模拟数据集中的样点做半方差云图和趋势分析，未发现明显的各向异性和趋势。

1.4 研究方法

（1）半方差模型和序贯高斯条件模拟。本文采用半方差模型描述连续的土壤属性的空间自相关结构，揭示可能影响数据分布的因素，并用序贯高斯条件模拟（Sequential Gaussian Conditional Simulation, SGCS）定量刻画土壤属性的变异性和不确定性，通过随机变量的系列结果来统计不确定性。

（2）精度评价。为了测试半方差模型的拟合效果以及验证高斯地统计模拟的精确度和可靠性，需利用一种客观的检验方法进行评价。对序贯高斯条件模拟的精度，采用验证数据集中样点的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含量空间变异结构特征

研究区土壤表层Hg含量的平均值为0.127 mg kg-1，变异系数达52.4%。半方差模型参数见表2，块金系数达0.504，属于中等强度变异，说明土壤Hg受到土壤母质、地形等结构因素和耕作、土壤改良等随机因素的共同影响。本文依据决定系数的大小，选择了拟合效果（R2=0.930）最佳的球状模型，并用此模型的参数进行克里格插值。

克里格插值产生的交叉验证报告，均方根预测误差，标准平均值，平均标准误差已达到最小，并且均方根预测误差与平均标准误差非常接近，符合一般情况下模型最优的准则。本研究将序贯高斯条件模拟方法的实现次数设置为100，E-type模拟的土壤Hg空间分布与简单克里格估值呈现了相似的趋势（图2），高值集中在研究区的东北部太湖渎区和西部的洼地圩区的。

表2 Hg含量空间变异结构特征Table 2 Structural characteristics of spatial variation of Hg content

2.2 空间预测精度

将序贯高斯条件模拟得到E-type型，与抽取的第1、25、50、75、100次实现，提取栅格上的数值到验证数据集中的846个点上，对比模拟值与真实测量值。随机抽取的5个单次实现的平均预测误差与均方根预测误差均较简单克里格插值大（相差在0.01以内），E-type型与简单克里格的平均预测误差和均方根误差相差在0.000 1以内，可见两种不同原理的方法在研究区尺度上运用的精度可以视作无差别。将阈值设定为当地土壤Hg的背景值0.2 mg kg-1和用于不确定性研究的0.15 mg kg-1的分类错误率范围分别为11.70%～12.88%和27.19%～30.14%。虽然5个单次实现的分类错误率均高于E-type，但是简单克里格估值与E-type型的分类错误率不存在显著差异，随机模拟与传统的地统计方法在精度上差异不显著。

简单克里格得到的空间分布更为平滑，E-type型会随着随机模拟指定次数的增加而更趋向于克里格估值（图2a、图2b）。平滑效应造成估值的数值值域范围较小［21］，克里格估值的取值范围为0.050～0.265 mg kg-1，E-type型模拟取值范围为0.048～0.276 mg kg-1。模拟数据集的峰度稍大，两种方法均会造成栅格数据向中值靠拢。在这种情形下，克里格插值的平滑效应尤为明显，数值集聚在0.121～0.18 mg kg-1。抽取的5个单次实现的模拟值均呈现正偏态，对栅格上的数值统计：0.091～0.12 mg kg-1的频数最大，0.061～0.09 mg kg-1和0.121～0.15 mg kg-1的频数次之，与输入数据的描述性统计一致。标准差在距离采样点较远的地方以及研究区边界有较大的误差方差，简单克里格得到的标准差在采样点半径500 m内较小，E-type型标准差（图2c、图2d）只在采样点局部较小，主要是因为克里格方差仅依赖于采样点的分布情况。

序贯高斯条件模拟借助Monto-Carlo方法，依据原始数据遵循的概率密度函数求解更具优势，模拟值的方差在任何时候均等于原始数据所遵循的高斯分布的方差，避免了平滑效应。此外，序贯高斯模拟条件将局部变异性重新添加到了其生成的表面中，克服了克里格方差与实际样本值无关的缺点［21］。克里格估值在每个位置得到的是最优估值，不存在任何可能的变化，误差表征的不确定性表现在不同的预测点位置之间。序贯高斯条件模拟的随机性服从条件，每次模拟结果再现输入数据的统计特征，即均值、方差、协方差、变异函数均与原始数据的均值、方差、协方差、变异函数保持一致。若干次模拟实现就在每个位置上产生了多个模拟结果，具备了以概率论为基础的不确定性评价，可以对该位置空间变量可能的取值区间做数理统计。

2.3 空间预测的不确定性评价

本文以评价标准0.15 mg kg-1作为土壤Hg含量的阈值为例进行序贯指示模拟，模拟次数设为1 000次，利用局部不确定性评价［17］计算各个点位超过该临界阈值的概率，对这些区域做超阈值的联合分布概率，可以评价“污染区”的风险程度。

本文以1 000次高斯序贯指示模拟为例探讨空间多个位置的不确定性，最大的标准差为0.053 8，落在了未采土壤样本的水域，其次较大的标准差0.045落在临界概率为55%～65%的区域内。标准差与临界概率两者比值的最大值可用作表征相对最大误差，最大比值为8.18%表明结果可靠。将指示栅格的大小设定为100 m2，用超阈值概率［17］表征不同置信度下的土壤Hg超阈值区域。按95%、90%、85%、75%的阈值概率分别有206个、307个、511个和1 838个栅格，再对这些区域进行超阈值的联合评价，由表5所示，单点阈值概率达到95%、90%、85%区域的联合概率分别为24.76%、16.61%和10.00%，多点联合概率较单点临界概率的评价标准更加严格。当单点的临界概率达到75%时，采用联合概率进行计算只有2.77%，说明尽管临界概率比较高，但置信度不足以用作划定依据。

图2 Hg空间预测图和标准差图Fig. 2 Spatial prediction map and standard deviation map of soil Hg content

图3 序贯高斯条件模拟的单次实现Fig. 3 Values distributions of SGCS

表层（0～20 cm）土壤Hg含量空间表达波动性大的区域，主要集中在该市距离太湖约3 km的新庄镇及和桥镇，同时该区域也是Hg含量较高的一带。结合土地利用类型和地图，新庄镇有五十多家灯具企业，江苏省首座含汞灯管处理装置在和桥镇投入运行，用于粉碎和酸洗回收金属汞和硫酸汞。固体废弃物如灯管的露天堆放及管理不善，易造成包膜破碎和污染物质向土壤渗透。污染物可通过大气沉降、径流等途径［22］，造成汞化合物对农田、水源的严重污染。经核实地图变更信息，宜兴市曾有六十余家涉铅、镉、汞、铬等重金属的化工企业，上百家涉汞灯具企业，数十家铅蓄电池企业。以及紫砂壶、玻璃、陶瓷、搪瓷、橡胶使用到的化工无机染料，农业有机肥的大量累积施用，均有可能造成土壤Hg含量的累积。土壤母质是土壤汞的最基本来源，但外界干扰是引起土壤Hg含量和空间波动性升高的重要因素［12］。结合高程DEM，土地利用方式和坡度是影响汞流失的重要原因，南部丘陵地区多个山庄和风景区的土壤Hg含量低，且空间表达的波动性很小。

图4 临界概率Pc条件下获取的土壤Hg含量＞0.15 mg kg-1的区域分布Fig. 4 Spatial distribution of Hg content＞0.15 mg kg-1 under the conditional of critical probability Pc

3 结 论

两种方法在本文的应用中精度无异；克里格估值具有明显的“平滑效应”，体现在空间分布的连续、估值值域范围的收束与数值分布向拟合面的集中，而序贯高斯条件模拟的数值服从概率密度函数，尊重了原始数据的分布，较好地保留了数据的变异特征；高斯序贯指示模拟可以评价空间任意位置上Hg含量大于该阈值的概率，衍生的不确定性研究结果有助于定量评估环境风险或健康风险、检测各种资源环境应用模型的稳健性以及误差传播等问题。