基于数据稳健性的农用地土壤重金属空间分异研究

刘霈珈，吴克宁，罗 明

（1.郑州大学公共管理学院，河南 郑州 450001；2.社会治理河南省协同创新中心，河南 郑州 450001；

3.中国地质大学（北京）土地科学技术学院，北京 100083；4.自然资源部土地整治重点实验室，北京100035；5.自然资源部国土整治中心，北京 100035）

1 引言

土地资源关乎亿万同胞之温饱，国家根基之稳固，然而随着中国城镇化、工业化与现代化进程的不断深入，土壤环境问题日益严重。伴随“镉大米”“重金属污染蔬菜”等一系列事件的曝光，从2014年《全国土壤污染状况调查公报》中指出全国土壤总超标率达16.1%、耕地土壤点位超标率为19.4%、以无机污染为主等土壤污染情况，到2018年8月1日起实施的《土壤环境质量 农用地土壤污染风险管控标准（试行）》《土壤环境质量 建设用地土壤污染风险管控标准（试行）》，表明土壤重金属污染问题已成为社会各界关注与研究的热点，其中土壤重金属的空间分布、变异规律等则是土壤环境质量时空信息表达的重要研究方向。

土壤重金属空间变异研究始于20世纪60年代，初期多以传统统计学方法研究土壤物化性状的空间变异结构[1]。随着地统计学的诞生与引入，研究对象由纯随机变量向区域化变量扩展、研究目的向揭示变量空间相关性转化、研究方法和研究内容都得到长足发展。土壤重金属空间分布常用克里金插值（Kriging）、样条插值（Spline）、反距离权重插值（IDW）等空间插值法，根据有限、离散的采样点信息预测区域范围内土壤重金属含量[2]。Von Steiger等早在1996年利用Kriging法绘制了瑞士东北部土壤中的Cr、Cu、Pb、Zn 4种重金属的分布图[3]。YASREBI 等比较了普通克里格法（OK）和IDW法对土壤化学性质的插值精度，认为OK模型要优于IDW[4]。但不同模型对数据结构要求有很大不同，这也会对后期面状信息模拟造成影响。如克里金插值法要求区域变量要满足平稳假设才能给出其最佳线性无偏估计（某点处的确定值）[5]。而数据稳健性处理在数据过程控制研究中常用的方法有Box—Cox转换法、Johnson转换法[6-11]。但始终缺乏一套以数据自身结构特征为出发点，更合理的重金属空间分异分析思路与方法，使得后续土地资源的合理评价、安全利用潜藏风险。

本文选取长江三角洲经济开发区、江苏省南部某市为例，以数据稳健性为突破口，尝试基于传统统计分析、数据稳健性分析、空间变异分析和空间插值分析构建农用地表层土壤重金属的空间分异研究的方法体系。依此分析研究区As、Cd、Cu、Hg、Pb、Zn这6种农用地表层土壤重金属数据特征，探讨数据稳健处理方案，分析重金属空间变异结构，对比最佳空间插值方案和空间分布差异，尝试为土壤环境质量空间分异及其后续相关研究补充更严密的研究思路、更科学的理论与技术方法、更合理的数据支撑。统计分析、稳健性分析、空间变异分析和空间插值分析构建土壤重金属空间分异研究的方法体系（图1）。

2 研究区概况与样品采集

2.1 研究区概况

研究区地势南高北低，总面积19.97万hm2，人口约108万人，其中农用地面积11.85万hm2，建设用地面积3.23万hm2，其他土地面积4.89万hm2，是名副其实的人口大市。该市气候温和、四季分明、热量条件好、降水丰沛、境内河网发达，地貌形态多样，土壤类型多样，农作物一年可两到三熟。境内已形成了陶瓷、纺织、机电、化工、轻工、建材、工艺品等富有地方特色、门类齐全的工业体系。但同时也给环境带了超负荷的压力，对土壤环境的破坏日趋严重。

2.2 样品采集

土壤样品采样区域为江苏省某市全域。根据《土地质量地球化学评估技术要求（试行）》，土壤采样点按网格化布置，农用地共布设884个点，采样密度约2个/km2，采样深度为0～20 cm，避开外来土和新近扰动过的土层，并去掉表面杂物和土壤中的砾石等五点采样法取土，四分法保留1 kg土样装聚乙烯自封袋采集样品。样品于实验室自然风干，拣出石块等杂物后送检土壤pH值、TOC、CEC和Cd、Hg、As、Pb、Cu、Zn这6种重金属。

