简介：Cesium 与克里金插值算法的组合，为离散温度观测点向连续温度场的三维渲染提供了完整落地路径；这份示例针对 GIS 前端与气象可视化开发者，逐步展示了从站点数据到平滑色带的关键流程。资源包共 585 个文件，大小约 7.27MB，其中 321 个 JS 脚本承担核心逻辑、29 个 JSON 保存数据与配置，141 张 PNG 和 53 张 JPG 提供贴图素材，27 个 CSS 负责界面样式，另有 SVG、WASM 等辅助资源。目前已有 167 人学习浏览，代码按初始化、插值、渲染和样式等模块组织，便于定位关键入口，也适合对照理解变异函数拟合、网格生成与颜色映射。读者可获得可直接运行的 Cesium 温度场示例，并可将插值思路扩展到降雨量、空气质量等连续面状数据的表达上。

1. 克里金插值遇上 Cesium：把离散温度计变成连续色斑的一次完整落地

做气象数据可视化的前端基本都遇到过同一个尴尬：地图上只有十几个自动气象站的温度读数，用户却追问「整个区域的温度分布长什么样」。单点图标说不清楚，空间插值才是答案。这个 Cesium Kriging 项目解决的正是这件事——把离散温度采样点通过克里金插值推算成连续色斑，渲染到三维地球上。项目最值得拆的地方在于：它没有依赖现成的插值库，而是把克里金从变异函数拟合、权重矩阵求解到网格渲染的整条链路在 Cesium 里完整实现了。适合想搞懂插值原理并落到前端渲染的人，也适合需要在气象、环境监测可视化里做类似功能的地图开发者。下面按「选型 → 拆解 → 排坑 → 改造」的顺序展开。

2. 技术选型与项目结构：插值算法、坐标投影与渲染方案的三层决策

2.1 克里金插值的算法选型：普通克里金为什么是默认答案

克里金最早是矿冶领域的地质统计学方法，核心假设是空间上距离越近的点，数值相关性越强。它跟反距离加权（IDW）最大的区别在于：IDW 的权重是人为指定的幂函数，而克里金会根据半变异函数（Semivariogram）从数据本身反推出「多近算近、多远算远」，并且能给出每个预测点的方差估计。在温度插值这个场景里，这个方差就是业务上常说的置信区间——站网稀疏的区域，插值结果方差大，预警阈值就得谨慎处理。

算法家族里常见三个变体。简单克里金假设均值已知；普通克里金假设均值未知但恒定，通过拉格朗日乘子在求解权重时同时估计均值；泛克里金处理带趋势场的空间数据。对气象温度场这种变化平缓但没有严格恒定均值的场景，普通克里金是工程上的默认答案。代码里只需要在矩阵里加一行约束（权重和为 1），复杂度增量几乎为零。这也是我第一次实现时踩过的坑——一开始图省事写了个自定义距离权重版本，结果站点稀疏区域插值结果和直觉完全对不上，后来老老实实换回普通克里金才稳定。

2.2 工程目录与数据流：从 JSON 文件到图层上屏的路径

项目结构刻意保持了「数据、算法、渲染」三层分离：

cesium-kriging-demo/ ├── index.html # Cesium 容器与页面入口 ├── config.js # 网格步长、色带、插值范围等配置 ├── src/ │ ├── data/ │ │ ├── stations.json # 温度采样点：lon/lat/temp │ │ └── loader.js # 数据解析与坐标预处理 │ ├── kriging/ │ │ ├── variogram.js # 半变异函数计算与模型拟合 │ │ ├── solver.js # 线性方程组求解（高斯消元） │ │ ├── kriging.js # 权重求解与插值入口 │ │ └── grid.js # 规则网格生成 │ └── render/ │ ├── colorRamp.js # 数值到颜色的映射 │ └── cesiumLayer.js # Cesium 图元渲染

数据流是单向的：loader 解析 JSON 并做坐标换算 → variogram 计算变异函数 → kriging 求解权重并逐点插值 → grid 生成网格 → cesiumLayer 渲染。每一步都是纯函数，单独传数据就能出结果，这给后面做自动化测试留了很大余地。我一般会在每层入口打一个console.debug输出中间结果，排查问题时能直接定位到是哪一环出了问题。

3. 核心模块拆解：从经纬度采样点到屏幕色斑的完整链路

3.1 数据预处理：距离计算必须统一投影基准

