农产品产地环境风险防控

研究土壤污染物的空间分布及对人体健康的影响，首先应保障各种插值技术的正确应用，保障估计土壤污染物空间变异性的精度。如果无法正确评估土壤污染物的分布情况及各区域估计值的准确性，就无法保障正确评估土壤污染对人体健康的影响。因此，提高土壤污染物空间变异性估计的精度，是评估土壤污染对人体健康影响的基础。

提高地质统计学估计精度是一项系统工程，应从准备工作到最终结果输出、验证等全过程加以掌控，其关键步骤如图3-2所示。

图3-2　地质统计学估计土壤污染物空间变异性精度控制步骤图

3.2.6.1　准备工作

（1）基础资料的准备　应用地质统计学研究土壤污染物空间变异，需做好充分的准备工作，首先要收集影响污染物空间分布的基础资料。这些资料的收集对于评估差插结果的精确性、正确性、分析污染物来源、解析空间变异的原因至关重要。

收集的资料主要包括以下内容：研究区的地形地貌信息、土壤信息，主要包括坡度、坡向、河流分布、流向、地质条件、成土母质情况、土壤类型、土壤肥力情况、土壤质地等；污染物、污染源情况，包括区域污染源类型、分布、排放的污染物类型、去向、数量、污染物排放历史、处理情况、污染物进入研究区的利用情况如污灌等；土地利用方式，如作物种植类型、耕作方式、收获方式、土地利用类型分布差异性等；农业投入品情况，如畜禽粪便使用情况，化肥、农药使用情况等；研究区的行政区划情况，如县界、乡界等；特征地物情况，如主要河流、居民地分布、影响样点数据的收集或者未来样点调查的物理或生物属性；遥感资料；气候、气象资料；研究区以往的环境质量调查、监测资料，已发表的与研究区相关的其他资料；其他可能影响研究结果的定性和定量信息。

只有了解了区域特点才能进行正确的布点采样，揭示污染物在不同尺度的变异特点，避免布点密度过稀无法分析微域变异性或布点过密浪费人力、物力。

（2）采样器具准备

工具类：铁锹、铁铲、圆状取土钻、螺旋取土钻、竹片及适合特殊采样要求的工具等。

器材类：GPS、罗盘、铝盒、样品袋、样品箱、手机及备用电池等。

文具类：样品标签、采样记录表、铅笔、资料夹等。

安全防护用品：工作服、工作鞋、安全帽、药品箱等。

3.2.6.2　布点采样分析

土壤环境质量调查中，所获取的信息量将随样本数的增加而增加，但信息量的增大是以成本的急剧增加为代价的。如何找到成本与信息量的最佳结合点，就目前的研究看来，尚无法给出明确的答案。但在布点采样前应把握以下几个问题：

一是了解调查区信息，具体需了解的情况可参见本节3.2.6.1准备工作的相关内容。

二是掌握对于数据质量的要求。①数据收集的目的。对于本研究而言，其目的就是了解土壤重金属的空间变异性，了解区域风险的差异。因此，采样方法的设计，采样密度等的设计都应围绕这一中心问题进行。②主要的评价指标。采集的样品数量应能满足实际分析需要，应具有实际统计意义，能满足评价指标的要求。③关心的统计参数。如平均值、方差、中值、百分位数，变异系数等。④对于精度的要求。精度要求高，则采样密度就高；精度要求低，则采集的样品密度相对较低。《土壤环境监测技术规范》（HJ/T 166—2004）给出了基本采样量的计算公式。⑤采样尺度的要求。如对于大尺度的研究，则采样点距离相对较大，而对于田块级的研究，其采样尺度就相对较小。

三是熟悉采样的限制因素。采样条件的限制，如采样人员的限制、技术条件的限制，时间及计划的限制，经费的限制，地理、地质条件的限制，其他的限制。

进行地质统计学布点采样应根据以下步骤设计采样方案：①研究监测要获得的结果。②收集相关的布点采样方案，研究各方案的优缺点及对监测任务的适用性。③用数据表达各布点采样方案的实行及成本核算。④根据布点采样精度及预算的要求，计算各采样方案需要的采样点数量。⑤选择最有效且经济可行的布点采样方案。⑥制定方案实施计划。