3 方法体系的构建

土壤重金属空间分异研究的一个重要理论依据是重金属在空间上的相关性和异质性。通常要依据有限的采样点数据，选用不同的空间插值方法推演得到面状信息。由此可见，空间插值结果的合理性就与插值样点数据情况和插值方法的选择密切相关。但由于人为活动对重金属含量干扰的逐步增强，采样点数据中常会出现偏高或偏低的异常值，应通过复查排除系统误差与偶然误差，保留非偶然误差的异常值，这往往是重金属研究的重点，不能随意舍去。同时还应尝试削弱因异常值存在而显著波动的数据结构特征，避免由此可能引发的局部夸大效应，为后续空间分异研究提供更合理的数据支撑。

因此，本文以数据稳健性为突破口，基于概括性

图1 空间分异方法体系Fig.1 Method system of spatial variation

（1） 传统统计分析旨在了解样点数据自身的结构特征。本文通过分析6种重金属数据的最大值（Maximum）、最小值（Minimum）、全矩（Range）、平均值（Mean）、标准差（Standard Deviation）、变异系数（Coefficient of variation）、峰度（Kurtosis）、偏度（Skewness）等变量来估计其集中程度、离散状况和分布情况[12]。其中，变异系数通过反映重金属数据的离散变异情况，在一定程度上反映其受人为活动影响的程度。根据区域内可比原则，若CV＜15%为弱变异；若15%≤CV＜30%为中等变异；若30%≤CV＜100%为强变异；若CV≥100%则为极强变异[13-14]。同时将Kolmogorov-Smirnov正态性检验（K-S检验）置信水平设置为ɑ=0.05；若检验的P值＜0.05，则否定原假设断定总体呈非正态分布[15]。特别的，假设数据符合正态分布，利用3σ阈值识别法分析消除异常比率了解数据整体异常情况。

（2） 稳健性分析旨在对样点数据进行正态化稳健处理，弱化因异常值的存在而显著波动的数据结构，尝试寻找一种既能充分保留原有数据结构特征又符合后续空间分异分析研究要求的稳健数据变换方法[16-17]，为空间分异特征的分析提供更合理的数据支撑。本文选用Normal Score Transformation（NST）对重金属数据进行稳健处理，并对比分析稳健数据（NST数据）与原始数据（YS数据）的结构特征。

该变换方法分为三步[17]：

第一步，对n个数据集Z（μα）进行升序排序：

式（1）中：k是数据集Z（μα）第n个数据的排序序号。

第二步，升序排列后的数据集Z（μα）的累积频率可用下式来计算：若（Zμα）有相同的权重这就认为样本数据的直方图代替了研究区域的数据。若（Zμα）有不相同的权重图2）。

本文在前一种假设的情况下进行数据的转换。

第三步，序列k的数据集常态得分变换和标准累积分布函数的p*k分位数相匹配：

图2 NST原理[17]Fig.2 Theory of Normal Score Transformation (NST)

（3）空间变异分析旨在了解土壤重金属的空间变异结构特征，总结其数据结构规律，为后续进行样点信息空间可视化提供更合理的理论依据。本文利用GS + 9.0软件对土壤表层重金属原始数据（YS数据）和稳健数据（NST数据）均进行半变异函数拟合[18-19]，分析理论模型中的块金值（nugget value,C0）、基台值（sill,C0+C1）、变程（range,a）和块金效应（C0/（C0+C1））。通常块金效应可以表明系统内数据的空间相关性程度。当其小于0.25时，表明数据间具有强烈的空间相关性，以结构性变异为主，受自然因素影响为主；当其在0.25～0.75时，表明数据具有中等程度的空间相关性，受自然因素和随机人为因素共同影响；当其大于0.75时，表明数据具有较弱的空间相关性，以随机性变异为主，受人为因素影响较大[20-21]。

半变异函数的基本公式为：

式（3）中：γ（h）是半变异函数；h是样点间距；N（h）是以为间距的样点数量；Z（xi）和Z（xi+h）是区域化变量Z（x）在xi和xi+h位置上的值。

