做空间插值的同学应该都有这种体会:同一份点数据,普通克里金偶尔会“翻车”——数据明明在某个方向有持续的抬升或者递减趋势,插出来的预测图却总是一块一块的,残差还带有明显的空间结构。这时候就该考虑泛克里金插值了。泛克里金(Universal Kriging)并不是一个新的“魔法算法”,它做的事情其实很朴素:把数据中隐含的区域性趋势当作可估计的成分先处理掉,再对剩下的平稳残差做克里金插值,最后把趋势面加回来。它解决的核心问题,就是当数据不满足普通克里金要求的平稳性假设时,如何处理那些缓慢变化的主导趋势。
这篇内容我按自己平时实际操作的顺序来写:先说为什么非要用泛克里金,再讲必须搞懂的几个核心概念,然后是 ArcGIS 里的完整操作流程,接着是参数调优和常见坑,最后聊一下它和其他插值方法的选型。适合 GIS 相关专业的学生,以及做土壤、气象、水文、环境调查这些需要把离散点转成连续栅格的从业者,内容包括可以直接复现的步骤和判断标准。
1. 为什么要用泛克里金:从普通克里金的局限说起
1.1 普通克里金的“平稳性”假设到底卡在哪里
普通克里金(Ordinary Kriging)是地统计插值里最常用的方法,它的核心假设是数据满足二阶平稳或内蕴平稳:空间上均值恒定,变量之间的相关性只取决于两点之间的距离而非绝对位置。但在实际工作里,这条假设经常不成立。举几个我做过的例子:一个采样区的地下水埋深沿河流方向逐渐变浅,这是确定性趋势;土壤有机质含量从坡顶到坡脚呈递减变化;城市热岛效应下气温从中心向外围均匀下降。这些数据都带有明显的“趋势项”,均值和位置强相关,直接上普通克里金,半变异函数会被趋势污染,算出来的块金值和基台值整体偏大,插值结果在趋势方向会出现系统性偏差。
泛克里金的思路就是把这层趋势显式地建模出来,而不是假装它不存在。它将区域化变量表示成两部分之和:一个确定性趋势项 m(s),一个平稳随机残差 ε(s)。趋势项用坐标的多项式来拟合,比如一阶趋势就是平面方程,二阶趋势就是二次曲面方程。先拟合趋势面,再从原始值中减去趋势,得到平稳残差,对残差做普通克里金插值,最后把趋势面加回来,得到最终预测面。这其实就是“回归 + 残差克里金”的组合,只是 ArcGIS 把它封装成了一个单独的地统计模型。
1.2 趋势项存在时用普通克里金会发生什么
我用一组模拟数据验证过这个问题。在数据里人为设置一个沿东西方向的线性趋势,再用普通克里金和泛克里金分别做插值。普通克里金的预测误差均值出现明显偏移,预测标准差图在趋势方向呈条带状拉长,交叉验证的标准化均方根误差明显偏离 1。原因倒也好理解:半变异函数是对所有方向上的差异取平均,趋势项会造成相邻点在趋势方向上的差异被系统性高估,于是模型把结构性的趋势误判成了随机性的空间变异,块金效应被抬高,空间自相关的真实强度被掩盖。
相比之下,泛克里金把趋势先拿掉,残差的空间自相关结构会干净很多。所以当你发现数据存在明显的区域趋势时,泛克里金是比普通克里金更合理的选择。但这也引出另一个问题:怎么判断趋势是否存在?我在第 4 节会专门讲,这里先记住一个结论——如果直方图正常、分布也不错,但交叉验证指标一直调不好,先把“是否存在趋势”这个问题查清楚。
2. 上手前必须搞懂的核心概念
2.1 平稳性假设与趋势/漂移
平稳性在地统计里分几个层次。严格平稳要求整个空间上的联合分布不随位置平移而改变,这在实际数据里几乎不可能成立;二阶平稳只要求任意两点的协方差只与距离有关,这是普通克里金的理论基础;再弱一点是内蕴平稳,只要求增量 z(s) - z(s + h) 的方差存在且只与距离有关。
泛克里金放松的是“均值恒定”这一条:它允许均值是空间位置的函数。这个均值函数就是“趋势”或“漂移”。ArcGIS 的 Geostatistical Wizard 里通过 Trend Removal Order(趋势移除阶次)来设置:一阶趋势对应线性平面,二阶趋势对应二次曲面。需要注意,泛克里金里的趋势项通常只做低阶多项式近似,不是用来拟合并预测复杂变异细节的——趋势拉得太高,残差会被掏空,反而得不偿失。
2.2 半变异函数:空间自相关的核心载体
泛克里金的残差部分仍然靠半变异函数来建模。半变异函数表达的是空间上两点差异的方差随距离的变化,定义为:
γ(h) = (1 / (2N(h))) · Σ [z(sᵢ) − z(sᵢ + h)]²
其中 h 是两点之间的距离,N(h) 是距离为 h 的点对数量。实际计算时,样本点不可能正好在每个距离上都有点对,所以要把距离分组,按滞后距离(Lag Size)划分成若干个区间,每个区间内计算平均半方差。ArcGIS 会用一个拟合曲线去近似这些经验点,这个拟合曲线就是理论半变异函数模型。
理论模型里最核心的四个参数:块金效应(Nugget)表示距离趋近于 0 时的空间方差,主要来源于测量误差和微尺度变异;基台值(Sill)表示空间方差趋于稳定的上限;变程(Range)表示空间自相关存在的最大距离,超过变程就可以认为两点之间互不影响;偏基台值(Partial Sill)等于基台值减去块金值,反映真实空间结构贡献的方差比例。
半变异函数模型本身有很多种:球面模型在原点附近线性上升,到变程后趋于平缓,适合土壤属性这类空间连续性适中的变量;指数模型从原点弹性上升后渐近逼近基台值,适合温度、降水这类渐变较慢的气象变量;高斯模型在原点附近非常平缓,适合高度连续、光滑的现象,但对采样误差很敏感。实际选哪个模型,不要只看哪个拟合优度高,交叉验证的误差小才是硬指标。
2.3 泛克里金与普通克里金的本质区别
普通克里金在估计权重时,只要求权重之和等于 1,以此保证无偏性;泛克里金除了要求权重之和为 1 之外,还要求权重与趋势项的多项式基函数正交。这在实际计算里意味着:普通克里金把区域均值当作未知常量,泛克里金把均值当作坐标多项式的函数来估计。最直观的差别是,普通克里金的预测结果在数据范围外会迅速回归到区域平均值,而泛克里金在数据范围外还会延续趋势方向的变化趋势。
ArcGIS 的 Geostatistical Wizard 里,Kriging Method 下拉菜单中可以直接选择 Universal Kriging,这是一种显式的泛克里金实现;同时 Ordinary Kriging 下面也提供了 Trend Removal 的选项。两者在做趋势移除时的数学处理可以说互为表里,但 Universal Kriging 把趋势项与克里金权重放进同一个方程组中联合求解,理论上更严谨。我个人的习惯是:明确判断出趋势后用 Universal Kriging 更省心,如果只是想测试性对比一下,用带趋势移除的 Ordinary Kriging 也能跑,只是后续解释起来要格外注意。
3. ArcGIS 泛克里金插值完整流程实操
3.1 数据准备与预处理
泛克里金虽然是地统计模型,但操作前要先检查数据质量,这一步决定后续所有结果的可靠性。
第一,确认数据是点要素且包含数值型字段。面要素和多部件要素不能用;只有文本字段也不能用。
第二,确认数据使用投影坐标系。克里金方法依赖真实距离计算,如果你拿到的采样点还是经纬度,必须先用 Project 工具投影到合适的投影坐标系。否则半变异函数的步长和变程实际上使用的是“度”而不是“米”,结果完全没意义。这里强烈建议按数据所在纬度选择合适的投影,不要偷懒用 Web Mercator。
第三,检查重复点和异常值。多个完全重合的样点会让半变异函数在距离接近 0 处的点对数量暴增,块金效应被严重高估;异常值对半变异函数的影响更是灾难性的。我习惯先用直方图和 QQ Plot 检查数据分布,再看 Voronoi 图找离群点。样本量方面,个人经验是至少要有 30 个有效点,少于这个数趋势项拟合和变异函数估计都会不稳。
3.2 打开地统计向导并配置模型
ArcMap 和 ArcGIS Pro 的入口稍有不同,但核心流程一致。在 ArcMap 里先确认自定义 > 扩展模块中的 Geostatistical Analyst 已勾选启用,然后从 Geostatistical Analyst 工具栏下拉箭头选择 Geostatistical Wizard。Pro 里则直接在分析工具箱搜索栏输入“Geostatistical Wizard”即可。
在向导的第一个面板中:
- 左侧方法列表选择 Kriging / Co-Kriging。
- 数据集选择你的点要素图层,目标字段选择要插值的数值字段。
- 在 Kriging Method(类型)下拉菜单里选择 Universal Kriging。
- Trend Removal Order 选择 First(一阶趋势)或 Second(二阶趋势)。
这里要强调一点:Trend Removal Order 不是越高越好。一阶趋势只能描述线性倾斜面,二阶趋势可以描述弯曲面,但阶次太高容易把真实的空间变异也吸收进趋势项,残差变成纯噪声。通常先做一次趋势分析工具,看数据的趋势面形态再决定用一阶还是二阶,后面我会给具体判断方法。
然后点击下一步进入 Semivariogram/Covariance Modeling 面板。这一步是泛克里金插值的核心,也是最需要人工介入的一步。ArcGIS 默认会帮你自动拟合模型参数,但自动拟合往往基于全局最优化,不一定适合你的数据,建议勾选手动调整后再点“Optimize model”看自动拟合结果,再根据交叉验证情况微调。
3.3 关键参数设置与交叉验证
在半变异函数与协方差建模面板中,需要确认四项内容:
滞后大小(Lag Size):ArcGIS 默认用最大距离除以步数来自动计算。如果你觉得云图上的点太稀疏或者太拥挤,可以手动调整。经验上,滞后大小取样点间平均距离的 1/2 到 1 倍比较合适,过大损失细节,过小后期点对太少曲线不稳定。
步数(Number of Lags):默认 12 步左右。总步数乘以滞后大小基本覆盖了点对距离分布的绝大部分范围。如果变程远超最大滞后距离,说明步数太少或滞后太小。要调整到半变异函数曲线呈现“先上升后平稳”的形态,这样后续拟合才有意义。
半变异函数模型:在下拉列表里换着试。不要只点了一个模型就结束,我会把球面、指数、高斯这几个常用模型都跑一遍,最后看交叉验证指标再定。
各向异性(Anisotropy):如果四个方向上的半变异函数曲线差异明显,考虑勾选各向异性并设置旋转角度和缩放系数。处理地下水埋深这类具有明确水流方向的数据时,这个选项几乎是必开的。
完成参数设置后,点击“Cross Validation”标签页查看交叉验证结果。这里有几个核心指标:
| 指标 | 理想值 | 说明 |
|---|---|---|
| 预测误差均值(Mean Predicted Error) | 接近 0 | 偏差的总体大小,过正则存在系统偏差 |
| 均方根误差(Root Mean Square Error) | 尽可能小 | 预测值与实测值的平均偏差幅度 |
| 平均标准误差(Average Standard Error) | 与 RMSE 接近 | 模型预测不确定性的平均值 |
| 标准化均方根误差(Root Mean Square Standardized Error) | 接近 1 | 大于 1 说明低估预测不确定性,小于 1 说明过度平滑 |
| 标准化平均误差(Mean Standardized Error) | 接近 0 | 标准化后的系统偏差 |
我通常优先看标准化均方根误差和预测误差均值。标准化均方根误差落在 0.9~1.1 之间都算可以接受,落出去的话就回退调整模型参数。
3.4 输出预测栅格与报表解读
交叉验证满意后,在向导里点击完成,ArcGIS 会在内容列表生成一个地统计图层。在内容列表中右键这个图层,选择 Data > Export Raster 可以输出预测栅格。导出时注意点:
- 输出像元大小(Cell Size)不要贪小,设得太小不仅文件大,而且在小尺度上预测结果会因为搜索邻域内样本不足而出现大量无效像元;我一般直接取采样点平均间距的 1/3 到 1/2。
- 输出范围默认跟随数据范围,如果希望贴合研究区边界,可用扩展数据范围中的 Clip 工具先切好掩膜,再在导出栅格时把分析范围设为图层范围。
- 别忘记同时导出预测标准差栅格。只有预测图没有不确定性图,在正式报告里说服力会差很多,审稿人通常都会问这一句。
4. 参数选择与模型优化的实战经验
4.1 如何判断数据是否存在趋势
判断趋势存在与否,我一般用两个办法,两个结果互相印证。
第一个是使用地统计模块里的“趋势分析”工具。在 ArcMap 的 Geostatistical Analyst 工具栏下拉箭头里选择 Explore Data > Trend Analysis,在弹窗里可以通过旋转视角观察:如果绿线(代表南北方向)或蓝线(代表东西方向)出现明显的高低起伏、倾斜或弯曲,就说明对应方向存在趋势。如果所有方向的线都基本平直,则可以优先考虑普通克里金。
第二个办法更定量一些:把点要素转成表格,用字段计算或脚本导出 x、y 坐标和 z 值,在 Excel 或统计软件里做 z 对 x、y 的多元线性回归,看 x、y 回归系数是否显著,R² 高不高。如果回归显著且 R² 大于 0.3,数据中的趋势成分已经不小了,用泛克里金是合理的。
另外还有一个辅助判断指标:如果普通克里金的预测误差图呈现明显的“香蕉形”或“条带状”分布,那通常意味着残差里还有结构性的趋势信息没有被提取干净,这时候泛克里金往往能改善结果。
4.2 半变异函数模型选型的取舍
半变异函数模型的选择没有绝对的对错,但有规律可循。球面模型是默认选项,本质上是空间相关性随距离线性衰减然后归零,适合大多数环境变量;指数模型衰减速度前期快后期慢,适合温度、降水这类尺度较大的气象变量;高斯模型在原点处曲线平滑,适合高程这类变化极为平缓的连续表面。问题在于,真实数据往往没有这么理想,比如土壤数据在短距离上可能因为采样误差出现明显的块金效应,这时候高斯模型会严重高估短距离的相关性。
我的做法是在保持趋势阶次和搜索邻域不变的前提下,把球面、指数、高斯、Matern 这几个模型各跑一遍,比较交叉验证的标准化均方根误差和均方根误差。有一次做土壤盐分数据,球面模型的标准化 RMSE 是 1.2,换成指数模型后降到 0.95,预测误差均值也从 -0.35 降到接近 0。所以不要迷信默认参数,花十分钟换模型试一遍,结果差异往往很明显。
4.3 搜索邻域设置对结果的影响
搜索邻域决定了预测点由哪些样本点参与计算。ArcGIS 里默认的 Standard 邻域会按扇区划分,默认是四扇区,并设置最大和最小邻居数。如果邻居数太少,预测结果容易围绕样本点出现“靶心效应”,也就是以样点为中心的圆圈状假象;如果邻居数太多,权重会被平均化,局部细节被平滑掉。
经验上,对于几百个点的数据,最小邻居数设为 2~5,最大邻居数设为 10~20 是比较稳的起点。如果样点分布极不均匀,可以考虑勾选 Use sector rotation 或增加扇区数量,确保每个方向都有足够的样本参与权重分配。注意一点:泛克里金的趋势项参数是在全局范围内估计的,但搜索邻域设置会直接影响局部预测的平滑程度,两者需要配合调整。
4.4 步长与滞后分组参数经验值
步长和滞后分组这两个参数经常被忽略,但它们直接决定半变异函数云的形态。ArcGIS 自动计算的步长通常偏大,因为它是用整个数据集的最大距离除以固定步数来定的。最大距离里可能包含大量相隔很远的点对,而这些点对在空间自相关上毫无意义,还会把曲线后半段压低。
我自己的习惯是先用近邻分析工具计算采样点平均最近邻距离,把步长设为这个值的 1/2 到 1 倍。比如之前处理过的一组地下水位数据,130 个点,自动步长算出来是 7800 米,但样点平均间距才 1200 米左右,我把步长改到 800 米,步数设为 15,半变异函数云立刻好看了很多,交叉验证的 RMSE 也明显下降。此外还要保证半变异函数曲线的最大滞后距离覆盖到变程。如果变程值远大于图上能看到的距离,就说明参数设置有问题,需要增加步数或加大步长。
5. 常见问题与排查技巧实录
5.1 插值图出现“靶心效应”
表现:预测图上以采样点为中心出现一圈圈的圆形波纹,像靶子一样。出现这种问题通常有四个原因:搜索邻域的最小邻居数设置得过小;半变异函数模型的块金效应占比过高;样点分布过于聚集,某些扇区内实际参与计算的点数不足;趋势阶次设置过高导致趋势项把局部信息全部吸走,残差成了纯噪声。
排查顺序建议:先看交叉验证里的每个样本点,那些误差特别大的点是否正好是孤立点;然后检查搜索邻域面板,把最小邻居数调大,比如从 2 调到 5,同时增加扇区数;如果还是没有改善,回头检查半变异函数,看块金效应占基台值的比例是不是已经超过了 50%,是的话考虑更换模型或检查原始数据是否有异常值。最后,如果趋势阶次设成了二阶,可以试试降为一阶或改为普通克里金对比,泛克里金过拟合的现象在图上通常会表现为大范围的高频波动,容易被误判成靶心效应。
5.2 交叉验证指标不理想
交叉验证指标不好的原因很杂,但可以按逻辑拆解。标准化均方根误差大于 1,说明模型的预测不确定性被低估了,实际误差比理论误差大,通常意味着半变异函数的基台值偏小或者块金值偏小,模型太“自信”;标准化均方根误差小于 1,说明模型过度平滑,预测面被压缩得过于平坦。
如果是前者,先检查数据分布是否满足正态性。泛克里金在计算趋势项和残差时对正态分布有一定要求,数据严重偏态会导致交叉验证指标失真。ArcGIS 的数据变换(Transformation)选项里可以选对数变换或 Box-Cox 变换,变换后要重新跑半变异函数拟合,不能只看变换前的结果。如果是后者,考虑降低步长、增加步数,让半变异函数曲线更贴近真实的空间变异结构,同时检查搜索邻域的最大邻居数是否设置得过大而抹平了细节。
另外,如果泛克里金的预测误差均值偏离 0 比较多,先怀疑趋势阶次选错了:一阶趋势倾向于低估弯曲趋势的峰值,二阶趋势容易出现边缘区域外推过度。可以在趋势分析工具里先看清趋势面形态再回来调整。
5.3 报错与局部 NoData 的处理
ArcGIS 在做插值时遇到明显的坐标系统缺失或错误坐标系会直接报错,通常提示无法计算距离。这种情况先把数据的坐标系检查一遍,确保定义正确而且已经投影到适合的投影坐标系。
局部 NoData 很常见,通常发生在数据覆盖范围的边缘,或者是像元大小设置得比样点间距还小导致部分像元周围没有足够样本。排查方法是打开预测标准差栅格:NoData 区域的边缘往往对应标准差很高的区域。可以先检查搜索邻域的半径设置,默认的邻域最大半径在 ArcGIS 里叫 Max neighbors,如果邻域半径不够大,边缘就会有大片 NoData。另一个解决办法是导出栅格时设置合适的像元大小,前提是先了解采样点密度,不要一上来就用 1 米像元去插一个间距 100 米的样点数据集。
6. 泛克里金与其他克里金方法的选型对比
6.1 泛克里金 vs 普通克里金 vs EBK
聊到选型,我把泛克里金(UK)、普通克里金(OK)和经验贝叶斯克里金(EBK)放在一张表里对比,这样看更直观:
| 对比项 | 普通克里金(OK) | 泛克里金(UK) | 经验贝叶斯克里金(EBK) |
|---|---|---|---|
| 均值假设 | 未知常量 | 坐标多项式的趋势函数 | 局部均值,通过子集模拟 |
| 趋势处理 | 不支持 | 显式建模趋势项 | 自动处理非平稳性 |
| 参数估计 | 单一半变异函数拟合 | 趋势项 + 残差的半变异函数 | 基于模拟的子集参数分布 |
| 适合数据 | 无明显趋势、平稳性较好的数据 | 有明显区域性趋势的数据 | 非平稳性强、数据量大、参数难调的数据 |
| 操作复杂度 | 低 | 中等 | 低,但计算开销大 |
| 边缘外推 | 回归到均值 | 沿趋势方向延续 | 较稳健 |
EBK 是 ArcGIS 近些年主推的方法,它对新手更友好,能自动处理参数不确定性。但它对数据量有一定要求,样本太少时模拟结果不稳定,而且计算量明显大于前两者。泛克里金的优势在于对趋势的显式建模,做科研和写报告时解释性强,可以明确说出“该区域存在一个由东向西逐步递减的趋势”,这是 EBK 的黑盒思路无法直接给出的结论。
6.2 什么场景别用泛克里金
泛克里金不是万能的,下面这些场景我建议直接绕开:
样本量太少,少于 20 个点时不要用。趋势项需要额外估计参数,样本量太小不仅趋势参数估不准,残差的变异函数也会很不可靠。这种情况下宁可做样条插值或反距离权重。
数据存在断裂面或强烈的局部突变,比如河流两侧的土壤属性差异极大,或者断层两侧的地下水埋深完全不同。这类非连续现象本质上不能用平滑的趋势函数描述,泛克里金会把断裂面抹成一个斜坡,反而不如直接分区插值。
趋势形态不是多项式形态,比如趋势沿某个方向呈周期性波动,或者受人类活动影响出现剧烈的不规则变化。此时高阶趋势项会开始拟合噪声,低阶趋势项又无法描述波动,回归克里金或结合辅助变量的协同克里金往往是更好的选择。
6.3 泛克里金的扩展玩法:回归克里金与辅助协变量
如果你的数据趋势很强,但形态又不是简单的一阶或二阶多项式,可以考虑手动实现“回归克里金”:
- 用普通线性回归或 GWR 建立目标变量与辅助变量(如高程、距河流距离)的回归模型;
- 从原始值中减去回归预测值,得到残差;
- 对残差做普通克里金插值;
- 将回归预测面与残差插值面叠加,得到最终结果。
这个方法在实际项目里非常实用。我在做一个降雨插值项目时,先用高程做回归,解释掉约 60% 的方差,残差的空间相关性明显减弱,再做普通克里金,交叉验证的 RMSE 比直接用泛克里金降低了将近 18%。如果你在 ArcGIS 里用地理加权回归工具,再把回归残差导出去做克里金,整个过程可以完全在 GIS 内完成,并不复杂。
另外,泛克里金也可以通过向导左侧方法列表切换到 Co-Kriging 类型,引入辅助变量做协同克里金。如果辅助变量与目标变量相关性高,比如电导率和土壤盐分,协同克里金能显著改善预测精度。注意辅助变量的采样密度最好高于目标变量,否则增强效果有限。
我个人在实际操作中的体会是:泛克里金的门槛不在操作,而在对数据的前期判断。拿到数据先花十分钟看趋势分析图,再决定要不要用泛克里金、用几阶趋势,这一小步能省掉后面大把调参时间。参数层面,半变异函数模型千万别只跑一个就把结果定下来,每次我都至少用三个模型交叉验证对比一遍再选最优。最后提醒一句,成果图一定把预测标准差一起输出,这是我见过最容易被忽略、但审图时最常被问到的一张图。