在众多的土壤采样方法中，大致可以归纳为以下几种：

（1）专业判断布点方法　采样点的位置、采样数量及采样时间的确定基于对要调查对象特征的了解、判断及以布点采样者的专业知识。专业判断采样方法是一种基于专业知识而非统计学理论的主观采样方法。因此，用专业判断法获得的监测样本对于总体的代表性，完全依靠采样者或采样布点设计者对监测区的了解及专业判断能力。

专业判断采样方法主要适用于以下情况：①对于相对较小范围或尺度的调查。②采集及需分析的样点数量少。③对于调查区有较充分的了解和认识。④调查某种污染物在某区域内是否存在或超过某一阈值的情况，如果调查发现污染情况存在，则还应进行进一步的布点采样分析。

对于地质统计学而言，本方法一般适用于二次采样的情况，即通过地质统计分析后，对于估计方差较大的区域应进行二次布点采样分析时，专业判断采样法一般能满足要求。其优点为：一般能符合采样计划及经费限制；其缺点为无法对调查采样的不确定性进行评估，无法外推未采样区的情况，无法对总体进行统计意义上的无偏估计。

（2）简单随机采样　本方法采用随机数发生器选择采样点的位置，所有选择的位置都具有相同的概率，连续的采样位置的选择完全独立于前面的选择。适用条件为采样区条件较均一，且用于地质统计学多次采样的前期布点采样工作。该采样方法提供了统计学意义上的总体平均值和方差的无偏估计，便于理解和实施，采样点数量的计算及数据统计相对简单，可采用均方差和绝对偏差法或变异函数及相对偏差法估计采样点数量。但由于所有的采样点都是随机选择，因此有可能出现采样点团聚现象，且无法利用其他有利于提高采样精度的辅助信息。

（3）分区随机采样法　分区随机采样法是也称划分监测单元采样法或分层采样法，是对随机采样法的改进。本方法在采样前首先根据已掌握的调查区的基本资料，将调查区分为几个独立的采样单元，每个独立的采样单元其土壤的理化性质（如土壤类型、土壤质地等）、农作物种类、耕作制度、污染物含量及分布情况应相对均一。在每个采样单元内根据以往的调查结果及研究的实际情况确定采样点数量，采用随机布点采样法布点。分区随机采样法对调查区总体的估计精度一般较完全随机采样法要高。

适用条件为：调查区影响污染物分布的因素存在较大变异性；对调查区基本情况较了解，能正确划分监测单元；调查区范围相对较大；调查者希望对调查区内不同的子区情况分别进行统计分析；对调查区内不同的监测单元采用不同的布点采样方法。

该方法对不同的采样单元可以布测不同密度的监测点数量，从而提高了监测的精度，降低了成本，增加了监测的灵活性；可以突出对某一特征的监测，如对水污染引起农田污染的监测，大气污染引起农田污染的监测，固体废物污染引起农田土壤污染的监测等；可以对某一小尺度的问题单独划分监测单元，以突出其对总体的影响或对其进行更深入的研究。

但监测单元的划分对于正确的布点采样及分析结果的代表性有很大的影响，因此要求监测单元的划分者应对监测区有必要的了解，在划分监测单元前应充分收集监测区的相关资料。监测精度的提高与划分监测单元所需辅助信息高度相关，因此对于划分监测单元的辅助信息要求较高，从而在一定程度上增加了信息收集的难度。与随机采样一样，采样者实际的采样位置与设计的采样位置有可能有一定的偏差，从而影响了监测结果的精度。分区（或划分监测单元）随机采样法监测结果的统计分析工作相对较烦琐，一般不易掌握。在一个监测区内，为便于统计结果的计算，监测单元数量一般不应超过6个。

（4）系统/网格布点采样法　系统采样法也可以称为网格采样法或规则采样法，在这一采样法中，样本的位置或是处于规则网格的中心，或是在网格线的节点上，或是在网格内部随机布点。网格可以具有不同大小和形状，最普通的形状是正方形、长方形、六角形和三角形网格，一般认为三角形网格较其他方法更有效。如果被监测的污染物空间变异性存在各向异性时，可考虑采用长方形网格。

