Cesium 克里金插值完整指南：从原理到气象热力图渲染

简介这是一套基于Cesium三维地球库与克里金插值算法实现温度数据三维渲染的完整资源包面向网页地理信息系统开发者、三维可视化工程师以及气象与地学数据分析人员。克里金插值通过对离散温度点的变异函数分析将不规则站点数据转化为平滑连续的温度场直观呈现空间温度分布与变化趋势解决离散气象站点数据难以直接生成连续三维温度场的问题。压缩包共585个文件大小约7.27MB包含321个JS脚本核心插值逻辑与Cesium功能模块、27个CSS样式表、29个JSON配置以及大量贴图资源、SVG图标、WASM与HTML等辅助资源目录结构清晰与标准Cesium工程一致。已有167人浏览学习。通过使用这套资源可快速掌握克里金插值在三维GIS中的实践方法并直接复用其数据解析、插值渲染与交互设计代码。1. 从站点温度到连续热力面Cesium Kriging 插值渲染的整体思路如果你手头只有几十个气象站的温度观测点却想在三维地球上渲染出一整片连续的温度云图用普通的反距离加权或最近邻插值出来的色块会像打补丁一样生硬。克里金插值Kriging的价值恰恰在于它不只根据距离远近给权重还引入空间自相关性——距离近的点权重更高但点与点之间的空间结构变差函数会参与估算因此插值结果更平滑、更符合气象场的连续变化规律。我拆过的这个 Cesium Kriging 项目就是把克里金算法从一个黑匣子变成一张可以贴到地球表面的纹理图整个过程涉及数据格式整理、半方差计算、权重矩阵求解、canvas 热力色带映射以及 Cesium 图元更新。适合想用纯前端方案做气象可视化、又不愿意把插值丢给后端服务的开发者。2. 克里金插值原理与选型为什么选普通克里金数据格式怎么备2.1 克里金家族怎么选普通克里金是绝大多数可视化场景的默认起点克里金并非单一算法而是一族基于变异函数Variogram的地统计插值方法。简单克里金假设区域化变量的均值已知且恒定普通克里金假设均值未知但恒定这更贴近气象温度场——你不可能先验地知道整个区域的平均温度。泛克里金则允许均值随坐标漂移趋势面适合有明显高程分带或纬度分带的场景但代价是模型复杂度上升前端求解不稳定。实际落地上普通克里金Ordinary Kriging是可视化渲染的默认起点求解相对稳定、参数少、插值结果视觉上平滑且对站点数量不敏感。选型时还有一个被忽视的考量克里金权重矩阵求解本质上是在解一个 n 元线性方程组站点数量超过 300 个后矩阵求逆的计算量暴涨前端 JS 会明显卡顿。我之前做过一次 500 个站点的温度场实验纯 JavaScript 解方程花了近 2 秒直接卡掉用户交互。所以如果站点数超过 300要么分块插值要么改用插值半径裁剪。这个我在后面排查章节里再展开。2.2 数据格式准备GeoJSON 站点数据和面状边界一步到位Cesium Kriging 渲染的输入是两个部分站点观测数据和渲染边界。站点数据规范做法是用 GeoJSON 的 FeatureCollection每个 Feature 的 geometry 是 Pointproperties 里带温度值。边界数据则用于裁剪插值结果——克里金插值本身是铺满整个包围盒的如果不裁剪海洋区域也会被强行插值产生毫无意义的色块。// GeoJSON 站点数据示例温度单位摄氏度 { type: FeatureCollection, features: [ { type: Feature, properties: { temperature: 23.5 }, geometry: { type: Point, coordinates: [116.4, 39.9] } } ] }这里有个关键点properties 里的字段名必须和插值脚本里读取的字段名一致。否则你在插值代码里传feature.properties.temperature而数据里字段叫temp脚本一眼看过去不会报错结果值全是NaN。这类问题很难查我一般拿到数据先做一次字段名合法性检查——直接遍历 features 打印所有 key。边界数据建议用 Polygon GeoJSONWeb Mercator 投影或 WGS84 经纬度坐标皆可Cesium 底层能处理重投影。但有一个坑需要注意边界 Polygon 的坐标环必须闭合即首尾坐标点相同。很多从 GIS 软件导出的边界数据首尾没闭合克里金裁剪阶段绘制的 canvas 路径就会异常最终出现半个面被填充的怪相。2.3 变差函数Variogram模型球状、指数、高斯还是线性克里金的核心是半方差函数用来刻画两个点之间空间差异随距离变化的规律。半方差计算出来后需要用理论模型去拟合。常用的有四种模型公式特征适用场景前端计算成本球状Spherical有基台值有变程超过变程后稳定气象温度、降水插值最常用低指数Exponential渐近逼近基台值变程为无穷空间连续性强、过渡平缓的场低高斯Gaussian抛物线型拟合短距离内变化极缓高程、污染浓度等平滑物理场中线性Linear无基台值最简快速测试、站点少最低在 Cesium Kriging 这个场景下我一般直接采用球状模型原因只有一个简单气象要素场在空间上具有有限变程超过一定距离后温度就彼此独立了这个物理意义对应球状模型最自然。而且球状模型的求解最稳定给前端用的库也默认支持它。如果你发现插值结果出现了明显的“牛眼”一个个同心圆色斑大概率是变程range设太小找到合适的变程值需要多试几次。3. 克里金算法核心实现从半方差到权重矩阵3.1 半方差计算与经验变差函数算法落地的第一步是把所有站点两两配对计算每对站点之间的距离和半方差。半方差的定义是两点温度差平方的一半。这一步类似于求解一个距离矩阵复杂度是 O(n²)。站点数为 100 时计算量是 4950 对JavaScript 毫秒级就能完成500 个站点就到了 12 万对开始有压力。function computeSemivariance(points) { const n points.length; const pairs []; for (let i 0; i n; i) { for (let j i 1; j n; j) { const dx points[i].lon - points[j].lon; const dy points[i].lat - points[j].lat; const dist Math.sqrt(dx * dx dy * dy); const semivariance 0.5 * Math.pow(points[i].value - points[j].value, 2); pairs.push({ dist, semivariance }); } } // 按距离分组计算每个距离桶的平均半方差 const buckets {}; pairs.forEach(p { const key Math.round(p.dist / 0.5) * 0.5; // 0.5 度距离桶 if (buckets[key]) { buckets[key].sum p.semivariance; buckets[key].count; } else { buckets[key] { sum: p.semivariance, count: 1 }; } }); const variogram Object.keys(buckets).map(key { return { dist: parseFloat(key), semivariance: buckets[key].sum / buckets[key].count }; }).sort((a, b) a.dist - b.dist); return variogram; }这段代码做了三件事计算距离和半方差、按距离分桶、对每个桶求平均。分桶的目的是把散点拟合的经验变差函数离散化避免每个距离都有唯一值导致后续拟合不稳定。桶宽0.5 度是经验值——经纬度坐标下 0.5 度约 55 公里对于省级范围的温度场够用。如果桶宽太大变程信息会被抹平太小则每个桶内样本数不足拟合曲线抖动厉害。3.2 球状模型的拟合最小二乘法求基台和变程经验变差函数计算出来后需要拟合球状模型参数。球状模型公式γ(h) C0 C · S(h)其中 h 是距离C0 是块金值nugget代表测量误差或微观尺度变异C 是偏基台值partial sillC0C 是基台值sill变程range是模型到达基台值的距离。拟合原理是最小二乘法初始化一组参数计算模型预测值求出与经验值的残差平方和然后迭代搜索参数空间让残差平方和最小。这里不用自己实现完整的非线性优化简化做法是变程固定为经验变差函数曲线开始持平的距离附近基台值固定为经验半方差的均值上限然后用迭代逼近块金值。function fitSphericalModel(variogram) { // 经验基台值取所有点半方差的最大值 const maxSemi Math.max(...variogram.map(v v.semivariance)); // 变程找到第一个接近基台值的距离点 const range variogram.find(v v.semivariance 0.95 * maxSemi)?.dist || 10; // 初始参数块金值取最小值偏基台值 基台 - 块金 let nugget variogram[0] ? variogram[0].semivariance : 0; const sill maxSemi; let partialSill sill - nugget; // 网格搜索微调在前 80% 范围内逐步优化局部参数 let bestResidual Infinity; let bestRange range; for (let r 0.5; r range * 2; r 0.5) { const residual variogram.reduce((sum, v) { const modelValue nugget partialSill * sphericalModel(v.dist, r, sill); return sum Math.pow(modelValue - v.semivariance, 2); }, 0); if (residual bestResidual) { bestResidual residual; bestRange r; } } return { nugget, partialSill, range: bestRange }; } function sphericalModel(h, range, sill) { if (h range) { return sill * (1.5 * h / range - 0.5 * Math.pow(h / range, 3)); } return sill; }参数说明range决定了插值的影响半径这个值对最终渲染效果影响极大。range 太小插值面破碎出现大量局部斑点range 太大插值面过于平滑温度梯度细节全丢。网格搜索步长 0.5 度是一个平衡点追求精细可以调到 0.1但迭代次数会翻好几倍。实际项目中我发现 0.5 步长配合二次精调已经足够——先用大步长找大致范围再在小范围内用更细步长精调。3.3 权重矩阵构建与线性方程组求解克里金插值的核心方程组是一个特殊的线性系统对角线之外的元素是两点间的半方差值由拟合好的模型算出维度是n1×n1其中 n 是参与插值的站点数。最后一列和最后一行是拉格朗日乘子项用来保证权重之和为 1。解方程组的标准做法是高斯消元。这里有一个关键性能问题每次插值一个目标点都要重新构建矩阵并求解一次。假设渲染一张 500×500 的网格图就是 25 万个目标点每个目标点都要解一个n1元方程组——这在纯前端是不可行的。所以实际工程中不会对角点级别去求权重而是用现有的开源库实现——这个 Cesium Kriging 项目里用的是 kriging.js 的实现方式。它把求解过程拆成两步先生成所有站点两两之间的克里金矩阵并做 Cholesky 分解这一步和具体目标点无关然后对每个目标点只需要代入矩阵分解结果做一次前向/回代。这样矩阵分解只做一次后续每个目标点的求解复杂度从 O(n³) 降到了 O(n²)。// 使用 kriging.js 风格的实现流程 function createKrigingModel(lon, lat, values, variogramModel) { const n lon.length; const distanceMatrix []; // 构建距离矩阵 for (let i 0; i n; i) { distanceMatrix[i] []; for (let j 0; j n; j) { distanceMatrix[i][j] distance(lon[i], lat[i], lon[j], lat[j]); } } // 构建克里金矩阵用模型计算每个距离对应的半方差 const krigingMatrix []; for (let i 0; i n; i) { krigingMatrix[i] []; for (let j 0; j n; j) { krigingMatrix[i][j] variogramModel(distanceMatrix[i][j]); } } // 添加拉格朗日约束行/列 for (let i 0; i n; i) { krigingMatrix[i][n] 1; krigingMatrix[n] krigingMatrix[i]; } krigingMatrix[n][n] 0; // Cholesky 分解矩阵注意矩阵必须是正定的 const lower choleskyDecomposition(krigingMatrix); return { lower, lon, lat, values }; } // 后续对每个目标点只需做一次快速求解 function predict(lon, lat, model) { const n model.lon.length; const b []; for (let i 0; i n; i) { b[i] variogramModel(distance(lon, lat, model.lon[i], model.lat[i])); } b[n] 1; // 拉格朗日约束 const weights solveWithCholesky(model.lower, b); return weights.reduce((sum, w, i) sum w * model.values[i], 0); }这里的逻辑重点在于 Cholesky 分解要求矩阵正定。如果数据中有两个站点坐标完全相同——这种情况在合并多源数据时非常常见——矩阵就会奇异Cholesky 分解直接失败或产生 NaN。遇到这种数据不要硬解先做坐标去重再跑模型。4. Cesium 渲染实现canvas 热力纹理与多层叠加4.1 插值结果转 canvas色带映射与透明度处理克里金计算完成后输出通常是一个二维数组网格点的预测值。要将它渲染到 Cesium 上中间有个关键的转换步骤把网格数组映射成 canvas 像素图。这一步需要处理数值到 RGBA 颜色的映射以及无效区域的透明化。function renderToCanvas(gridData, colorScale, width, height) { const canvas document.createElement(canvas); canvas.width width; canvas.height height; const ctx canvas.getContext(2d); const image ctx.createImageData(width, height); for (let y 0; y height; y) { for (let x 0; x width; x) { const index (y * width x) * 4; const value gridData[y][x]; // 无效值 / 超出边界区域的像素设为全透明 if (value null || isNaN(value)) { image.data[index] 0; image.data[index 1] 0; image.data[index 2] 0; image.data[index 3] 0; continue; } // 将数值映射到色带上的位置如蓝-绿-黄-红 const color interpolateColorScale(value, colorScale); image.data[index] color.r; image.data[index 1] color.g; image.data[index 2] color.b; image.data[index 3] 200; // 半透明便于叠加底图 } } ctx.putImageData(image, 0, 0); return canvas; }这里最容易翻车的是色带映射。温度值的范围是连续的但不同地区不同季节的温差很大。如果把色带范围写死为 -20℃ 到 40℃夏季南方城市的插值结果整个区域都落在色带的一个小范围内看起来全是一片深红色温度梯度完全看不出来。正确处理方式是需要根据实际数据动态计算色带上下限我一般用数据的最小值和最大值再加 10% 的边距防止极值点溢出。4.2 Cesium 实体加载图元坐标绑定 canvas 纹理canvas 生成后需要把它作为 Cesium 的表面纹理绑定到地球或者一个矩形区域上。常用做法有两种用 RectangleGeometry 生成矩形图元配合 Material 动态纹理或者用 Entity 的 polygon 配合 MaterialProperty。我推荐用 entity 方式因为 entity 支持数据源监听后续数据刷新时可以直接改 property不需要删了重建。// 用 entity custom material 把 canvas 纹理贴到地球表面 function addHeatmapLayer(viewer, bounds, canvas) { const material new Cesium.Material({ fabric: { type: Image, uniforms: { image: canvas } } }); const entity viewer.entities.add({ rectangle: { coordinates: Cesium.Rectangle.fromDegrees( bounds.west, bounds.south, bounds.east, bounds.north ), material: material, // 避免与地形遮挡冲突 classificationType: Cesium.ClassificationType.BOTH } }); return entity; }重点参数是bounds——它定义了纹理在地球上的贴图范围必须与插值计算时使用的包围盒完全一致否则纹理和站点位置对不上。另一个坑是 Cesium 的 Material 对 canvas 尺寸有要求canvas 宽高必须是 2 的整数次幂如 256×256、512×512、1024×1024。不满足这个条件时纹理在部分设备上会出现边缘模糊或采样异常。插值的网格密度决定了 canvas 分辨率网格 64×64 对应 canvas 64×64不是 2 的幂时需要补齐网格 100×100 则建议直接冲刷成 128×128 的 canvas多出来的区域留白。4.3 多层叠加与透明度站点层、插值层、边界层如何协同实际渲染时往往需要不止一层。我对这个项目的处理方案是三层叠加底层读 Cesium 自带影像地形 影像底图提供地理上下文。中间层克里金插值热力面透明度控制在 180220 之间既能透出底图地形又能清晰显示温度色块。顶层站点标注层Biz 图标 温度数字用 BillboardCollection 或 Entity 的 label 实现。站点层和插值层存在视觉遮挡问题——label 标签如果设置disableDepthTestDistance: Number.POSITIVE_INFINITY标签永远浮在最上方不受地形遮挡但这也会导致热力面被标签大面积盖住。我的习惯做法是label 只在鼠标悬停或点击时显示平时只显示小圆点。这样可以避免密密麻麻的标签把热力图糊死。三层叠加还有一种常见的视觉污染色带的绿色和底图的植被颜色撞在一起导致温度场的边界模糊。解决方式不是去改底图而是给色带加一个偏暗或偏亮的前景色去区分。我顺手做了一个测试对比表灰度底图 蓝-红温度色带的辨识度最高彩色底图 半透明温度色带其次彩色底图 不透明温度色带最差。5. 避坑与常见问题排查插值结果与渲染效果不符时先查这五个地方5.1 现象温度云图出现明显的圆圈波纹像水面涟漪原因变程range设置过小导致每个目标点只能被附近几个站点影响局部权重过度突出插值面不平滑。这本质上是变差函数拟合没有收敛到真实的变程值。解决打开调试日志输出经验变差函数散点图和拟合曲线。如果散点图在某个距离后明显水平而拟合模型的变程远小于这个距离手动把变程range改为散点图水平处的距离值。在代码里也可以做一个自适应变程取经验变差函数 95% 基台值对应的距离作为变程下限。5.2 现象Cesium 上热力图位置与站点标记错位偏移了约半个画面原因canvas 纹理的宽高比和Rectangle.fromDegrees定义的矩形宽高比不一致。Cesium 材质贴图默认拉伸填充整个矩形如果纹理是正方形、矩形区域是 2:1那么图像在水平方向被压缩导致视觉偏移。这不是克里金插值的错是纹理映射参数设置的问题。解决canvas 宽高比必须与矩形范围经纬度比例匹配。比如矩形范围东西跨度 10 度、南北跨度 5 度地形比例约为 2:1canvas 也应该用 2:1 的尺寸如 512×256。不要直接使用 1:1 canvas 去贴非 1:1 的区域。5.3 现象插值结果正确但某些区域显示为黑色撕裂块或整体呈暗色原因canvas 创建时没有使用全透明背景填充或者putImageData之前没有把像素初始化。canvas 默认背景是透明的 blackrgba(0,0,0,0)但如果你在绘制时遗漏了部分区域这些区域会透出黑色底。另外色带范围如果最小值不是从 0 开始低温区域映射出来的颜色可能偏深看起来像黑块。解决插值前先把 canvas 的每个像素初始化成全透明。校验色带映射结果在开发模式下允许控制台打印温度最小值、最大值以及对应颜色值看是否落在合理的色带范围。5.4 现象站点数量多时交互卡顿鼠标拖动地球明显掉帧原因插值重算被放在了鼠标移动事件里每次移动都触发一次 O(n²) 的半方差计算和 O(n³) 的矩阵分解。用户拖动 Cesium 相机时几十次重算直接把主线程塞满。解决把重算逻辑做防抖 异步。先在 mousemove 事件里设置一个标记位相机停止移动 500 毫秒后才触发重算。重算本身放在requestAnimationFrame或 Web Worker 里执行主线程只负责传递数据。另外站点数超过 300 时务必做分块插值把一个大区域切成 4 块每块分别插值后再拼接性能提升明显。5.5 现象边界外的区域被插值填满海洋上全是颜色原因克里金计算本身不知道你的边界在哪它只会按矩形网格输出边界外依然有预测值。前面提过要裁剪但很多人因为裁剪代码粘贴错误而跳过这一步。解决使用 canvas 的clip路径方法——先在 ctx 上描绘边界 Polygon 的路径然后调用ctx.clip()再进行drawImage或putImageData。这样裁剪后边界外保持全透明。注意Polygon 的坐标必须是经纬度转换到 canvas 平面坐标的转换公式要和纹理映射的 bounds 一致否则裁剪会错位。6. 进阶技巧交互采样验证插值可靠性用临时色带辅助体验最后一个值得分享的技巧是插值结果不能只看颜色要靠“交互采样”来验证。做法并不复杂在 Cesium 上绑定一个 click 事件点击点位的经纬度通过viewer.camera.pickEllipsoid获取然后把这个经纬度送到克里金模型里做单点预测同时把预测值与距离最近的实际站点值做差值显示在页面上。如果差值大于 3℃说明插值在这一片区域的可信度偏低你需要回头检查站点密度或变程设置。我在实际项目里通常给页面加一个“验证模式”开关开启时点击地图任意位置显示克里金预测值、最近站点实际值、两者偏差等偏差绝对值全部稳定在 2℃ 以内才认为这次插值配置是合格的。这个习惯帮我排查掉了好几次“看起来渲染很美观、实际数据失真”的问题。function addSamplingValidator(viewer, krigingModel, disabled) { const handler new Cesium.ScreenSpaceEventHandler(viewer.scene.canvas); handler.setInputAction((movement) { const cartesian viewer.camera.pickEllipsoid(movement.position, viewer.scene.globe.ellipsoid); if (!cartesian) return; const cartographic Cesium.Cartographic.fromCartesian(cartesian); const lon Cesium.Math.toDegrees(cartographic.longitude); const lat Cesium.Math.toDegrees(cartographic.latitude); const predicted krigingModel.predict(lon, lat); // 找到最近的站点并输出偏差 const nearest krigingModel.findNearestStation(lon, lat); console.log(预测值: ${predicted.toFixed(2)}℃, 最近站点: ${nearest.value.toFixed(2)}℃, 偏差: ${(predicted - nearest.value).toFixed(2)}℃); }, Cesium.ScreenSpaceEventType.LEFT_CLICK); }这段代码的价值不只是验证它顺带起到了“透明度调节参考”的作用。当鼠标在一个颜色边界附近采样你会发现即使经纬度相差很小插值结果也不是平滑过渡——如果这个“阶梯”出现在温度梯度大的区域比如山区反而说明插值真实反映了地形影响如果出现在平坦地带就要怀疑变程参数是否合理了。另外再提一个体验细节给热力层加一个“临时透明度”滑块不用每次改代码。做法是把Material的 alpha 属性暴露出来绑定一个input typerange的change事件。这样在展示给非技术人员看的时候你可以一键调整透明度找到“底图可见又不过度抢眼”的平衡点。从那次做某跨平台系统的气象模块之后我每次做插值可视化都把“采样验证”作为必走流程——不验证就发布迟早会被用户指出某个点位和实际气温对不上。希望你也能把这一步养成习惯少走弯路。本文还有配套的精品资源点击获取