本文介绍基于ArcMap软件，实现普通克里格、回归克里格方法的空间插值的具体操作。
目录
- 1 背景知识准备
- 2 回归克里格实现
- 2.1 采样点与环境变量提取
- 2.2 子集要素划分
- 2.3 异常值提取
- 2.4 土壤有机质含量经典统计学分析
- 2.5 回归方程求取
- 2.6 残差提取
- 2.7 残差普通克里格求解
- 2.8 土壤有机质含量回归克里格求解
- 2.9 回归克里格精度评定
- 2.10 回归克里格专题地图制作
- 3 普通克里格实现
- 3.1 普通克里格精度评定
- 3.2 普通克里格专题地图制作
- 4 两种插值方法对比
- 4.1 精度对比
- 4.2 插值结果对比
- 5 一些值得讨论的问题
- 5.1 范畴型变量求解
- 5.2 ArcMap崩溃
- 5.3 环境要素提取零值处理
- 5.4 相关性分析与回归方程结果对比
- 5.5 Box-Cox变换
- 5.6 回归克里格结果错误
- 5.7 回归克里格结果异常值
- 5.8 研究区域水体辨别
- 参考文献
1 背景知识准备
前期几篇博客分别基于地学计算的基本概念与相关操作，进行了详细的讲解：
地统计学的基本概念及公式详解：地学计算主要概念及公式推导全解；
MATLAB计算变异函数并绘制经验半方差图：空间数据变异函数计算与经验半方差图绘制。
全局多项式（趋势面）与IDW逆距离加权插值：MATLAB代码：空间数据全局多项式插值法（趋势面法）与逆距离加权（IDW）法插值与结果分析。
本文所涉及相关操作的原理与对应代码等，都在上述三篇博客中，本次就不再赘述。
结合以上相关基础知识与基本操作方法，本次我们就将通过回归克里格与普通克里格这两种方法，基于ArcMap、MATLAB、SPSS等软件，计算土壤空间属性插值的数值。其中，由于本文所用的土壤采样点空间数据集并不是我的，因此遗憾不能将这一数据一并提供给大家；但是依据本篇博客的思想与对操作步骤的详细解释，大家用自己手头的数据，可以将相关操作与分析过程加以完整重现。
本文所用的软件：AcrMap 10.2软件、SPSS 25软件、MATLAB R2018a软件、Microsoft Office Excel 2019软件。
本文所用的数据：研究区域的矢量地图与空间数据集，其中空间数据集包括湖北省荆门市沙洋县288个土壤采样点位置数据与土壤有机质含量数据，以及对应区域90米分辨率高程（DEM90）、长度坡度系数（Length-Slope Factor，LS）、水流强度指数（Stream Power Index，SPI）、径流指数（the Wetness Index，WTI）等13个环境变量[1, 2]。其中，本文将上述13种环境变量分别简写作CH、CV、HillShade、LS、m、QFD、SOS、SPI、DEM90、WTI、α、β、Ω。
2 回归克里格实现2.1 采样点与环境变量提取
本文中，我所使用的288个初始土壤采样点数据为矢量图层文件，后续对其分析时需要提取为具体列表的形式；另一方面，本文13个初始环境变量数据均为栅格图层，呈面状分布；为进行后续回归方程求解、残差计算等操作，我们同样需要将各环境变量在采样点上的数值加以提取。上述步骤均在AcrMap 10.2软件中实现。
在AcrMap软件中，首先将采样点矢量图层导入，并选择“系统工具箱”→“Conversion Tools”→“Excel”→“表转Excel”，配置输入表、输出路径等信息后点击“确定”。这一步骤亦可于环境变量提取后进行。
随后，利用“系统工具箱”→“Spatial Analyst Tools”→“提取分析”→“多值提取到点”模块，对13个环境变量分别进行采样点对应位置的数值提取。提取时发现，若不由图层列表选择环境变量，而是由资源管理器直接选择，可以进行较为方便的多选操作。
提取完毕后，可以在采样点图层的属性表中查看288个点的位置信息、土壤有机质含量及其各自分别对应的13个环境变量的数值。随后发现，提取出的环境变量数值具有某个点多为0的情况。这一问题随后通过调整操作得以解决，具体附于本文5.3部分。
选择“系统工具箱”→“Conversion Tools”→“Excel”→“表转Excel”，配置输入表、输出路径等信息后点击“确定”，即可将包含有13个环境变量信息的采样点属性表导出。
2.2 子集要素划分
由于后期我们需要借助验证数据集对回归克里格与普通克里格方法进行对比，因此需要划分子集要素。
选择“系统工具箱”→“Geostatistical Analyst Tools”→“工具”→“子集要素”，配置相关属性后点击确认。
其中，设置训练要素为全部采样点中随机选取占比80%的部分，设置测试要素为全部采样点中随机选取占比20%的部分。
可以分别用两种不同颜色表示训练要素对应的点与测试要素对应的点，如下所示。
后期将用上述训练要素与测试要素数据对回归克里格与普通克里格插值结果加以对比。
2.3 异常值提取
首先，将分布于环境变量图层外的一个采样点剔除。这一部分的具体说明附于本文5.3部分。
此外，由于288个土壤采样点对应有机质含量数据为实测数据，其在采样记录、实验室测试等过程中，可能具有一定误差，从而出现个别异常值；而异常值的存在会对后期空间插值效果产生较大影响。因此，结合本文开头所提及的博客以及数据的实际情况，我们选用“平均值加标准差”方法的“3S”方式，借助MATLAB软件对本次作业异常数据加以筛选、剔除。
下图即为通过MATLAB筛选出的训练要素中的异常值点号信息。异常值需要在Excel与ArcMap中同时剔除。
2.4 土壤有机质含量经典统计学分析
运用SPSS软件对剔除异常值之后的285个采样点有机质含量数据加以统计分析。
结果如表1所示。
由表1可得，上述采样点土壤有机质含量最小值为7.51g/kg，最大值为34.66 g/kg；其平均值与中位数，相差不大，说明有机质分布整体较为均匀；其变异系数为0.23，呈现出中等变异性。其偏度与峰度分别为0.024与-0.669，说明其分布陡缓程度与正态分布较为类似，但稍偏尖顶峰状态；其斜偏度较之正态分布具有一定差距，呈现左偏。
2.5 回归方程求取
数据预处理结束后，即可开始利用训练要素中各点有机质含量及其对应的13个环境变量数值，求取其间的回归方程。
首先需要将上述数据导入SPSS软件。
在SPSS软件中，选择“分析”→“相关”→“双变量”，利用皮尔逊相关系数在双尾检验条件下，对土壤有机质含量与其它13个环境变量之间的相关性加以探究。
相关性结果如表2所示。
随后，选择“分析”→“回归”→“线性”，将“方法”设置为“步进”，并分别设置进入和除去的概率为0.05、0.10，0.10、0.11与0.15、0.20。将上述三种概率对应的回归模型分别编号为1、2与3，并依次得到三种模型对应的回归方程、R^2、F与显著性。
对应回归模型及其参数如表3所示。
2.6 残差提取
得到上述三种回归模型后，分别利用已有环境变量栅格图层对土壤有机质含量以及计算。
在ArcMap软件中，利用“系统工具箱”→“Spatial Analyst Tools”→“地图代数”→“栅格计算器”模块，分别生成三种回归模型对应的土壤有机质含量图层。
求出上述三个土壤有机质含量图层后，利用“系统工具箱”→“Spatial Analyst Tools”→“提取分析”→“多值提取到点”模块，提取出训练要素中各点对应的有机质含量数值。
提取后导入到Excel软件，并利用Excel分别求出上述三种回归模型对应训练要素各点的残差。
利用MATLAB软件对三种回归模型对应残差进行正态分布检验，发现模型三残差原始数据符合正态分布，因此此处暂选择模型三进行后续计算。随后发现其所得结果出现异常值，且考虑到模型三变量显著性水平选取过高，因此最终借助模型二对应残差重新执行下述操作，求得最终有机质含量回归克里格结果。
2.7 残差普通克里格求解
将上述求得的训练要素各点对应残差数据导入ArcMap软件训练要素点图层中，利用“Geostatistical Analyst Tools”→“地统计向导”→“克里金法/协同克里金法”模块，对残差加以普通克里格插值。
试验变异函数散点图及其拟合如下。
配置相关属性，完成普通克里格插值。其中，由于残差数据自身已通过正态分布检验，因此未在上述步骤对数据加以变换。同时，依据计算得出试验变异函数散点图，选择球状模型对其加以拟合。
得到结果图层如下所示。
将沙洋县区域图层同时叠加，可以看到普通克里格插值结果并未完全覆盖沙洋县整个区域。
利用结果图层的“图层属性”→“范围”模块，将插值结果扩大至完全覆盖沙洋县。
确认无误后，利用结果图层的“数据”→“导出数据”模块，将残差克里格插值结果图层转为栅格面图层，并依据沙洋县边界范围将得到的栅格面图层加以剪裁。
剪裁后，即可得到沙洋县土壤有机质含量残差的普通克里格插值完整结果图层。
2.8 土壤有机质含量回归克里格求解
在ArcMap软件中，利用“系统工具箱”→“Spatial Analyst Tools”→“地图代数”→“栅格计算器”模块，将上述残差克里格插值结果图层与前述模型三所对应回归方程计算结果图层相加，从而得到利用回归克里格插值方法计算的沙洋县土壤有机质含量插值结果图层。
得到结果如下所示。
可以看到其中具有-140左右的异常值；尝试运用前述模型二进行同样的求解操作，发现异常值不再存在。针对这一问题的思考与结论，具体附于本文5.7部分。
2.9 回归克里格精度评定
随后，利用“系统工具箱”→“Spatial Analyst Tools”→“提取分析”→“多值提取到点”模块，将土壤有机质含量的回归克里格插值结果提取至测试要素对应点位置。
利用MATLAB软件，求取测试要素中各点土壤有机质含量实测值与回归克里格计算结果之间精度衡量指标。
各精度衡量指标结果如表4所示。
2.10 回归克里格专题地图制作
首先，将计算过程中的残差普通克里格插值结果、回归方程计算结果等制作为专题地图。其中，由于暂未找到相关精确资源，因此下列展示结果未对沙洋县境内的水体加以区分。这一问题亦在本文5.8部分有所讨论。
最后，将基于回归克里格方法得出的沙洋县土壤有机质含量插值结果制作为专题地图。
3 普通克里格实现
基于与前述内容类似的步骤，直接将经过异常值剔除的训练要素各点对应的土壤有机质含量进行普通克里格插值。
其中，由于土壤有机质含量均为正数，因此在进行克里格插值时采取了Box-Cox变换。
试验变异函数散点图及其拟合如下。
3.1 普通克里格精度评定
利用测试要素对普通克里格插值结果加以精度衡量，所得各项精度衡量指标如表5所示。
3.2 普通克里格专题地图制作
将基于普通克里格方法得出的沙洋县土壤有机质含量插值结果制作为专题地图。
4 两种插值方法对比
结合上述回归克里格与普通克里格所得结果，由精度与实际结果等角度对两种方法加以对比、分析。
4.1 精度对比
将前述回归克里格插值与普通克里格插值结果精度衡量指标加以对比，探究两种方法的效果。
指标对比如表6所示。
由表6可知，由整体角度观之，普通克里格插值方法与回归克里格方法对应平均误差、平均绝对误差与均方根误差均较小，即二者对于土壤有机质含量的插值计算结果较为准确。
由两种方法结果对应精度衡量指标观之，普通克里格方法对应平均误差较小于回归克里格方法；而普通克里格方法平均绝对误差与均方根误差均较大于回归克里格方法，且其相关系数小于后者。
其中，由于平均误差在计算时未对误差求得绝对值[3]，因此其符号可以反应计算结果与实测值的大小关系，其数值大小则表明所得结果的可靠程度；而平均绝对误差和均方根误差均可以分别较好反应出计算值与实测值之间误差的实际情况与残差的样本标准差；相关系数则可以评价计算值与实测值之间的线性相关程度。
由此可以看出，两种克里格方法求得的土壤有机质含量数据均较大于实测数据，但二者计算结果的整体准确性较高；其中，回归克里格方法求得的测试要素点土壤有机质含量较之普通克里格方法，平均绝对误差降低3.27%左右，均方根误差降低5.97%左右。
综上所述，回归克里格方法相对普通克里格方法所得计算结果更加准确，在一定程度上提高了空间插值的效果。
4.2 插值结果对比
借助上述普通克里格与回归克里格插值方法所得结果图，对其各自特点加以对比、分析。
此处值得注意的是，本文前述两幅克里格插值结果图层中，我均采用了“拉伸”色带式图例。而为更好对二者插值结果图加以对比，此处将二者结果图层的图例加以统一处理，运用相同的图例表示方式，即使得两幅图中同一颜色代表相同的数值含义。
对比上述两幅图。由空间区位分布角度观之，普通克里格方法与回归克里格方法所得插值结果整体趋势一致，呈现出土壤有机质含量自沙洋县中、西部地区向东部地区递减的变化特征。即在沙洋县中部、西北部，包括南部地区，土壤有机质含量整体较高，而在东部地区含量整体较低；其中，含量最高区域主要分布于沙洋县中部夹堰村一带与西部常家湾一带，含量最低地区则主要分布于沙洋县东部新村、登家台一带。
由空间聚集分布角度观之，普通克里格方法所得插值结果整体较为平缓，多呈现块状分布，形成较类似于等高线状的局部极大值或极小值中心分布趋势；而回归克里格方法所得插值结果较之前者更加分散、零碎，空间变异较为复杂，很少出现结果聚集分布的态势；很好展现出一些局部空间、细节地区的土壤有机质含量情况。
此外，可以看到回归克里格方法插值结果可以明显呈现出沙洋县河流的分布区位，以及其南部地区明显可以看到长湖的分布位置；而在普通克里格方法插值结果中则看不出这些细节。之所以出现这样的原因，我认为是因为回归克里格中起到主要决定作用的趋势项造成的。如本文第二部分所述，在回归克里格方法计算过程中，由于趋势项在通过回归方程求解时，融入了大量环境要素（其中就包括与河流关系较为密切的SPI），因此可以较为明显的看出上述河流、湖泊等环境要素的区位信息。
此外，由于上述两幅图使用同一图例，或暂时看不出普通克里格方法与回归克里格方法插值结果的数值范围差异。而若结合前述两幅专题地图，则可以看出普通克里格的插值结果区间为9.5339至30.9705g/kg，而回归克里格的插值结果区间为6.6434至34.0612g/kg；即回归克里格插值结果范围区间较之普通克里格大。
这一结果亦符合相关文献的研究成果[3]。这一现象的原因为，普通克里格插值由于更多依赖于临近点的实测值，通过变异函数等加以待插点数值的计算，因而各点所受临近采样点的影响较大，从而导致各点插值结果趋于平缓，从而进一步“磨灭”了数据的极值，导致最大值被过低估计，而最小值被过高估计；因此产生这一现象。
5 一些值得讨论的问题
将一些在上述操作过程中遇到的，或者我所想到的问题，在这里统一记录与思考。
5.1 范畴型变量求解
本文前述各环境变量均为连续型变量，相关操作均为对具体的连续型数值信息加以处理，属于“硬插值”[1, 7]；这种方法在一定程度上忽视了非数值类型变量对土壤相关属性的影响，如土壤类型、土壤质地等[7]。因此，结合考虑了数据不确定性的“软数据”，将非数值类型变量纳入计算范围，可以较为有效地提高计算精度[7]。参考相关资料，结合数学建模相关知识，尝试由研究方法与实际处理两个角度针对范畴型变量展开讨论。
在研究方法层面，一是可通过变换研究方法方式对其加以处理。有学者通过贝叶斯最大熵（Bayesian Maximum Entropy，BME）方法，将相关环境属性数据预测分为先验阶段与后验阶段两个部分。其中，先验阶段依据经典统计学指标、自然法则等一般知识数据计算先验概率密度函数；而在后验阶段，将“硬数据”与“软数据”结合，进行贝叶斯条件化，获取后验概率密度函数，并进而完成计算[7, 8]。
二是可对部分操作所选用方法加以调整。例如，在相关性分析层面，本文选择对13个相关环境变量与土壤有机质含量进行皮尔逊相关性分析；而针对范畴型变量，例如对两个无序分类变量进行相关性检验，则需要借助卡方检验（Chi-square Test），由理论频数和实际频数的吻合程度对变量的相关性加以定量分析。在回归分析层面，本文选择利用线性逐步回归对模型加以求解；而针对范畴型变量进行回归方程求解，如因变量为一个无序分类变量，自变量为多个无序分类变量或多个二分变量（多分类变量亦可，但二分变量运用较多）与连续变量结合，则可使用Logistic回归方法。
上述讨论在研究方法层面对包含范畴型变量的数据处理加以探究。而在更为具体的实际处理层面，可通过如下方式对范畴型变量加以计算。
一是可对范畴型变量自身加以变换。可选择直接将某一多分类变量对应各类分别转换为若干二分变量，进而直接利用上述转换后的二分变量加以相关处理与计算。
其中，以不同土壤质地类型的划分为例，探讨范畴型变量变换这一方法。文献[2]将土壤质地类型分为Light loam、Loam、Sandy loam、Sand、Fine sand与Medium loam等6类；则我们可分别使用x_1、x_2、x_3、x_4、x_5、x_6等6个二分变量对其加以表示：
其中，6个二分变量同时仅可有一个为1，其余均为0。
随后，可将回归方程写作：
其中，a、b、c、d、e、f均为系数，α为其它环境变量对应项。通过这一公式，即可将范畴型变量进行定量化计算。
二是可建立范畴型变量与连续型数据变量之间的关系。文献[7]依据土壤有机质含量将训练集中各点分为n组，并分别对各组统计组内属于各土壤质地类别的点数，并将各点数除以该组点数总和，作可作为有机质与土壤质地之间的定量关系。
其中，R_(〖OM〗_k )表示第k个有机质组对于土壤质地的定量关系，〖ZD〗1表示第1个土壤质地类别，〖cou〗(1,k)表示第k个有机质组中属于第1个土壤质地类别的点个数，〖cou〗_k表示第k个有机质组中全部点个数，m表示土壤质地类别总数。
随后，依据待预测点Z_ij所属土壤质地类别，依据上式反求出该点有机质取值概率，即其对第k个有机质组的相似度。
将该点对各有机质组的相似度归一化，则可获得该点对应有机质含量的模糊分布。至此，即将上述土壤质地这一范畴型变量转换为具有一定不确定性的“软数据”，并进一步进行后续计算。
5.2 ArcMap崩溃
本文操作部分初期，若直接将原始288个土壤采样点导入个人电脑ArcMap 10.2软件中，会出现软件强制停止运行的现象；多次尝试无果。多台电脑尝试后，发现部分电脑同样具有这一问题。
随后，通过多次尝试发现，若在导入这一初始采样点图层前，任意添加一幅“.png”格式图片，软件则将不会报错，继续正常运行。
结合后续操作，个人认为这一问题可能是原始采样点数据文件中部分内容与ArcMap 10.2软件不兼容导致。
5.3 环境要素提取零值处理
在对环境要素进行“多值提取至点”时，我原本将“点位置处的双线性插值（可选）”选项选中。但这样得到的提取结果中会出现个别点多数环境变量数值均为0的情况。
因此，尝试在ArcMap软件中将提取后环境变量多为0的点（共四个）一同选中，查看其分布位置。可以发现，上述出现数值多0的点均分布于沙洋县边界区域。
即使用插值后，边界点受到沙洋县范围外空白像元的影响，从而导致提取的环境变量为0。因此，取消选中上述插值功能，并再次执行“多值提取至点”操作；但随后发现，这样执行操作后依然仍有一点的各项环境变量数值提取结果多为0。
同样在ArcMap软件中选中这一点，并放大。
可以看到，该点位于沙洋县边界图层范围内，而不位于环境变量栅格图层范围内——即处于二者之间。因此，该点在提取时无论是否使用插值方法，均会导致最终的提取结果多为0。
由此，我选择将该点作为异常值点，并在后续异常值剔除步骤中一并将其删去。
5.4 相关性分析与回归方程结果对比
在执行逐步线性回归前，我先利用皮尔逊相关性分析方法，对13个环境变量数据与土壤有机质含量之间的相关性加以衡量；随后，利用逐步回归方法求出环境变量与土壤有机质含量的回归方程。具体相关性分析结果与回归方程结果已附于本文前述表2、表3处。
由两表所示结果可知，相关性分析所得结果与回归方程所得结果具有一定不相符之处。例如，相关性分析显示，与土壤有机质含量具有较显著相关关系的环境变量包括QFD、SOS、SPI、DEM90与β等，其余环境要素与有机质含量之间相关性不显著；而回归模型中，却包含了部分相关性不显著的环境变量，如模型三中包括了m、α与WTI等。
随后通过互联网等资源获知，相关性分析仅针对某单一自变量进行其与因变量之间相关性的分析，即各环境变量在进行相关性分析时为独立状态；而逐步回归则是将模型内已有全部变量结合，作为整体来衡量模型与因变量之间的拟合关系。因此，针对二者不同的结果，应以逐步线性回归方法最终纳入回归模型中的环境变量为准。
5.5 Box-Cox变换
本文中，涉及到残差的普通克里格与土壤有机质含量的普通克里格等。其中，后者由于均为正数，因此可适合大部分正态分布转换方法；而前者由于存在部分负值，从而受阻。
其中，在针对不同模型对应残差进行普通克里格计算时，考虑到Box-Cox变换依据幂参数选取的不同，具有两种形式；因此误认为其可以适用于存在负值的数据。而随后的实践中发现这一想法错误。
因此，借助MATLAB软件“boxcox”函数的官方说明，进一步确定Box-Cox变换仅适用于正数情况。
5.6 回归克里格结果错误
在进行回归克里格计算时，前期计算结果总和所得普通克里格土壤有机质含量插值结果有着较大差异。经过对比，发现所得回归克里格结果在整体趋势上与普通克里格所得结果相反——回归克里格插值结果较大区域，其普通克里格结果反而较小；反之亦然。
上图即为出现问题的回归克里格插值结果。可以看到，这一结果与个人最终结果所呈现的有机质含量空间分布趋势正好相反。
因此对这一问题进行多次重复操作等处理，但结果均表现出同样的错误。随后，考虑到所得结果与普通克里格结果正好为相反趋势，因此对残差加以检查；从而发现由于残差在计算时，错将实测值与回归方程值之差误写作回归方程值与实测值之差。将这一问题改正后重新执行回归克里格插值，得到与普通克里格一致的结果。
得到上述结果后进一步思考，进而引出又一问题——如前所述，一般地，往往认为回归克里格中，对结果起到较为重要作用的因素为空间趋势项，而非残差项。而为何在此处，由于残差计算的错误导致最终回归克里格结果出现如此之大的差异？
随后，个人结合相关中、外文献与本文前述土壤有机质含量经典统计学分析结果（表1），尝试探究并发现问题所在。首先，回归克里格中起到重要影响作用的成分的确是空间趋势项（这即是其较之普通克里格方法结果更加零碎、更好展示出细节的原因），而其影响具体有“多么”重要，则需结合研究区域实际情况等加以判，其随着目标变量的空间变异性特征而发生变化。其次，上述错误相当于将残差符号取反，即造成的影响为残差原有比重的两倍，进而对最终结果造成较大影响。
5.7 回归克里格结果异常值
利用前述回归模型三进行回归克里格插值，所得结果中包含-140左右的负值。结合ArcMap软件“识别”工具，对地图中较低值区域加以数值获取，发现其数值均为15左右，并未找出具体-140值所在。
由此，针对这一异常值问题，个人判断是由于回归方程纳入了一些具有较低或较高值的环境变量，从而导致极个别回归模型结果图层像素出现异常值。由于该值较低，导致即使加入残差项后亦并无明显改进。随后，尝试利用模型二对应回归方程进行后续操作，问题得以解决。
另一方面，前述模型三对应逐步回归过程中，所选用进入和除去的概率均较大，故其部分环境变量系数显著性水平偏高（即数值较大）；选用模型二在避免异常值像素出现的同时，同样避免了这一显著性水平问题。
5.8 研究区域水体辨别
在本文中我们可以注意到，所研究区域沙洋县境内水体较多——如其东北部边界紧邻汉江，南部区域为大面积的长湖，且境内包含金鸡水库、安洼水库、潘集水库等小型水库，以及各类小型河流；因此，本次在利用普通克里格方法或回归克里格方法对研究区域土壤有机质含量进行插值计算时，水体部分自然需要剔除。
一般地，为提取水体等区域，需要结合数字化后的、精细的土地利用现状图或相关遥感卫星影像等。而在本次操作中，个人暂未找到方便获取、质量较好，且可在相关地学软件中处理的资料等。因此，如前所述，在运用两种克里格方法进行插值计算时，我并未将水体剔除。
在本文操作部分进入尾声、进行专题地图制作时，我曾分别考虑利用相关图片编辑软件对水体中面积较大的长湖进行区分或利用实验材料中NDVI数据提取水体，分别如下左、右图所示。
上述第一种方法，一方面由于水体为手动编辑，十分不严谨，不符合科研的严密性，使得实验结果可靠性严重下降；另一方面，其仍未对沙洋县其它水体加以区分。上述第二种方法，一方面水体与云层、雪等NDVI数值均较小，在不了解对应遥感图像成像气候条件情况下提取结果可靠性不高；另一方面，提取过程中所选NDVI阈值具有主观性，阈值对应结果因而均有不确定性。因此，最终结果中放弃了对水体进行提取操作。而在其它可获得相关资料的情况下，还需尽可能避免上述将水体纳入土壤属性插值范围的情况。
此外，可以借助遥感影像，结合监督分类等方式对水体区域加以提取。监督分类的相关操作大家可以参考ENVI+ERDAS实现Hyperion叶绿素含量反演：经验比值法、一阶微分法。
参考文献
[1]Zhang C, Yang Y. Can the spatial prediction of soil organic matter be improved by incorporating multiple regression confidence intervals as soft data into BME method?[J]. CATENA. 2019, 178: 322-334.
[2]Zhang S, Huang Y, Shen C, et al. Spatial prediction of soil organic matter using terrain indices and categorical variables as auxiliary information[J]. Geoderma. 2012, 171-172: 35-43.
[3]曹祥会，龙怀玉，周脚根，等. 河北省表层土壤有机碳和全氮空间变异特征性及影响因子分析[J]. 植物营养与肥料学报. 2016, 22(04): 937-948.
[4]李龙，周飞，田杰，等. 地形因子对半干旱地区土壤有机碳含量的影响[J]. 北方园艺. 2019(16): 104-109.
[5]张国峰，杨立荣，瞿明凯，等. 基于地理加权回归克里格的日平均气温插值[J]. 应用生态学报. 2015, 26(05): 1531-1536.
[6]姜勇，梁文举，李琪. 利用与回归模型相结合的克里格方法对农田土壤有机碳的估值及制图[J]. 水土保持学报. 2005(05): 99-102.
[7]张若兮，杨勇，张楚天. 基于范畴型变量和贝叶斯最大熵的土壤有机质空间预测[J]. 土壤通报. 2015, 46(02): 312-318.
[8]李爱华，柏延臣. 基于贝叶斯最大熵的甘肃省多年平均降水空间化研究[J]. 中国沙漠. 2012, 32(05): 1408-1416.