适用范围为：对监测区进行全面的调查，尤其当被监测的污染物质存在空间相关性时，系统采样法较随机采样法具有更高的精度；评价污染物随时间及空间的变化趋势时；调查监测区内的“热点”区域。

该方法可以覆盖整个监测区；布点、采样设计直观、简单；可以对不同尺度的网格进行套合，以便对重点区域进行更深入的研究；如果被监测的污染物存在时间或空间的连续性，而采样尺度小于被监测的污染物的时、空变异范围，则系统采样法对污染物的时、空变异性的反映优于其他采样方法；采样者不必对监测区有深入的前期了解。但对于监测区的范围、面积及采样点数量必须清楚。由于采样行或列可能与所关心污染物的变化特征重合，从而引起对总体特征的不正确估计。如果采用的是完全等间距采样，当总体有某种周期性的规律，而采样间隔又恰好与这种规律周期吻合的时候，则采样将因无法反映规律性而失去代表性，因此系统采样法应尽量避免采用行或列、等间距采样方式，而应采用划分网格在网格内随机采样的方式。如果对监测区有较多的了解，熟悉监测区污染物的分布特征及分布规律等，系统采样法有可能不是首选的方法，因为此方法采样点数量要求较大，消耗的人物、物力、财力较其他方法多。如果网络布设的尺度较大，有可能对于一些“热点”地区无法捕捉。

（5）划分监测单元的系统采样法　划分监测单元的系统采样法将监测单元采样法与系统采样法进行了融合，首先对研究区域根据以往的调查资料划分监测单元，然后分别在监测单元内采用系统采样法进行布点，但布点网格的设计可根据监测单元的情况不同而不同。本方法集中了监测单元采集法与系统采样法的优点，适用于进行大区域土壤环境质量调查。

采用划分监测单元进行土壤样品的采集，虽然在整个研究过程中会降低成本，但其问题在于此方法会使土壤采样点代表的面积有一定的差异。如对于已知的污染区，土壤监测样点的数量较一般没有污染的区域要多，从而污染区样点的代表面积就会小于没有污染区样点的代表面积，如果这种代表面积的差异较大，则在应用地质统计学前应对监测数据进行降聚处理，从而消除丛聚效应对估计结果的影响。

消除数据丛聚效应的方法有很多，如泰森多边形法或称格网降聚法（Thayer，2003；刘鑫等，2011；于振汉等，2011），本方法核心就是给丛聚在一起的地质数据分配较小的权值，给稀疏分布的数据分配较大的权值。如图3-3所示，图中面积的大小代表监测点权重的大小。

图3-3　消除监测数据丛聚效应的泰森多边形图

权重的计算公式为：

式（3-54）中，wi为分配给监测点i的权重系数；ai是监测点的泰森多边形面积；A是区域总面积。

由于克里金法本身具有分团效应的特点，因此如果监测数据的代表面积不是差异很大时，也就是说监测点分布比较均匀时，可不必对原始数据进行降聚处理。

3.2.6.3　实验室分析质量控制及分析质量保证

实验室分析质量控制服从于研究目的，大范围的产地环境质量普查及调查是在一个较大空间内研究土壤污染物的含量及其分布行为，研究土壤污染对农产品质量的影响，对样品的分析的准确性要求有足够的稳定性和可比性。因此，在分析质量控制的程序构成上，主要包括实验室内部分析质量控制和实验室之间的分析质量控制，并包括协作系统试验而规定的统一分析方法。

一般实验室土壤环境质量普查和农产品质量调查分析质量控制基本步骤如图3-4所示，实验室样品检测流程如图3-5所示。

3.2.6.4　探索性空间分析

探索性空间分析是对数据分布情况、数据空间相关情况、数据的比例效果等进行统计分析研究，以初步判断监测数据是否具有空间相关性，是否满足某种插值对数据自身要求的一种整体探索。