（4） 空间插值分析旨在依据重金属的空间相关性和异质性，预测未采样区域重金属含量，以较少的样点获得较准确的空间分布信息。本文选用反距离加权法（IDW）、径向基函数法（RBF）和普通克里金插值法（OK）三种空间插值法对重金属原始数据（YS数据）、稳健数据（NST数据）进行对比分析[22-23]。通过交叉验证两种数据三种方法的平均误差（ME）、均方根误差（RMSE）和平均相对误差（MRE）来对比分析其插值精度，为后续研究重金属空间分布提供基础[24]。

式（4）—式（6）中：Z（xi）为预测值；Z*（xi）为原始采样值。

4 土壤重金属空间分异分析

4.1 传统统计分析

利用SPSS 20对884个6种表层土壤重金属原始数据（YS数据）进行统计分析，结果表明（表1）。

（1）通过分析全矩R，6种重金属离散程度由大到小呈现如下排序：Zn＞Cu＞Pb＞As＞Cd＞Hg，其中Zn全矩高达2 998 mg/kg。

（2）通过分析标准差σ，6种重金属均值代表性由强到弱呈现如下排序：Hg＞Cd＞As＞Cu＞Pb＞Zn，其中Zn标准差σ高达91.8 mg/kg。

（3）通过分析变异系数CV，6种重金属变异程度由强到弱呈现如下排序： Cd＞Zn＞Cu＞Hg＞Pb＞As。所有元素CV系数均大于30%，属于强变异。其中，Cd和Zn的CV值分别高达125.54%和122.27%，属于极强变异。表明6种重金属的空间分布均受到很强的人为影响，不利于直接用于空间插值，应做进一步数据稳健分析与处理。

（4）通过分析峰度K，6种重金属分布陡峭程度由强到弱呈现如下排序： Zn＞Cu＞Pb＞Cd＞Hg＞As，其中Zn峰度高达948.66。

（5）通过分析偏度S，6种重金属分布形态均为右偏（即正偏），与正态分布相比其偏斜程度由大到小呈现如下排序：Zn＞Cu＞Pb＞Cd＞Hg＞As，其中Zn偏度高达29.61。

（6）通过对样本进行K-S检验P值均低于0.05，说明6种表层土壤重金属含量总体都不符合正态分布。

利用3σ阈值识别法识别6种表层土壤重金属异常值，结果表明数据结构合理：异常比率均小于1.5%。顺序如下：As（1.32%）＞Hg（1.23%）＞Cd（1.05%）＞Pb（0.79%）＞Cu（0.53%）＞Zn（0.18%）。Zn偏度最大，但其异常样本量仅2个，异常比率最小。As偏度最小，异常比率却相对最大，异常样本量最多。

由此可见，土壤污染研究中异常值不能随意舍去，进行数据稳健处理必不可少。

4.2 稳健性分析

利用Matlab编程运行NST数据稳健处理过程，并进行正态检验，结果表明（表1）。

稳健数据（NST数据）与原始数据（YS数据）相比，6种表层土壤重金属数据的标准差与平均值不变，变异系数CV的差异分别是：As（-0.003）、Cd（1.065）、Cu（-0.008）、Hg（0.01）、Pb（-0.004）、Zn（-0.014），表明变换后数据与原始数据变异程度几乎相同，保持了良好的数据结构特征；但峰度K均为-0.02，偏度S均为0，表明变换后数据比原始数据更近于标准正态分布，能为后续空间分异研究提供更稳定的数据支撑。

4.3 空间变异分析

利用地统计软件GS + 9.0 对原始数据（YS数据）和稳健数据（NST数据）进行半变异函数拟合，统计结果表明（表2）。

（1）通过分析决定系数R2和残差RSS，NST数据中 Cu、Pb、Zn三种重金属含量半变异函数拟合后RSS异常。通过对比原始数据拟合的R2、RSS、最大值、异常样本量和异常比率可知（表1，表2）：Cu、Pb、Zn全矩分别是578.6 mg/kg、412.9 mg/kg和2 998 mg/kg，最高异常值分别是其平均值的20倍、11倍和40倍。巨大的全矩和超高异常值与做理论半变异函数模型拟合的平稳假设不符，最终导致这三种重金属理论半变异函数拟合后RSS异常。同理，亦可解释原始数据理论半变异函数模型拟合后R2只有0.29和0.44。但经回查原始数据与实地调研结果均表明，异常值点位并非偶然误差，不可随意舍去。

