一、问题的提出

农田土壤污染状况详查是“土十条”提出的排在首位的重要任务。根据“土十条”提出的要求，2018年底前需要查明农用地土壤污染的面积、分布及其对农产品质量的影响。这是一项相当艰巨的任务，因为国土资源部也曾组织协调全国77家单位10万多人次，按照统一技术规范实施了全国土地地球化学调查，即多目标区域土壤地球化学调查，历时15年，仅完成全国耕地中面积的68%调查任务，即13.86亿亩耕地【1】。如果扣除调查后期数据成果整理时间，能用在全国土壤污染状况详查的野外采样和实验室检测的时间非常有限，如何在在不损失调查精度的前提下提高效率，是当下面临“时间紧、任务重”双重压力的全国土壤污染状况详查阶段所亟待解决的技术难题。

根据《全国土壤污染状况调查技术规定》，农田土壤污染状况调查采用网格布点方法，并结合土地利用类型适用不同网格布点密度，其中耕地调查网格布点密度为8km×8km=64km2/点位,林地和草地调查网格布点密度为16km×16km=256km2/点位,毫无疑问，详查阶段的网格布点密度（精度）将大于调查阶段，参照固体矿产勘查阶段划分及其技术规程，土壤详查阶段网格布点密度可能呈倍数增加，例如，农田土壤调查阶段的网格布点密度64km2/点位的倍数布点密度依次为32km2/点位,16km2/点位, 8km2/点位, 4km2/点位, 1km2/点位,所对应的调查点位数量（实物工作量）则呈几何级数增加，即相当于调查阶段工作量的2倍，4倍，8倍，16倍，64倍，毫无疑问，网格布点密度越大，其调查精度越高，但所需要的调查时间也越长。因此，如何平衡土壤详查阶段的精度和效率，是当下面临“时间紧、任务重”双重压力的全国土壤污染状况详查阶段所亟待解决的管理问题。

本质上，全国土壤污染状况详查是基于统计学理论的抽样实践，旨在通过有限样本对全国性土壤环境质量现状这一真实总体进行统计学意义的估计，毫无疑问，随着样本数量的增加，统计期望值与真实总体之间的偏差将缩小，即精度提高。但是，如何定量评价调查精度？例如，耕地土壤详查有5种可选网格布点密度，依次为 32km2/点位,16km2/点位, 8km2/点位, 4km2/点位和 1km2/点位，究竟应采用那种布点方案？其所对应的精度是多少？与增加的实物工作量相比较，那种布点方案更优？即网格布点优化，是当下面临“时间紧、任务重”双重压力的全国土壤污染状况详查阶段所亟待解决的理论问题。

事实上，国土资源部早在15年前就已经开展以农田土壤地球化学调查为主要内容的全国农业地质调查工作，并在全国范围采用统一网格布点方法，其中表层土壤（0-20cm）采样密度1km×1km，深层土壤（150-180cm）采样密度4km×4km，与全国土壤污染状况调查方法基本一致，如果以全国土地地球化学调查数据为参照，定量评价不同网格布点调查方案的相对精度将成为可能。本研究拟以浙江省杭嘉湖地区农业地质调查所获取的6335个点位土壤地球化学调查数据为基础，应用地形图缩编中高程点自动抽稀方法【2】，重构具有不同网格布点密度的数据集，进一步应用趋势分析和地统计学等地质数学方法，定量评价不同网格布点方案的相对精度，以期为全国土壤污染状况详查之网格布点方案优化提供决策依据。

二、数据准备

本研究拟采用“按范围控制自动抽稀”方法重构土壤调查数据集，即将杭嘉湖地区农业地质调查工作底图（1:25万）分割成等间距的网格，网格大小依次为256km2,128km2, 64km2,32km2,16km2，8km2和4km2 ，依次对原始数据（1km×1km）进行自动抽稀，即，每个网格只保留一个点位数据，其余点位的数据全部删除，经过七次自动抽稀，即可构建具有7种网格布点密度的土壤调查数据集，相当于采用不同网格布点方法对杭嘉湖平原进行了七次调查，与原始数据（1km×1km）进行比较，即可获得不同网格布点的相对调查精度。为保证统计数据的可靠性，本研究拟选择测试精度较高，区域性规律较明显的表层土壤全氮含量为研究对象，不同批次抽稀后表层土壤全氮含量空间分布如图1-图8所示。