（1）异常值的处理　监测数据含有异常值，对半变异函数的影响很重要，特别是在变程α之内的异常值影响实验半变异函数的计算，引起块金值和基台值的明显增加，从而影响估值的精度。影响半变异函数的异常值分为全局异常值及局部异常值两种。全局异常值是对全部监测数据而言的，这种异常值可通过常规异常值检验方法查到，如采用格鲁布斯检验、T检验、狄克逊检验等查找，如果采有非参数方法，也可以采用频率柱状图查找，在频率柱状图特别是偏态分布的频率柱状图中两端的数据点或孤立点，应做进一步分析，根据实际情况剔除或作掩膜处理。对于以上方法在使用中应注意以下几点：①在应用上述方法判断异常值前，应先检查原始记录是否按规定的要求填写完全、正确；若查明是过失或错误的数据，如采样、分析、测定操作错误或发生意外污染等，应予舍去。②以上检验方法仅适用于正态分布的总体样本，若来自对数正态总体，应先将数据取对数，然后用对数数据样本再实施以上方法。③狄克逊检验仅适用于样本容量大于25的样本，其他两种参数方法大、小样本都适用。④实际应用中，往往采用两种以上的方法综合判断数据的异常值。⑤当样本中异常值不止一个时，应逐个判断，逐个剔除。⑥异常值可能来源于外来污染，如果是由于土壤污染造成的异常值则不能剔除，而应采用掩膜的方式，即在计算实验半变异函数时将其剔除，而实际插值时考虑异常值。

图3-6是采用天津市东南郊砷元素监测数据制作的频率柱状图，可见图中的A、B两个柱所代表的监测数据可判断为异常值，但其剔除与否，还应根据实际情况决定。通过实际分析可知，这两个柱所代表的两个监测点是由于测试错误造成的，应剔除。

图3-4　实验室分析质量控制程序

而局部异常值是对局部区域而言的，局部异常值对于全局数据而言不表现为异常。因此，采用传统统计方法难以发现监测数据的局部异常值。局部异常值可通过H-Scatter法等分析。H-Scatter图可以分析数据和距离之间的关系。根据地质统计学的原理，位置近的点，性质就越相似，数据的变异性就越小，数据点就越易靠拢于斜率为1的直线，而如果相邻点变异性大，数据点就会偏离分散，表现为局部异常值。

图3-5　实验室样本检测流程图

由图3-7可见，与图中A、B、C相关的数据很有可能为异常值，应在分析中剔除或做掩膜处理。

全局异常值也可以通过半变异函数云图查找，如图3-8所示，全局异常值在半变异函数云图上常表现为与总体分离的点（图中小方框所示）。根据半变异函数的性质可知，随着距离的增加，半变异函数值也随之增加，半变异函数云图中的点值也随之增加，但如果在较近距离内有的点的半变异值也很高，则可判定为局部异常值。(www.xing528.com)

采用半变异函数云图法，H-Scatter图法和频率柱状图法相结合的方法判断监测数据的异常值，是推荐的方法。因为这样得出的结果相互验证，更少发生错误的判断。

图3-6　频率柱状图法检验全局异常值

图3-7　H-Scatter方法检验数据局部异常值

图3-8　半变异函数云图

除以上几种方法判断数据的异常值外，也可以采用局部莫兰指数（Moran’s I）法进行空间聚类与异常值分析。本方法可以对数据进行更全面、宏观、直接的展示，以分析异常情况。

给定一组要素（输入要素类）和一个分析字段（输入字段），聚类和异常值分析工具可识别具有高值或低值的要素的空间聚类。该工具还可识别空间异常值。为此，该工具计算局部Moran’s I值、z得分、伪p值和表示每个具有统计显著性的要素的聚类类型的编码。z得分和伪p值表示计算出的指数值的统计显著性，其计算公式为：

统计数据的ZIi得分的计算方法如下：

其中，

正值I表示要素具有包含同样高或同样低的属性值的邻近要素；该要素是聚类的一部分。负值I表示要素具有包含不同值的邻近要素；该要素是异常值。在任何一个实例中，要被视为具有统计显著性的聚类和异常值，要素的p值必须足够小。

要提高异常值、空间聚类的精确性，则应遵循以下两项基本准则：一是统计分析的数据点位数不应小于30。二是要选择适当的空间距离，基本上所有要素都应至少具有一个相邻要素，任何要素都不应将其他要素作为相邻要素，尤其是在值有偏斜时，每个要素都应具有8个左右的相邻要素。