采样数据长这样：每个站点一个经纬度和一个温度值。看起来简单，但直接拿经纬度算距离会出问题——经度 1 度在赤道附近和在高纬度地区的实际地面距离差了一倍。所以 loader 里第一步就是把经纬度换算成 Web Mercator 平面坐标，后续所有距离计算都基于这个平面坐标。

// loader.js —— 采样数据解析与坐标预处理的常见写法 export function preprocessStations(rawData) { return rawData.map(item => ({ lon: item.lon, lat: item.lat, temp: item.temperature, // 距离计算统一用平面坐标，避免经纬度球面距离的畸变 x: webMercatorX(item.lon), y: webMercatorY(item.lat) })); } function webMercatorX(lon) { return lon * 20037508.34 / 180; } function webMercatorY(lat) { const rad = lat * Math.PI / 180; const merc = Math.log(Math.tan(Math.PI / 4 + rad / 2)); return 20037508.34 * merc / Math.PI; }

这里的 20037508.34 是 Web Mercator 投影在赤道处的半周长，属于标准常量。换算之后站点间的距离就可以直接Math.sqrt(dx*dx + dy*dy)，半变异函数和高斯消元里的所有距离计算都走这套平面坐标。在市级范围的数据集上，Web Mercator 的投影误差可以忽略；如果你做的是省级以上大范围数据，我一般会改用以区域中心为原点的等距方位投影，代码只需替换webMercatorX/Y两个函数即可。

3.2 半变异函数拟合与权重矩阵求解

克里金插值的核心是半变异函数。对每对采样点，计算它们的距离 h 和半方差 γ（等于两点温度差平方的一半）。把所有点对按距离分箱取平均，得到离散的变异函数点，再用球状模型去拟合。

// variogram.js —— 半变异函数计算与球状模型拟合 export function computeBinnedVariogram(stations, binCount = 10) { const pairs = []; for (let i = 0; i < stations.length; i++) { for (let j = i + 1; j < stations.length; j++) { const dx = stations[i].x - stations[j].x; const dy = stations[i].y - stations[j].y; const h = Math.sqrt(dx * dx + dy * dy); const gamma = 0.5 * Math.pow(stations[i].temp - stations[j].temp, 2); pairs.push({ h, gamma }); } } const maxH = Math.max(...pairs.map(p => p.h)); const binSize = maxH / binCount; const bins = Array.from({ length: binCount }, () => ({ distances: [], values: [] })); for (const p of pairs) { const idx = Math.min(Math.floor(p.h / binSize), binCount - 1); bins[idx].distances.push(p.h); bins[idx].values.push(p.gamma); } return bins .filter(b => b.distances.length > 0) .map(b => ({ h: b.distances.reduce((s, v) => s + v, 0) / b.distances.length, gamma: b.values.reduce((s, v) => s + v, 0) / b.values.length })); }

分箱数量是第一个需要调的参数。站点少（少于 20 个）时点对本来就不多，分 10 箱可能有些箱子是空的；我一般按站点数量动态决定，Math.max(5, Math.floor(n * (n - 1) / 2 / 20))这样保证每个箱子平均有 20 个点对才可信。

得到分箱数据后拟合球状模型：γ(h) = c0 + c·(1.5h/a − 0.5(h/a)³)，h 超过变程 a 时 γ 恒等于基台值 c0 + c。这里的 c0 叫块金值，代表测量误差和微观变异，a 叫变程，代表空间相关性的有效范围。

// 球状模型拟合：返回带 gamma(h) 函数的变异函数对象 export function fitSpherical(binned) { const maxH = Math.max(...binned.map(b => b.h)); const maxGamma = Math.max(...binned.map(b => b.gamma)); // 初始值：块金取最小箱均值，基台取最大箱均值，变程取最大距离 const base = { model: 'spherical', c0: binned[0].gamma, sill: maxGamma, range: maxH }; // 实际实现里用梯度下降对三个参数做微调，目标是让模型值与箱均值残差平方和最小 const fitted = gradientDescentFit(base, binned); return { ...fitted, gamma: function (h) { if (h === 0) return 0; if (h > fitted.range) return fitted.sill; const r = h / fitted.range; return fitted.c0 + (fitted.sill - fitted.c0) * (1.5 * r - 0.5 * r * r * r); } }; }