原始数据自动抽稀结果显示，随着土壤采样密度变稀，杭嘉湖平原地区表层土壤采样点越来越少，采样点位个数由原始的6635个递减到12个（网格面积512km2）【3】。原始数据每抽稀一次，相邻两种网格密度所构造数据集的区域等值线分布图都有细微变化。但从大尺度的区域土壤全氮含量分布趋势看，随着网格布点密度变稀，其空间分布信息越模糊。当网格布点密度为1/8km2和1/16km2时，土壤全氮含量的空间分布特征与原始数据分布特征非常相似；当采样点为1/32km2和1/64km2时，得到的土壤N的空间分布信息丢失比较严重，仅能反映其趋势变化；而网格布点密度变稀至1/128km2以下，表层土壤全氮含量等数值等值线图开始有明显变化，高密度小尺度采样方法所得到的真实值远远高于低密度大尺度采样后的估计值，说明所得到的土壤N含量分布格局与实际情况具有明显的差异。

三、统计分析

1、数理统计

原始数据经过7次抽稀所构建的数据集点位数量差异显著，为便于开展点对点的精度比较，需要对不同数据集进行空间插值（克里格法），成对数据t检验结果如下表。

表1  不同采样密度下土壤全氮含量成对数据T检验表

统计检验结果表明，不同采样密度多对应的数据集间有显著差异，说明网格布点方案对总体估计的影响极其显著，不同批次调查数据进行比较时需要特别说明。

表2  不同网格布点密度表层土壤全氮含量方差分析

同网格布点密度所对应数据集方差分析结果如表2，计算结果表明，网格布点密度大于耕地土壤污染状况调查密度（8km×8km=64km2/点位）的网格面积依次为32km2,16km2, 8km2,4km2，其所对应区域土壤全氮含量均值依次为1436.5mg/kg,1448.4mg/kg,1405.9 mg/kg,1413.1 mg/kg, 8km×8km点位密度对应的全氮含量均值为1402.9mg/kg,不同网格布点密度所对应的土壤全氮含量均值差异不大，说明上述网格布点方案具有一定的可比性。

2、空间半变异函数

利用不同网格布点密度所对应数据集进一步构建空间空间变异函数，结果如图9所示。

图9  不同网格布点密度所对应的土壤全氮含量空间变异函数

图9结果显示，当网格布点密度大于8km×8km=64km2/点位,即网格面积为4km2、8km2、16km2和32km2所对应数据集的变异函数曲线几乎重合，说明这4种采样密度都能够很好的反映原杭嘉湖平原地区土壤全氮含量的空间分布情况。当网格布点密度为64km2/点位时，变异函数开始发生变化，其基台值增大，但仍然能用相同的指数模型拟合。但若网格布点密度小于64km2/点位，即网格布点密度降低至1/128km2、1/256km2和1/512km2时，其变异函数曲线开始发生较大偏离，不能用相同的理论模型进行拟合，无法反映原始数据的空间变异情况。由此可见，64km2/点位网格布点密度可能就是杭嘉湖地区表层土壤全氮含量空间变异的拐点。换而言之，杭嘉湖地区表层土壤调查网点密度不应低于8km×8km=64km2/点位。

3、趋势面分析

趋势面分析是拟合数学面的一种统计方法，即用数学方法计算出一个数学曲面以拟合数据中的区域性变化的"趋势"，趋势面和剩余面构成，前者反映宏观分布规律，是确定性变化的结果；后者反映微观变化，是随机变化的结果。杭嘉湖地区表层土壤全氮含量趋势面分析结果如图10-图16所示。