（2）数据正态性检验　普通克里金、简单克里金、泛克里金、序贯高斯条件模拟值法等方法大都建立在监测数据呈正态分布基础之上，特别是条件模拟值法要求数据严格服从正态分布，而普通克里金法要求数据服从近似正态性即可（析取克里金方法要求数据服从严格正态分布，指示克里金不要求数据正态性）。但数据分布过于偏态，会使半变异函数缺乏稳健性，影响分析精度。图3-9为天津市东南郊土壤Cu含量频率柱状图，图3-10为土壤Cu估值结果的交叉验证估计误差图。由图可见，Cu元素估计随Cu含量增加而增加，表明估计结果产生了系统误差。因此，在数据插值前，应对数据进行分布检验。采用条件模拟值法进行空间插值，可采用K-S法、峰度—偏度法等参数方法进行检验，采用克里金插值法进行空间插值，可采用直方图法、Q-Q图或P-P图等直观方法检验数据分布情况。采用条件模拟值法时如数据不符合正态性，应首先采用非线性变换法对原始数据进行正态变换，如果采用克里金插值法，则可通过对数变换、三角函数变换、平方根变换等方法使原始数据接近于正态分布。

（3）比例效应　比例效应可粗略地认为，当样品的平均值增加或减少时，样品的方差也随之增加或减少。判断比例效应是否存在主要是通过移动窗口法分析平均值与方差或标准差之间的关系。如果平均值和标准差之间存在明显的线性相关性，则比例效应存在，否则，比例效应不存在。消除比例效应的方法主要是通过对原始数据取对数或采用相对变异函数的方法（王政权，1999）。

（4）空间自相关分析　要计算区域化变量的半变异函数或所选择克里金插值及条件模拟值技术内插污染物空间分布表面，则应首先分析数据是否存在空间相关性。如果数据不存在空间相关性，则插值仅能用确定性插值方法，如全局多项式、局部多项式、反距离插值等。

图3-9　天津东南郊土壤Cu元素频率柱状图

图3-10　土壤Cu估值交叉验证估计值与估计误差散点图

空间自相关有正相关与负相关两种情况，如果相邻两点上的值均高或均低，则称其为空间正相关，否则称为空间负相关。空间自相关分析在地理统计学科中应用较多，现已有多种指数可以使用，但最主要的有两种指数，即Moran的I指数和Geary的C指数。与统计学上的相关系数相近，Moran的I指数其值变化于0和±1之间。当I＝0时代表空间无关，当I＞0为正相关，而I＜0时为负相关。Geary的C指数取值一般在0～2，大于1表示负相关，等于1表示不相关，小于1表示正相关。Moran的I指数和Geary的C指数呈负相关，基本上是线性的趋势。

区域化变量在某个滞后距离内空间自相关是否显著，还需要经过显著性检验才能确定。ArcGIS空间数据统计之分析模式工具中的空间自相关及DPS软件中的空间自相关分析都给出了详细的统计分析过程（唐启义等，2006）。还要注意的是，空间自相关（Global Moran’s I）全局统计量用于对数据的总体模式和趋势进行评估。如果空间模式在研究区域内保持一致，这些全局统计量最有效。局部统计量[如热点分析的（Getis-Ord Gi*）工具]用于在相邻要素的环境下对每个要素进行评估，然后将局部情况与全局情况进行比较。请考虑一个示例，当在计算一组值的平均值时，还要计算某个全局统计量。如果所有值都接近20，则平均值也将接近20，并且该结果可以非常好地表示或概括整个数据集。但如果一半值接近1，而另一半值接近100，则平均值将接近50。可能不存在任何接近50的数据值，因此，该平均值并不能很好地表示或概括整个数据集。创建数据值的直方图时，将全是双峰分布。类似地，当所估量的空间过程在研究区域内保持一致时，使用全局空间统计量[包括空间自相关（Global Moran’s I）工具]将最有效，这样所得到的结果将能很好地表示或概括总体空间模式。

3.2.6.5　实验半变异函数计算及拟合

