简介：这份资源面向从事Web三维地理信息开发的前端与GIS工程师，聚焦在Cesium场景中实现克里金插值并渲染连续温度场这一典型需求。它把离散温度测点经变异函数分析估算为平滑曲面，解决不规则采样数据难以直观呈现空间变化的问题，适合已具备Cesium基础、希望进阶空间插值与可视化能力的中级开发者参考。压缩包共585个文件、约7.27MB，以321个js脚本承载插值与渲染逻辑，配合141个png与53个jpg图像、29个json配置、27个css样式及svg、wasm等资源，构成可直接运行的完整示例工程。已有167人学习下载。读者可从中获得克里金算法与Cesium集成的实现思路、温度场平滑渲染的完整代码组织方式，以及样式与配置文件的目录结构参考，便于迁移到气象、环境监测等三维可视化项目中。
1. 克里金插值：为什么温度场渲染不能只靠“就近取色”
在地理可视化里，把几十个气象站的温度点铺到三维地球上，最直觉的做法是给每个点一个颜色，再让 GPU 在点与点之间做线性插值。结果往往是一片“色块拼接”的观感：站点密集处还行，站点稀疏处颜色过渡生硬，山脊两侧明明温差很大，屏幕上却是一条平滑的渐变带。这不是渲染管线的问题，而是插值模型的问题——线性插值只认距离，不认空间相关性。
克里金插值（Kriging）解决的正是这件事。它不把每个观测点当成孤立的采样，而是先估计整个区域的“空间变异结构”，再据此给每个未知点分配权重。换句话说，它回答的是“离我近的点该信多少、离我远的点该信多少、方向不同的点该不该区别对待”。把这套权重算出来的温度场贴到 Cesium 的球面上，色带才会跟着地形和气候格局走，而不是跟着站点分布的疏密走。搞懂这套机制，你才知道为什么同样的站点数据，换一种插值模型，渲染结果会差出一个量级。
2. 空间变异函数：克里金凭什么“知道”哪里该平滑、哪里该突变
2.1 从半变异函数说起：距离越远，差异该有多大
克里金的核心不是某个公式，而是一个先验假设：空间上相近的观测值，比相距较远的观测值更相似。这个假设被量化为半变异函数 γ(h)，它描述的是“相隔距离 h 的两点，其观测值之差的平方的期望的一半”。实际计算时用实验半变异函数估计：
# 实验半变异函数：对每一对距离落在 [h-dh, h+dh] 内的点，统计半方差 def experimental_variogram(coords, values, bins): gamma = [] for h in bins: pairs = [] for i in range(len(coords)): for j in range(i + 1, len(coords)): d = distance(coords[i], coords[j]) if abs(d - h) <= dh: # 距离落入当前滞后距区间 pairs.append((values[i] - values[j]) ** 2 / 2) gamma.append(sum(pairs) / len(pairs) if pairs else 0) return gamma这段代码的关键在(values[i] - values[j]) ** 2 / 2：取平方是为了让正负偏差都贡献变异，除以 2 是因为半变异函数定义上把成对差值的方差摊了一半。bins是人为划定的滞后距区间，区间宽度dh决定了曲线的平滑程度——太窄噪声大，太宽会抹掉短程结构。跑出来的散点还要拟合一个理论模型（球状、指数、高斯），拟合出的三个参数才是后续权重计算的输入。
2.2 块金、基台、变程：三个参数决定插值的“性格”
拟合出的理论变异函数有三个关键量。块金值（nugget）是 h 趋近 0 时的截距，代表测量误差和小于采样间距的微观变异；基台值（sill）是曲线趋于平稳的高度，代表总变异；变程（range）是曲线达到基台的距离，超过这个距离，两点就不再被认为空间相关。这三个参数直接决定插值场的“性格”：变程大，说明影响范围广，插值结果偏平滑；块金值大，说明短程噪声重，插值会更贴近局部观测点，甚至出现“牛眼”效应。
我一般会这么理解：变程是这套模型的“视野半径”，块金是它的“不信任阈值”。温度数据通常变程较大、块金较小，因为气温在几十公里尺度上是连续变化的；如果拟合出来块金异常大，多半是站点数据里混进了未质控的异常值，或者把不同海拔的站点硬塞进同一个平面模型。
2.3 普通克里金与泛克里金的取舍
普通克里金假设区域均值未知但恒定，权重完全由变异函数和点间几何关系决定；泛克里金则允许均值随坐标漂移，适合存在明显趋势面的场景。温度场恰恰常有趋势——纬度升高、海拔升高都会带来系统性降温。如果只用普通克里金，趋势会被当成随机变异吸收进残差，插值场在趋势方向上会偏。常见的实现思路是先拟合一个趋势面（比如温度的线性海拔梯度），对残差做普通克里金，再把趋势加回去。这一步在代码里体现为“去趋势—插值—回加”三段式，而不是直接对原始值调克里金。
3. 权重方程组：每个未知点的温度是怎么被“投票”出来的
3.1 普通克里金方程组的结构
给定 n 个观测点，要估计未知点 x0 的值，克里金把它写成观测值的加权和：Ẑ(x0) = Σ λi · Z(xi)。权重 λ 不是随便给的，而是通过求解一个带约束的线性方程组得到，约束是权重之和为 1（保证无偏）。方程组长这样：
import numpy as np def ordinary_kriging_weights(coords, x0, variogram_model): n = len(coords) A = np.ones((n + 1, n + 1)) b = np.ones(n + 1) for i in range(n): for j in range(n): A[i][j] = variogram_model(distance(coords[i], coords[j])) # 观测点之间的半变异 A[i][n] = 1.0 # 拉格朗日乘子列 A[n][i] = 1.0 b[i] = variogram_model(distance(coords[i], x0)) # 观测点到未知点的半变异 A[n][n] = 0.0 b[n] = 1.0 weights = np.linalg.solve(A, b) return weights[:n] # 前 n 个是各观测点权重，最后一个是拉格朗日乘子矩阵 A 的左上 n×n 块是观测点两两之间的半变异值，最后一行一列是约束条件引入的拉格朗日乘子。右端向量 b 的前 n 项是观测点到目标点的半变异值。解这个方程组，得到的就是“每个观测点该占多少话语权”。注意 A 的对角线是 0（点与自身的半变异为 0），这让矩阵在数值上容易接近奇异，实际实现里常加一个极小的正则项或改用 Cholesky 分解。
3.2 为什么权重可能为负，以及负权重意味着什么
普通克里金的权重并不保证非负。当某个观测点与目标点的半变异值小于它与其他点的平均半变异时，解出来的 λ 可能为负。负权重在数学上是为了满足无偏约束，但在物理上意味着“这个点的值要往反方向拉”。温度插值里出现负权重，往往说明变异函数模型拟合得偏了，或者观测点分布极度不均。工程上常见的处理是对权重做截断再归一化，但这会破坏克里金的最优性，属于用精度换稳定性的取舍。
3.3 从权重到温度场：逐像素求解的代价
三维地球上要渲染的是一整片区域，不是单个点。如果对每个像素都解一次 n+1 阶方程组，计算量会随分辨率平方增长。实际实现通常分两步：先在规则网格上做克里金插值（网格分辨率远低于屏幕像素），把结果存成一张二维温度栅格；再把栅格作为纹理贴到 Cesium 的几何体上，由 GPU 做双线性采样。这样克里金的计算量只跟网格大小和站点数有关，跟相机分辨率解耦。网格步长的选择是个经验活：步长太大，色带会糊；步长太小，站点稀疏区会出现插值伪影，而且纹理上传开销上升。
4. 落到 Cesium：温度栅格怎么变成球面上的色带
4.1 从经纬度网格到纹理：数据布局的转换
克里金算出来的是一张以经纬度为索引的二维数组。Cesium 的材质系统吃的是纹理，所以要把数组映射成 RGBA 像素。温度值本身不直接当颜色通道用，而是先归一化到 [0,1]，再通过一张色带查找表（colormap）转成颜色。常见的实现思路是把色带做成 256×1 的一维纹理，在片元着色器里用归一化温度做纹理坐标采样，这样换色带不用重算克里金，只换一张查找表即可。
// 片元着色器片段：用归一化温度采样色带纹理 uniform sampler2D u_colormap; varying vec2 v_st; // 网格纹理坐标 varying float v_temp; // 插值后的归一化温度，来自顶点属性或数据纹理 void main() { float t = clamp(v_temp, 0.0, 1.0); vec4 color = texture2D(u_colormap, vec2(t, 0.5)); // 一维色带按 t 采样 gl_FragColor = color; }这里v_temp的来源有两种常见做法：一是把温度栅格直接作为数据纹理传入，在片元里采样得到温度再查色带；二是把温度编码进顶点属性，由光栅化插值。前者适合栅格分辨率固定的场景，后者适合网格顶点数可控的规则网格。两种做法的共同点是：克里金的重计算发生在 CPU 侧或预计算阶段，GPU 只负责采样和着色，避免每帧重算权重。
4.2 地理配准：矩形纹理怎么贴到曲面上
经纬度网格是等间距的矩形，但地球是曲面，直接贴上去会在高纬度产生拉伸。Cesium 的Rectangle几何体配合GeographicTilingScheme可以按经纬度范围贴图，内部会处理投影变形。关键参数是矩形的西、南、东、北四个边界，必须和克里金网格的经纬度范围严格对齐，否则色带会和地形错位。我踩过的坑是网格按像素中心对齐、矩形按像素边界对齐，差半个像素，放大后能看到色带整体偏移一条缝。
4.3 高度与遮挡：温度场该贴在地表还是悬空
温度是近地面量，通常贴在地形表面。做法是把矩形几何体的height设为 0 并开启heightReference: CLAMP_TO_GROUND，让每个顶点贴合地形。但克里金网格本身是平面的，贴到起伏地形上会出现“色带穿过山体”的观感——山脊处的色带被地形切断。更细的做法是把温度栅格作为地形材质的叠加层，或者用自定义 Primitive 在片元里按地形高度做混合。这属于渲染层的取舍：贴地简单但视觉上不连续，叠加层自然但实现复杂度高。
5. 避坑与排查：克里金温度渲染里最容易翻车的几处
5.1 色带整体偏冷或偏热：归一化区间选错了
现象是整片区域颜色集中在色带一端，看不出层次。原因多半是归一化用了全局固定区间，而当前区域的温度跨度远小于全局跨度。解决方法是按当前数据集的实际最小最大值做动态归一化，或者用分位数截断（比如取 2% 和 98% 分位）避免极端值把区间撑爆。注意动态归一化会让不同时刻的同一温度对应不同颜色，做时间序列动画时要固定区间，否则色带会“呼吸”。
5.2 站点稀疏区出现同心圆：块金值被低估
现象是插值场在孤立站点周围出现一圈圈“牛眼”。原因是变异函数拟合时块金值给得太小，模型过度信任单个观测点，导致该点权重畸高。解决办法是重新检查实验半变异函数的短程段，适当抬高块金值，或者对孤立站点做邻域限制——只让一定距离内的站点参与权重计算，避免远距离点被赋予异常权重。
5.3 权重方程组求解报奇异矩阵：重复坐标或共线点
现象是np.linalg.solve抛LinAlgError。原因是存在坐标完全相同的站点（数据合并时没去重），或者多个站点近似共线导致矩阵秩亏。解决方法是插值前做坐标去重，对距离小于阈值（比如 1 米）的点做合并取均值；矩阵求解改用带正则的lstsq或加对角扰动。这个坑在站点数据来自多个来源拼接时特别常见。
5.4 色带与地形错位：矩形边界和网格范围没对齐
现象是放大后色带边缘和海岸线、山脊对不上。原因是克里金网格的经纬度范围按格点中心计算，而 Cesium 矩形按格点边界计算，两者差半个格距。解决办法是统一约定：网格范围取格点边界，矩形直接用同一组边界值；或者在生成纹理时把采样点偏移半个像素。这个偏差在小比例尺下看不出来，放大到城市级就会暴露。
5.5 帧率骤降：每帧重算克里金
现象是相机移动时卡顿。原因是在渲染循环里调用了克里金插值。克里金是 O(n³) 级别的计算（解方程组），绝不能放进每帧逻辑。正确做法是把插值结果缓存成纹理，只在数据更新或时间步变化时重算；相机移动只触发纹理采样，不触发权重计算。如果确实需要动态更新，把克里金放到 Web Worker 里异步算，主线程只负责上传纹理。
6. 设计取舍：克里金不是越精确越好
把克里金接进 Cesium 的过程中，最值得反复权衡的不是公式推导，而是“精度花在哪里”。变异函数拟合得再精细，如果站点本身只有几十个，插值场的可信度上限就摆在那里；网格步长缩到屏幕像素级，计算和纹理开销翻几倍，视觉收益却很快饱和。我一般的做法是先把变异函数和网格步长做成可配置参数，用一组已知真值的模拟温度场做交叉验证，看不同参数下插值误差和渲染帧率的曲线拐点，再定默认值。
另一个取舍在实时性上。克里金天然是批处理模型，适合静态或低频更新的温度场；如果要做逐秒刷新的动态场，要么降低网格分辨率，要么退回到更轻量的反距离加权做过渡，把克里金留给关键帧。这不是妥协，而是承认不同插值模型各有适用边界——克里金的价值在于它把“空间相关性”显式建模了，代价是每次求解都要解一个稠密方程组。
从那以后我每次接插值类渲染需求，都强制先画一张实验半变异函数散点图再动手写渲染代码，因为那张图会直接告诉你数据到底支不支持你想要的平滑度。希望帮到你。
本文还有配套的精品资源，点击获取
