简介这份资源围绕 Cesium 与克里金插值展开面向需要在三维地球场景中呈现连续温度场的 WebGIS 开发者与地理信息专业学习者。它解决的核心问题是如何把离散、不规则分布的温度观测点通过变异函数分析插值成平滑连续的温度表面并叠加到 Cesium 的交互式三维地球上适用于气象可视化、环境监测与空间分析类项目。压缩包共 585 个文件约 7.27MB以 321 个 js 脚本、141 个 png 与 53 个 jpg 图像、29 个 json 配置、27 个 css 样式为主另含少量 svg、wasm、html 等基本覆盖 Cesium 运行所需的脚本、贴图、样式与数据配置。目前已有 167 人学习下载。读者可从中获得克里金插值在 Cesium 中的完整实现思路、温度场渲染的代码组织方式与资源目录结构便于对照搭建自己的三维温度可视化应用并在此基础上扩展其他空间插值场景。1. 克里金插值为什么它比反距离权重更适合温度场渲染在 Cesium 里做温度数据渲染很多人第一反应是反距离权重插值——简单、快、代码短。但如果你真拿它渲染过气象站数据会发现一个很别扭的现象插值出来的温度场在站点附近还行一到站点稀疏的区域就开始出现明显的“牛眼”伪影等值线像被站点拽着走。这不是代码写错了而是反距离权重只考虑距离不考虑空间相关性。克里金插值解决的就是这个问题。它的核心结论很反直觉插值权重不是靠距离公式算出来的而是靠数据自身的空间结构反推出来的。换句话说它先让你从已知站点里“学”一个变异函数描述温度随距离怎么变化再用这个函数去决定每个未知点的权重。温度场天然有空间自相关性——近处比远处更像——克里金恰好把这件事数学化了。在 Cesium 这种三维地球场景里温度数据渲染的难点从来不是画点而是把离散站点变成连续场。克里金插值配合 Cesium 的 ImageryLayer 或 GroundPrimitive能把站点数据变成一张可叠加、可查询的温度栅格。搞懂它的底层机制你才能判断什么时候该用它、什么时候该退回反距离权重以及为什么你的插值结果会出现负温度或者条带状伪影。2. 变异函数克里金插值的“黑匣子”到底在算什么2.1 从半方差到空间自相关克里金插值的起点不是权重而是半方差函数。对任意两个站点计算它们温度差的平方的一半再按距离分组求平均就得到实验变异函数。这个函数描述的是距离为 h 的两点温度差异的期望有多大。常见的实现思路是先把所有站点两两配对按距离分箱每个箱里算平均半方差。距离近的箱半方差小说明近处温度相似距离远的箱半方差大说明远处温度差异大。当距离超过某个阈值后半方差趋于平稳这个阈值叫变程代表温度场的空间影响范围。// 计算实验变异函数输入站点数组每个点含 lon, lat, value function experimentalVariogram(points, binWidth 0.5) { const pairs []; for (let i 0; i points.length; i) { for (let j i 1; j points.length; j) { const dx points[i].lon - points[j].lon; const dy points[i].lat - points[j].lat; const h Math.sqrt(dx * dx dy * dy); // 距离 const diff points[i].value - points[j].value; pairs.push({ h, gamma: 0.5 * diff * diff }); // 半方差 } } // 按距离分箱求平均 const bins {}; pairs.forEach(p { const key Math.floor(p.h / binWidth) * binWidth; if (!bins[key]) bins[key] { sum: 0, count: 0 }; bins[key].sum p.gamma; bins[key].count 1; }); return Object.keys(bins).map(k ({ distance: parseFloat(k), gamma: bins[k].sum / bins[k].count })).sort((a, b) a.distance - b.distance); }这段代码里binWidth控制距离分箱的粒度太小会导致每个箱里样本太少、变异函数抖动太大则丢失细节。gamma就是半方差它衡量的是温度差异的平方期望。注意这里用的是平方差不是绝对差因为克里金的理论基础是二阶平稳假设平方差才有解析性质。2.2 理论变异函数模型为什么必须拟合实验变异函数是一堆散点直接拿它算权重会不稳定因为距离稍远一点可能就没有样本对了。所以必须拟合一个理论模型常见的有球状模型、指数模型、高斯模型。球状模型最常用因为它有明确的变程超过变程后半方差恒定。球状模型的公式是当 h a 时γ(h) C0 C * (1.5 * h/a - 0.5 * (h/a)^3)当 h a 时γ(h) C0 C。其中 C0 是块金值代表随机误差C 是拱高代表空间结构带来的方差a 是变程。拟合方法通常用最小二乘把实验变异函数的点代入模型优化 C0、C、a 三个参数。这一步是克里金插值里最“玄学”的地方——初值选不好拟合结果可能完全跑偏。我一般会先把块金值设为 0变程设为最大距离的一半拱高设为实验变异函数的最大值然后再迭代。// 球状模型拟合用简单的最小二乘网格搜索 function fitSpherical(variogramPoints) { const maxDist Math.max(...variogramPoints.map(p p.distance)); const maxGamma Math.max(...variogramPoints.map(p p.gamma)); let best { C0: 0, C: maxGamma, a: maxDist / 2, error: Infinity }; // 粗网格搜索 for (let C0 0; C0 maxGamma * 0.3; C0 maxGamma * 0.05) { for (let C maxGamma * 0.5; C maxGamma * 1.5; C maxGamma * 0.05) { for (let a maxDist * 0.2; a maxDist; a maxDist * 0.05) { let error 0; variogramPoints.forEach(p { const h p.distance; let model; if (h a) { model C0 C * (1.5 * h / a - 0.5 * Math.pow(h / a, 3)); } else { model C0 C; } error Math.pow(model - p.gamma, 2); }); if (error best.error) { best { C0, C, a, error }; } } } } return best; }网格搜索虽然慢但胜在稳定不会像梯度下降那样因为初值问题翻车。实际项目中如果站点数量上千可以先用粗网格定位再用局部细化。拟合出来的三个参数直接决定后续权重矩阵的形态所以这一步值得多花时间验证。2.3 普通克里金方程组权重到底怎么解有了变异函数模型就可以对每个未知点建立克里金方程组。假设有 n 个已知站点未知点 x0方程组是对于 i 1 到 nΣ λj * γ(xi, xj) μ γ(xi, x0) 以及约束条件Σ λj 1其中 λj 是权重μ 是拉格朗日乘子γ 是变异函数值。写成矩阵形式就是 [Γ 1; 1^T 0] * [λ; μ] [γ0; 1]。解这个线性方程组就能得到权重。// 求解普通克里金权重 function solveKrigingWeights(stations, target, model) { const n stations.length; const A []; const b []; // 构建矩阵 A 和向量 b for (let i 0; i n; i) { const row []; for (let j 0; j n; j) { const dx stations[i].lon - stations[j].lon; const dy stations[i].lat - stations[j].lat; const h Math.sqrt(dx * dx dy * dy); row.push(gamma(h, model)); // Γ 矩阵 } row.push(1); // 拉格朗日乘子列 A.push(row); // b 向量目标点到站点的变异函数值 const dx stations[i].lon - target.lon; const dy stations[i].lat - target.lat; const h Math.sqrt(dx * dx dy * dy); b.push(gamma(h, model)); } // 最后一行约束条件 const lastRow new Array(n).fill(1); lastRow.push(0); A.push(lastRow); b.push(1); // 高斯消元求解 return gaussianElimination(A, b).slice(0, n); } // 球状模型变异函数 function gamma(h, model) { if (h 0) return 0; if (h model.a) return model.C0 model.C; return model.C0 model.C * (1.5 * h / model.a - 0.5 * Math.pow(h / model.a, 3)); }这个方程组的规模是 (n1) x (n1)n 是站点数量。如果站点有几百个直接解方程没问题如果上千就需要用邻域搜索只选目标点附近的站点参与计算否则矩阵太大、解算太慢。高斯消元的时间复杂度是 O(n^3)所以邻域半径的选择很关键——既要保证精度又要控制计算量。3. 从插值结果到 Cesium 渲染栅格化与贴图的实现链路3.1 插值网格的生成策略克里金方程解出来的是单个点的温度值要渲染到 Cesium 地球上必须先把目标区域栅格化。常见做法是定义一个矩形范围按固定分辨率生成网格对每个网格点调用一次克里金解算。分辨率决定了渲染精度和计算量0.01 度大约对应 1 公里0.05 度大约 5 公里。// 生成插值网格 function generateGrid(bounds, resolution, stations, model) { const grid []; for (let lon bounds.west; lon bounds.east; lon resolution) { for (let lat bounds.south; lat bounds.north; lat resolution) { const target { lon, lat }; const weights solveKrigingWeights(stations, target, model); let value 0; for (let i 0; i stations.length; i) { value weights[i] * stations[i].value; } grid.push({ lon, lat, value }); } } return grid; }这里有个性能陷阱如果对每个网格点都重新解一次方程组计算量是网格数乘以站点数的三次方。实际项目中我一般会先用邻域搜索把参与解算的站点限制在 20 到 50 个这样单次解算的矩阵规模可控整体速度能提升一个数量级。3.2 栅格数据转 Cesium 纹理网格数据生成后需要转成 Cesium 能识别的纹理。最直接的方式是生成一张 Canvas把每个网格点的温度值映射成颜色再用 Cesium.SingleTileImageryProvider 或 Cesium.ImageryLayer 贴到地球上。// 网格转 Canvas 纹理 function gridToCanvas(grid, bounds, resolution) { const width Math.ceil((bounds.east - bounds.west) / resolution); const height Math.ceil((bounds.north - bounds.south) / resolution); const canvas document.createElement(canvas); canvas.width width; canvas.height height; const ctx canvas.getContext(2d); const imageData ctx.createImageData(width, height); grid.forEach(point { const x Math.floor((point.lon - bounds.west) / resolution); const y Math.floor((bounds.north - point.lat) / resolution); const index (y * width x) * 4; const color temperatureToColor(point.value); // 温度映射颜色 imageData.data[index] color.r; imageData.data[index 1] color.g; imageData.data[index 2] color.b; imageData.data[index 3] 180; // 透明度 }); ctx.putImageData(imageData, 0, 0); return canvas; }温度映射颜色通常用蓝到红的渐变但要注意温度范围可能因为插值出现异常值。比如站点最低温是 -10 度但插值结果可能出现 -15 度这时候颜色映射表要留出余量否则会出现颜色截断。我一般会把映射范围设为实际数据范围的 1.2 倍避免边界溢出。3.3 动态更新与性能取舍温度数据如果是实时更新的每次更新都重新解算整个网格会非常慢。常见的优化思路是只对变化区域重新插值或者降低非关注区域的分辨率。Cesium 的 ImageryLayer 支持动态替换纹理所以可以在后台线程算好新网格再一次性替换。另一个取舍是插值精度和渲染帧率。如果网格分辨率太高Canvas 尺寸会很大上传 GPU 的纹理也会很大低端设备可能卡顿。我一般会把网格分辨率控制在 256x256 到 512x512 之间既能保证视觉效果又不会给 GPU 太大压力。4. 避坑与排查克里金插值在 Cesium 里最容易翻车的几个点4.1 插值结果出现负温度或超高温现象温度场渲染出来某些区域出现 -50 度或者 60 度的异常值明显偏离实际。原因克里金插值在站点稀疏区域外推时权重可能变成负数导致插值结果超出原始数据的范围。这是克里金法的固有特性不是代码 bug。解决在插值后加一个裁剪步骤把结果限制在合理范围内。或者改用对数正态克里金它对正值数据更稳定。如果只是渲染用直接 clamp 到 [min, max] 也能接受。4.2 变异函数拟合不收敛现象拟合出来的变程特别小或者特别大插值结果要么全是噪点要么一片平滑。原因实验变异函数的分箱宽度选得不好或者站点分布极不均匀导致某些距离区间没有样本对。解决先画实验变异函数的散点图肉眼确认趋势。分箱宽度建议设为最大距离的 1/10 到 1/15。如果站点集中在城市可以先用网格重采样让站点分布均匀一些再拟合。4.3 Cesium 纹理贴图错位现象温度场整体偏移或者上下颠倒。原因Canvas 的 Y 轴方向和地理纬度的方向相反纬度越高 Y 越小但 Canvas 的 Y 是从上往下增加的。解决在计算 Y 坐标时用(bounds.north - lat) / resolution而不是(lat - bounds.south) / resolution。这个坑我踩过不止一次每次都要重新推导一遍坐标变换。4.4 站点数量多时页面卡死现象站点超过 500 个后浏览器直接无响应。原因每个网格点都解算全站点矩阵复杂度是 O(网格数 * 站点数^3)。解决引入邻域搜索每个网格点只选最近的 30 个站点参与解算。邻域半径可以用变程的 1.5 倍作为初始值再根据效果调整。另外可以把插值计算放到 Web Worker 里避免阻塞主线程。4.5 颜色映射表不生效现象温度值算出来了但渲染出来全是同一个颜色。原因温度值没有归一化或者颜色映射函数的输入范围不对。解决先统计插值结果的最小值和最大值用这个范围做归一化。如果数据范围很小比如 20 到 25 度要手动扩大映射范围否则颜色变化不明显。5. 设计取舍什么时候该用克里金什么时候该退回去克里金插值的计算成本比反距离权重高一个数量级这是它最大的代价。反距离权重只需要算距离克里金要拟合变异函数、解线性方程组。如果你的温度站点有上千个而且要求实时渲染克里金可能不是最优选择。我一般会先问三个问题站点数量多少、更新频率多高、精度要求多严。如果站点少于 200 个更新频率低于每分钟一次精度要求是气象级克里金值得上。如果站点上千、每秒都要刷新反距离权重或者径向基函数更合适。还有一种折中方案用克里金生成一张基准温度场后续更新只做局部修正这样既能保留克里金的空间结构又能控制计算量。另一个取舍是变异函数的拟合方式。自动拟合省事但遇到异常数据容易跑偏手动调参费时但可控性强。我一般会先用自动拟合出一组参数然后画实验变异函数和理论模型的对比图肉眼确认拟合质量。如果偏差太大就手动调整变程和拱高直到曲线贴合为止。验证插值效果的方法也很直接留出几个站点不参与插值用它们做交叉验证。计算预测值和实测值的均方根误差如果误差在可接受范围内说明变异函数拟合合理。这个方法比看渲染图靠谱得多因为渲染图的颜色映射会掩盖很多细节。从那以后我每次做克里金插值都强制走一遍交叉验证和变异函数拟合图检查不再凭感觉调参。希望帮到你。本文还有配套的精品资源点击获取