（1）实验半变异函数计算应遵循的原则　用来计算实验变异函数值的数据量应足够大，一般要求在30个数据点以上；在每个分离距离h上用来计算实验变异函数值的数据对数n（h）最少应大于20，如果用于计算半变异函数步长值的样点对太少，则对于短距离的变异估计会缺乏代表性。如图3-11所示，对于第一个半变异函数值而言，仅有一对样点值参加了计算（如图3-11B所示），显然这样计算的半变异函数其代表性也较差，且估计方差增大（张仁铎，2005）。

图3-11　半变异函数值及计算样点对图

注：右图方框内的值为异常值，未参与运算。

只用分离距离|h|≤L/2的实验半变异函数值来拟合模型。L是研究区内沿一方向的最大尺度。主要是因为局部变异函数的涨落方差与研究域的大小有关，当|h|≥L/2时，相对涨落方差约为200％，造成估计误差增加。

在实际工作中，因观测点并不完全呈规则的矩形网格分布，所以在求某个方向上的实验变异函数时，由于观测点不一定完全位于这个方向的同一条直线上，计算半变异函数时采用角度允许误差限和距离允许误差限的方式进行调整。角度允许误差限在数据处理程序中给定。若角度允许误差限为10°，要计算90°方向上的实验变差函数，则从某一点出发，对位于80°～100°扇形区域内的任一点都可以看成是在90°方向上的点。距离允许误差限取基本滞后距的一半，若基本滞后距为400m，要求两点间的距离为1 000m的实验半变异函数值，则接点相距800～1 200m范围内的任一点都可以近似看成是这两点间的距离为1 000m。

如果区域化变量不平稳造成实验半变异函数存在趋势，应首先将趋势剔除后，再计算剩余的实验半变异函数，而在插值时将剔除的趋势反加进来。分离距离h应尽可能小，以反映小尺度的变异性。如果实验半变异函数存在各向异性，则在计算各向同性实验半变异函数的同时，还需计算各向异性实验半变异函数。计算得到的实验半变异函数应避免点值之间较大的起伏，实验半变异函数曲线尽量平滑。实验半变异函数应尽量反映小距离范围内的变异性。实验半变异函数尽量避免翘头显现，以免影响实验半变异函数的拟合，特别是采用自动拟合时，这种影响更显著，如有可能将球形模型拟合成高斯模型等。

（2）半变异函数平稳性检验　根据线性平稳半变异函数的性质可知，如果半变异函数满足二阶平稳假设或本征假设，则r（∞）＝C（0），C（0）为半变异函数先验方差。此性质说明，如果半变异函数平稳，实验半变异函数最大值不应超过先验方差。在实际经验半变异函数常出现3种情况，一是实验半变异函数值先增加然后达到一稳定值（不大于先验方差），以后随着距离的增加，半变异函数值不再增加，这种情况说明半变异函数符合二阶平稳假设或本征假设。二是实验半变异函数值先增加然后达到一稳定值（不大于先验方差），以后随着距离的增加，半变异函数值又开始增加，表现为全局不平稳性。G.马特隆研究表明，如果在变程α范围内有足够多的数据点，就可以不考虑全局不平稳性，半变异函数符合准二阶平稳假设或准本征假设，在半变异函数拟合时只考虑平稳部分计算结果而忽略不平稳（漂移）部分。三是实验半变异函数值随着滞后距离的增加超过先验方差无法实现平稳而一直表现出增加趋势，这种情况说明半变异函数不平稳。

（3）实验半变异函数拟合　在拟合过程中，应将更多注意力集中在较小距离的实验变异函数值上；如果区域化变量Z（x）在变程α之外表现为不平稳，则半变异函数拟合时可只考虑变程α之内的半变异函数值，变程α之外的漂移可忽略。但如果区域化变量Z（x）在变程α之内表现为不平稳，则应先剔除不平稳量（即趋势）后再用正则定模型拟合剩余的实验半变异函数或用协方差函数描述区域化变量的空间变异性。对于各向异性的实验半变异函数或较难用一个正则定模型拟合的实验半变异函数，可通过相关模型的套合实现半变异函数的拟合。