梯度下降在这个参数只有 3 个、目标函数是凸函数的问题上收敛很快，几十次迭代就够。如果你不想自己写拟合，直接把c0取分箱第一箱的均值、sill取最大箱均值、range取最大距离，效果也不会差太多——拟合只是让色斑边界更平滑。

接下来进入克里金求解。普通克里金要解一个 (n+1)×(n+1) 的线性方程组：左上角是采样点两两之间的变异函数值，最后一行和一列是拉格朗日乘子约束（权重和为 1），右侧向量是每个采样点到目标点的变异函数值。

// kriging.js —— 普通克里金主流程 export function krigingValue(stations, variogram, x, y) { const n = stations.length; // 扩展矩阵：最后一维是拉格朗日乘子 λ const A = Array.from({ length: n + 1 }, () => new Float64Array(n + 1)); const b = new Float64Array(n + 1); for (let i = 0; i < n; i++) { for (let j = 0; j < n; j++) { A[i][j] = variogram.gamma(distance(stations[i], stations[j])); } A[i][n] = 1; b[i] = variogram.gamma(distance(stations[i], { x, y })); } for (let j = 0; j < n; j++) A[n][j] = 1; b[n] = 1; const weights = gaussianSolve(A, b); let value = 0; for (let i = 0; i < n; i++) value += weights[i] * stations[i].temp; return { value, weights }; }

求解完成后，目标点的温度就是各采样点温度的加权和。注意权重可能为负，这不是 bug——克里金的权重不要求非负，负权重意味着该点对预测值起「修正」作用。高斯消元的时间复杂度是 O(n³)，站点数在 500 以内完全没问题；超过这个量级，我一般会换 Cholesky 分解，因为克里金矩阵是对称正定的（加块金后），分解只需 O(n³/3) 且数值更稳定。

3.3 网格生成与 Cesium 图元渲染

插值本身只针对单个点，要在 Cesium 上画出色斑，需要把目标区域切分成规则网格，对每个网格中心做一次krigingValue。

// grid.js —— 在指定范围生成规则网格并逐点插值 export function generateGrid(bounds, step, interpolateFn) { const cells = []; const lonCount = Math.ceil((bounds.east - bounds.west) / step); const latCount = Math.ceil((bounds.north - bounds.south) / step); for (let i = 0; i < lonCount; i++) { for (let j = 0; j < latCount; j++) { const lon = bounds.west + i * step; const lat = bounds.south + j * step; const result = interpolateFn(lon, lat); cells.push({ lon, lat, value: result.value }); } } return cells; }

step是网格步长，0.01 度大约是 1 公里。步长决定了两件事：色斑的细腻程度和计算量，网格数是步长的平方级增长，0.01 度在 1°×1° 范围内就是 100×100 = 10000 个插值点，每个点都要解一次 n+1 维线性方程组。所以插值函数要尽量轻量，变异函数对象的gamma方法里不要做对象分配，能用数字就用数字。

拿到网格后，最简单的渲染方案是每个网格单元建一个 CesiumRectangle实体：

// cesiumLayer.js —— 网格渲染为 Cesium 色块 export function renderTemperatureLayer(viewer, grid, colorStops, step) { const entities = []; for (const cell of grid) { entities.push(viewer.entities.add({ rectangle: { coordinates: Cesium.Rectangle.fromDegrees( cell.lon, cell.lat, cell.lon + step, cell.lat + step ), material: Cesium.Color.fromCssColorString(rgbColor(cell.value, colorStops)), classificationType: Cesium.ClassificationType.BOTH } })); } return entities; }

classificationType: BOTH让色斑既能贴在地形上也能盖在 3D Tiles 建筑上，这是温度可视化最常见的诉求。但这个方案有明显的天花板——网格超过 5000 个单元时，Entity 数量过多会导致相机拖动卡顿。项目里保留了第二种渲染方案：把整张网格画到离屏 Canvas 上，生成一张带地理范围的图片，用ImageryLayer直接贴到地球表面，切换开关在config.js里配。5000 以下用 Entity，以上用 Canvas 贴图。

4. 避坑指南：矩阵奇异、色斑错位与渲染卡顿的排查记录

这一章是我反复折腾出来的血泪经验，每条都按「现象 → 原因 → 解决」记录。

4.1 插值结果出现离谱值：权重矩阵病态

现象：个别网格点的插值温度突然冲到 200°C 或者变成负值，色斑上出现刺眼的异常块。