（2）通过分析C0/（C+C0），除Pb和Zn外其他4种重金属原始数据与NST数据都表明其具有中等程度的空间相关性，受自然因素（如成土母质、地形、气候等）和随机人为因素（如工业企业、施肥、种植制度等）的共同影响。从两种数据对比可看出，Pb和Zn块金系数明显增大，由单纯结构性变异转向中等程度变异，表明受到自然因素与人为因素的共同作用。这也再次证明两种重金属异常值的存在，极有可能与这些异常点位附近的随机人为因素（如工业企业、施肥、种植制度等），即与潜在污染源密切相关。

表1 6种表层土壤重金属原始数据和稳健数据的统计特征Tab.1 The descriptive statistics of six heavy metals’ original data and normal score transformation data in topsoil

表2 6种表层土壤重金属理论半变异函数模型相关参数Tab.2 Semivariogram parameters of six heavy metals in topsoil

（3）通过分析全矩A，除Pb和Zn外其他4种重金属原始数据与NST数据都呈现如下排序：Cu＞As＞Cd＞Hg，且都分别是采样间距（500 m）的10倍以上，说明这4种重金属在较大范围内具有空间相关性，因此采样布点时可采用较大尺度。在既保证采样精度的情况下最大限度的节省成本，提高采样效率。尽管Pb和Zn原始数据与NST数据差异较大，却也都大于采样间距，数据分析合理。基于理论半变异函数拟合的RSS高异常，可相对适当增大采样密度，以进一步观测高异常区影响范围。

4.4 空间插值对比分析

利用GS+9.0和ArcGIS 10.2反复尝试、验证与优化，确定6种表层土壤重金属原始数据（YS数据）和稳健数据（NST数据）利用IDW法、RBF法的最佳插值参数、OK法的理论半变异函数模型后进行空间插值。

从平均误差（ME）角度分析可知，三种插值方法对6种表层土壤重金属YS数据和NST数据空间插值得到ME绝对值大小有如下逻辑顺序：

As元素：IDW-YS＞OK-YS＞RBF-YS＞RBFNST＞OK-NST＞IDW-NST

Cd元素：IDW-YS＞OK-YS＞RBF-YS＞OKNST＞RBF-NST＞IDW-NST

Cu元素：RBF-YS＞IDW-YS＞RBF-NST＞OK-NST＞OK-YS＞IDW-NST

Hg元素：IDW-YS＞OK-NST＞RBF-YS＞OK-YS＞RBF-NST＞IDW-NST

Pb元素：IDW-YS＞RBF-NST＞RBF-YS＞OK-YS＞OK-NST＞IDW-NST

Zn元素：IDW-YS＞OK-YS＞RBF-YS＞OKNST＞RBF-NST＞IDW-NST

因此，采用IDW法利用其NST数据进行空间插值得到的ME最趋近于0，在各元素中表现出明显的优越性，这种方案下插值后预测值与交叉验证样本平均值的误差量正负相抵消后较小，插值精度相对较高。

从均方根误差（RMSE）角度分析可知，三种插值方法对6种表层土壤重金属NST数据空间插值得到RMSE比对6种表层土壤重金属YS数据空间插值得到RMSE小，有如下逻辑顺序：

As元素：IDW-YS＞RBF-YS＞OK-YS＞OKNST＞RBF-NST＞IDW-NST

Cd元素：RBF-YS＞IDW-YS＞OK-YS＞RBFNST＞OK-NST＞IDW-NST

Cu元素：RBF-YS＞IDW-YS＞OK-YS＞OKNST＞RBF-NST＞IDW-NST

Hg元素：RBF-YS＞IDW-YS＞OK-YS＞OKNST＞RBF-NST＞IDW-NST

Pb元素：RBF-YS＞IDW-YS＞OK-YS＞RBFNST＞OK-NST＞IDW-NST

Zn元素：RBF-YS＞IDW-YS＞OK-YS＞OKNST＞RBF-NST＞IDW-NST

RMSE主要反映预测值与真值之间的偏差，该指标对一组样本量中的极大或者极小值有非常好的识别与表达，因此，RMSE常作为反映测量精度的重要指标。通过以上分析，利用IDW法对6种表层土壤重金属NST数据进行空间插值得到RMSE最小，表明该法插值精度相对较高。

从平均相对误差（MRE）角度分析可知，三种插值方法对6种表层土壤重金属数据空间插值得到MRE的大小有如下逻辑顺序：

As元素：OK-YS＞RBF-NST＞OK-NST＞RBFYS＞IDW-YS＞IDW-NST

Cd元素：OK-NST＞RBF-NST＞IDW-YS＞OK-YS＞RBF-YS＞IDW-NST