拟合模型的方法一般有实验法（也称目试法）、广义最小二乘法、加权最小二乘法、规划法、遗传算法、自动基台值拟合法等。这些拟合方法各有优缺点，各个相关模型所选用的拟合方法也不尽相同。Surfer软件选择的是最小二乘法，vairance 2D软件采用广义最小二乘法，而ISATIS软件采用基台值拟合法。基台值拟合法较适用于各向同性的情况，对于各向异性的情况还需采用实验法进行人工矫正。

实验半变异函数拟合是一个复杂的过程，虽然各国科研工作者提出了一些较好的拟合方法模型，但这些模型都有其局限性，无法做到完全的自动拟合。因此，采用人工干预下的实验半变异拟合模型拟合方法，应该是实验半变异函数拟合最有效的方法。

ArcGIS地质统计分析工具、Surfer、IASTIS、vairance 2D等软件在实验半变异函数拟合上都有一定的局限性。因此，在实际应用中，最好是多个软件一起使用，取各自的长处，弥补相应的不足之处，才能达到满意的效果。

（4）半变异函数计算方法的选择　实验半变异函数的计算方法除了常用公式外，还有其他很多种计算方法如madogram公式法、rodogram公式法、ralative virogram公式计算法等。这些公式只有在特定的条件下才适用，如原始数据存在比例效应或为了消除原始数据中大值的影响等情况。但由于这些计算方法都不是经典的计算方法，在应用时应特别小心，尽量不用。

3.2.6.6　数据插值方法的选择

（1）根据研究目的选择插值方法　如果要对研究区未采样点进行最优无偏线性估计，则线性克里金插值技术是首选方法。如果要研究区域化变量的波动性及结构性，克服数据估值中的平滑效应，保持未采样点的空间结构性，则条件模拟值法是首选方法。如果要研究区域化变量超过某一阈值的概率，计算子域平均估计值，则非参数克里金法是首选方法。

（2）根据各种方法自身的特点选择插值方法

①一般克里金法。对于一般克里金法而言，它要求数据完全符合二阶平稳假设，无全局及局部趋势，区域化变量平均值为已知常数，区域化变量为正态分布或近似正态分布随着计算机技术的起步，数据变换方法逐步被大型软件采用，如新版ArcGIS软件，默认的数据插值方法已经由普通克里金法转为一般克里金法。这一条件对于实际研究显然很难满足，因此，一般克里金插值技术在实际应用中不多见。

②普通克里金法。普通克里金法也要求数据符合二阶平稳假设，但根据普通克里金的计算方法可知，普通克里金要求区域变化量的平均值为m，但m可以是一未知数（也就像泛克里金一样，只要知道漂移的形式，而不必求出漂移就能估计未采样点数值一样），又因为普通克里金法采用移动邻域的方法确定估值点和计算所需样点，这种算法本身就蕴含处理数据不平稳性的特性。因此，普通克里金可用于变化的平均值而平稳的协方差函数的情况，可称为准二阶平稳或准本征的情况。也就是说，如果区域化变量Z（x）在变程α之内有足够多的数据点，且半变异函数平稳，而变程α之外存在漂移，则半变异函数拟合时可只考虑变程α之内的半变异函数值，变程α之外的漂移可忽略，此时插值方法依然可有用普通克里金插值技术。但如果区域化变量全局平稳，而在局部小的范围内存在漂移，此时也不能用普通克里金法，此时计算的半变异函数会存在波动性，如图3-12所示。

普通克里金法要求监测数据符合正态性分布或近似正态分布，如果监测数据为偏态分布，则在使用此算法前应采用对数变换、平方根变换等方法使数据符合或近似符合正态分布。

③泛克里金方法。泛克里金方法应用于区域化变量不符合平稳假设存在漂移的情况，表现在半变异函数上为半变异函数值一直处于上升趋势，超过先验方差后依然在增加，而没有基台值，如图3-13所示。或半变异函数值存在较大的波动性。泛克里金的漂移一般用到1次或2次漂移，再高次的漂移一般不用。区域化变量的平稳与不平稳是一个相对概念，它与研究尺度有关，如有的区域化变量在大尺度时表现为平稳，而在小尺度则不平稳，有的区域化变量则相反（王家华，1999）。泛克里金要求监测数据符合正态或近似正态分布。