原因：采样点中存在距离极近甚至重复的站点，或者多个站点近似共线，导致克里金矩阵接近奇异，高斯消元解出的权重出现 -8、12 这样的极端值。

解决：预处理阶段做站点去重（距离小于 50 米的只保留一个），同时在矩阵对角线加一个极小的块金值，避免严格奇异：

// solver.js —— 高斯消元前对矩阵做稳定性处理 const NUGGET = 1e-8; for (let i = 0; i < n; i++) { A[i][i] += NUGGET; // 块金相当于允许测量误差，数字越小越接近原始方程 }

块金加得越大，插值结果越平滑，但会牺牲对采样点的拟合精度。1e-8 这个量级是我在十几个站点上测出来比较稳的值，站点多时可以适当调大。

4.2 色斑整体偏移，和站点图标对不上

现象：把渲染出来的色斑叠加到底图上，色斑的冷暖分区和实际站点位置明显错位，比如站点明明 20°C，所在的网格块却渲染成了 15°C 的颜色。

原因：这是最典型的坐标混用事故。变异函数和克里金矩阵里用的是平面坐标，但网格生成时用的还是经纬度，两者混在一起算距离，导致权重计算全部失真。

解决：网格点的位置也要先换算成平面坐标，插值完成后再用原始经纬度去建 Rectangle。我在generateGrid里把经纬度转平面坐标的步骤提前，保证所有进krigingValue的坐标都是同一个投影基准。

4.3 相机一拖动就卡成幻灯片

现象：插值结果是出来了，但每次视角旋转都要等两三秒才能跟上，内存占用持续上涨最终浏览器崩溃。

原因：网格 100×100 时生成了 1 万个 Entity，Cesium 的 Entity 系统在数量过万后调度开销非常大，这是 Entity API 本身的性能边界。

解决：把 Entity 方案切换成 Canvas 离屏渲染 +ImageryLayer贴图方案；网格步长从 0.01 度放宽到 0.02 度，视觉差异很小但计算量直接降到原来的 1/4。另外，给渲染函数加一个「数据没变化就不重算」的缓存，用网格哈希做键，避免同一份数据反复构建图层。

4.4 色斑出现明显分层断层

现象：温度从 20°C 到 21°C 之间颜色过渡不自然，能看出明显的条带状色阶，像等高线一样。

原因：色带函数用离散阈值分段，比如value < 18 ? '#00f' : value < 20 ? '#0af' : ...，造成相邻网格块之间颜色跳变。

解决：把色带定义成连续插值函数，在 RGB 空间里对相邻色标做线性插值，色斑过渡就自然了。我习惯在 HSL 空间插值，色相过渡比 RGB 线性插值视觉上更均匀，不会出现中间发灰的问题。

5. 二次开发：实时数据接入、色带配置与验证习惯

拿到这个项目后最值得改的第一个点，是把静态 JSON 换成实时数据流。温度观测一般是分钟级更新，做法是 WebSocket 推送新站点数组，到达后先做一次快速检查：所有站点温度均值和上轮相比变化小于 0.1°C 就直接跳过，否则重新走一遍「拟合 → 网格 → 渲染」流水线。

// 实时更新的常见写法：防抖 + 变化阈值双重过滤 ws.onmessage = debounce((payload) => { const stations = preprocessStations(payload.stations); if (meanTemp(stations) - lastMean < 0.1) return; const variogram = fitSpherical(computeBinnedVariogram(stations)); const grid = generateGrid(config.bounds, config.step, (lon, lat) => { const p = toMercator(lon, lat); return krigingValue(stations, variogram, p.x, p.y); }); viewer.imageryLayers.removeAll(); renderAsImagery(viewer, grid, config.colorStops); lastMean = meanTemp(stations); }, 1000);

第二个改造点是色带配置化。项目里温度映射目前是固定色标，改成读取config.js里的colorStops数组后，业务方调整配色就不用动代码了。验证方式我很看重——插值算法是黑匣子，出了问题不好定位，所以我给流水线写了一个最小的测试用例：取 5 个站点构建已知温度场，插值后断言目标点的值与手工解方程的值误差在 1e-6 以内，数据格式变了先跑一遍测试再上界面。

从那以后我每次接入新数据源，都强制走一遍 loader → variogram → grid → render 的单测流程，再放上界面看效果，省掉了至少一半的联调时间。希望帮到你。

本文还有配套的精品资源，点击获取