图10 杭嘉湖地区表层土壤全氮含量之随机面与趋势面（网格布点密度=1/4km2）

图11 杭嘉湖地区表层土壤全氮含量之随机面与趋势面（网格布点密度=1/8km2）

图12 杭嘉湖地区表层土壤全氮含量之随机面与趋势面（网格布点密度=1/16km2）

图13 杭嘉湖地区表层土壤全氮含量之随机面与趋势面（网格布点密度=1/32km2）

图14 杭嘉湖地区表层土壤全氮含量之随机面与趋势面（网格布点密度=1/64km2）

图15 杭嘉湖地区表层土壤全氮含量之随机面与趋势面（网格布点密度=1/128km2）

图16 杭嘉湖地区表层土壤全氮含量之随机面与趋势面（网格布点密度=1/256km2）

杭嘉湖地区表层土壤趋势面分析结果表明，随着网格布点密度变稀少，估计值与实际值误差变大，趋势面发生形变，剩余面变化趋于平缓，局部异常和土壤全氮含量的微小变化难以发现,换而言之，详查网格布点密度过稀，土壤污染（局部异常）难以查明，若是土壤污染详查，其所对应的污染面积可能被夸大，也不利于定位污染位置。

4、相关分析

为定量研究不同网格布点密度所对应的调查精度，本研究拟以1/4km2网格布点密度为参照，将不同网格布点密度所对应的数据与之相对应，并进行相关分析，结果如图17。

相关分析结果表明，随着网格布点密度变稀，插值得出的土壤全量N值与原采样尺度下的值之间的线性相关系数R2逐渐减少，从采样密度为1/8km2的0.7643递减到1/256km2的相0.1823。若以R2＝0.50作为抽稀后采样密度下土壤N值与原采样密度土壤N值相关性紧密程度的判断标准[4]。网格布点密度1/8km2、1/16km2、1/32km2和1/64km2的土壤全氮含量与实测数据（网格布点密度=1/4km2）的相关系数R2均大于0.5，说明上述四种网格布点密度所对应的采样数据与实测数据高度相关，具有较高的精度，反之，网格密度1/128km2和1/256km2以及1/512km2的采样密度的数据结果已经与原采样结果不存在线性关系，说明上述3种网格密度所对应的土壤全氮量数据在较大程度上已经偏离集真实值，不能反映土壤全氮含量的空间变化趋势。由此可见，64km2可能是区域土壤详查网格布点密度的临界值，在平衡精度与实务工作量的前提下，8km2/点位可能是土壤详查的较佳点位密度。

四、结论

考虑到全国土壤污染状况调查中耕地土壤布点密度为8km×8km=64km2/点位，可将杭嘉湖平原原始数据抽稀后1/64km2所对应数据集视为“调查数据”，大于该网格布点密度的抽稀数据集视为详查数据，网格面积32km2,16km2, 8km2, 4km2所对应的“详查数据”均值变化幅度依次为2.40%，3.24%，0.21%，0.73%，该变化幅度可视为“详查精度”，其中网格布点密度8km2/点位的精度最高，其对应的实物工作量增加8倍。若以8km2/点位密度为界，以32km2和16km2为网格布点密度所对应的调查精度平均值为2.82%，以8km2和 4km2为网格布点密度所对应的精度平均值为0.47%，与前者比较，精度提高6倍，实物工作量增加4倍。

参考文献

1、中国耕地球化学调查报告（2015），中国地质调查局.

2、李乃良.基于Autocad的地形图缩编中高程点自动抽稀方法的探讨.城市勘查[J]，2007，3：81～82.

3、陈有鑑，杭嘉湖地区浅表层地质环境综合监测项目成果报告（2007），国家级地质环境监测专项工作，中国地质环境监测院.

4、[15]王政权.地统计学及在生态学中得应用.科学出版社[M].1999.

备注：该文首发于【中国环境修复网】http://www.hjxf.net/2017/0329/3076.html，首发日期：2017-03-29 10:49:34