图3-12　整体平稳局部漂移的一维随机函数

图3-13　存在全局趋势的半变异函数

注：图中标注为数据样点量。

④指示克里金法。指示克里金法用于评估区域化变量超过某一或一组阈值的概率时的计算，本方法要求区域化变量符合平稳假设，但不要求数据正态性，因此，如果研究区某区域化变量不符合正态分布，或异常值较多而又不能剔除时，可用此方法。指示克里金方法可用于离散变量的估计。

指示克里金方法在计算过程有可能存在有序性问题，本问题可通过以下方法解决：一是把对应于所有门槛值的克里金方程组合并成一个大系统，使得估计方差的和为最小；二是取指示克里金算法所给出的结果，拟合一个分布函数，使得它与最优指示克里金解的差值平方的加权和最小；三是顺序的矫正有续关系。

⑤析取克里金法。析取克里金法适用条件与普通克里金法相似，但本方法要求监测值严格符合正态分布。因此，在应用本方法时，应首先采用非线性方法对原始数据进行转化。析取克里金方法在原理上较普通克里金法先进，理论上讲其计算结果较普通克里金更精确。但析取克里金法计算过程复杂，计算过程中的数据的取舍及近似过程往往使计算结果误差增大。因此，在进行克里金估值时选择析取克里金法还是普通克里金法应根据数据计算结果确定并把握简单化的原则。

⑥条件模拟值法。条件模拟值计算方法可分为高斯算法和指示算法两类。高斯条件模拟值法能更好地再现监测数据中极值的离散性，而指示模拟算法则更忠实于极值的空间变异性。如果区域化变量为二值变量，则应采用指示模拟值算法。如果区域化变量为连续性变量，模拟过程中不强调极值的空间变异性，则可用高斯模拟算法，但如果模拟过程中强调极值空间变异性的可重现性，则应将连续变量进行编码处理。高斯条件模拟值法、矩阵分解法适用于模拟网格节点数不大于1 000而要求的模拟次数较高的情况。序贯高斯条件模拟值法没有模拟节点数量的限制，一般情况下都可应用，但缺点是其计算速度相对较慢。

3.2.6.7　插值结果精度评估

克里金插值技术优点之一是可以采用交叉验证法对插值结果的优劣性进行评估，从而评估实验半变异函数及拟合模型、插值方法选择的优劣。一般认为计算的实验半变异函数及其拟合模型，选择的插值方法应使克里金估计值产生以下较好的统计结果：平均误差尽可能接近于0；均方差尽可能小；平均克里金方差尽可能小；标准克里金方差尽可能接近于1；估计值与估计误差的相关系数尽可能小；估计值与实验值的相关系数尽可能的大。交叉验证是一个渐进和摸索过程，需要不断调整各模型的参数，并根据实验半异函数的各向异性及各向同性性质多次设置邻域搜索方向及邻域实验数，以达到尽可能最好。

条件模拟值结果精度的评估，相关研究还没有给出完美的解决方案，但国内外科研工作者进行了大量的相关研究工作。如张泽浦（1998）、王学军等（2003）通过对模拟结果基本统计信息的描述及通过模拟值产生的半变异函数与实验半变异函数的相似性分析了模拟结果的正确性。徐英等（2004）给出了选择最终数字模拟地图的3条原则等。由于以往的研究没有给出最终模拟结果选择的数学表达方法，本研究在前人研究的基础上，提出了选择最终模拟数字地图的相关系数法。其基本思路是根据条件模拟值法的基本要求，其模拟产生的最终结果应与实际监测的数据有最大的相似性，即模拟结果数据的离散趋势、中心化趋势、分布趋势、百分位值应与监测数据有最大的相似性。当模拟结果数据与实际监测数据的相关系数最大时，模拟结果也最能代表实际数据的分布状况。本文采用代表数据离散趋势的标准差、方差、最大值、最小值、极差，中心化趋势的平均值、中值、众数，分布趋势的峰度、偏度，第10至第90分位值等指标来表征数据的总体特征。如果某一模拟结果的指标项与实际监测结果的指标项相关系数最大，则本次模拟结果最能代表实际数据的分布情况（赵玉杰等，2011）。