Cu元素：OK-NST＞RBF-NST＞IDW-YS＞OK-YS＞RBF-YS＞IDW-NST

Hg元素：OK-YS＞RBF-YS＞IDW-YS＞OKNST＞RBF-NST＞IDW-NST

Pb元素：RBF-NST＞OK-NST＞IDW-NST＞IDW-YS＞OK-YS＞RBF-YS

Zn元素：RBF-NST＞OK-NST＞IDW-NST＞IDW-YS＞OK-YS＞RBF-YS

MRE主要反映误差相对于实测值的大小，通过上述分析可知：采用IDW法对As、Cd、Cu和Hg元素的NST数据空间插值精度最高，Pb和Zn元素的NST数据用IDW法空间插值得到MRE并不低，但这也与其存在较大的异常值，数据存在突变有极大的关系。

图3 6种表层土壤重金属含量空间分布Fig.3 The spatial distributions of six heavy metals in topsoil

图4 6种重金属原始数据和NST数据三种插值组图Fig.4 Interpolation maps of six heavy metals in topsoil by using the original data and NST data

本文选用YS数据和NST数据，利用普通克里金法（OK）、径向基函数法（RBF）和反距离加权插值法（IDW）得到重金属空间分布情况（图4），并将其分别与6种重金属含量的空间分布（图3）进行对比发现结构性变异特征明显，空间分异趋势一致：As相对含量较高的集中分布在研究区西北与南部部分地区，Cd相对含量较高的集中分布在研究区西部、西南和东南部分地区，Cu相对含量较高的集中分布在研究区西北和南部部分地区，Hg相对含量较高的集中分布在研究区中部和东北部部分地区，元素相对含量较高的集中分布在研究区西部和南部部分地区，Zn相对含量较高的集中分布在研究区南部和东北部部分地区。对比结果与误差指标分析一致：利用NST数据采用IDW法插值区域识别效果最佳。

5 结论与讨论

5.1 结论

6种表层土壤重金属原始数据都不符合正态分布，均右偏且为强变异，异常比率都小于1.5%，数据合理。其中Zn的离散程度最大、均值代表性最差、峰度高度948.66，偏度最大，但是异常样本只有2个。NST变换后数据均符合正态分布，且与原始数据变异特征几乎相同。Cu、Pb、Zn三种重金属含量半变异函数拟合后RSS异常，主要原因是其数据存在巨大的全矩和超高异常值而不符合理论半变异函数模型拟合的平稳假设。但通过对比这三种重金属的块金系数发现其由单纯结构性变异转向中等程度变异，表明其受到自然因素与人为因素的共同作用，与As、Cd和Hg一致。异常值的存在与实地调研数据一致，异常点附近可能有潜在的污染源。这对重金属污染研究有重要价值，不可随意舍去。从变程来看，采样数据合理，但Pb和Zn基于理论半变异函数拟合的RSS高异常，可相对适当增大采样密度，以进一步观测高异常区影响范围。通过空间插值对比发现，使用NST稳健处理后的数据利用IDW空间插值法能较好的模拟研究区表层土壤重金属的空间分布情况。

因此，本文以数据稳健性为突破口，基于传统统计分析、数据稳健性分析、空间变异分析和空间插值分析构建的农用地表层土壤重金属空间分异研究方法体系，不仅补充了常规土壤重金属空间分异研究中对局部异常值的处理与分析方法，更试图为后续该类研究提供较为系统的研究思路。

5.2 讨论

前人诸多方面的研究常为了使数据集符合正态分布将异常值按照一定规则进行剔除，这在一定程度上破坏了数据集本身的结构特征，对于表层土壤重金属空间分异这类极易受外界影响的研究而言，异常值的处理应更谨慎，不能按照常规做法任意的抛弃和修改。

数据稳健处理方法多样，但要以能合理保障、客观反映数据集自身结构特征为根本依据。本文中NST变换法虽可将数据进行完美的正态化处理，但却具有只适用于只有正值的数据集的局限性。因此，异常值的甄别与处理、数据稳健处理方法、空间变异分析思路等方面的研究决定了后续重金属调查、评价与安全利用研究的基础数据之准确性与合理性，对表层土壤重金属空间分异和分布研究至关重要，还有待更进一步深入。

致谢：感谢江苏省地质调查研究院土地质量评价与污染防治应用研究中心为本文撰写过程提供的帮助。