简介面向Web前端与三维可视化开发者的Cesium克里金插值示例包帮助理解如何在Cesium中结合geostat-js实现空间数据插值并生成动态3D点或数据表面。资源以HTML页面为入口搭配3个JavaScript库文件与1个GeoJSON数据文件另有1个RAR归档文件整体仅57KB结构精简适合已有Cesium基础、想深入Kriging算法应用的读者快速上手。资源包共6个文件浏览人数达1354人值得作为空间插值可视化的参考样板。通过阅读示例代码可掌握从准备经纬度与属性值数据、调用克里金模型预测、到在Cesium场景中渲染彩色点实体的完整链路还能了解不同变异性模型对插值精度和结果平滑程度的影响便于后续将这一能力扩展到地形拟合、环境监测、地质建模等真实场景。1. cesium克里金插值示例散点监测数据怎么变成三维热面拿到几十个离散站点的监测数据想做区域连片分析最常搜到的就是这个标题“cesium克里金插值示例html,三维开发实例 前端开发”。它背后其实是一整套前端三维开发写法在一个 HTML 页面里完成数据读取、空间插值、渲染上地球全程不依赖桌面 GIS。克里金是空间插值里最常用的地统计方法浏览器完全跑得动。适合气象、水文、环境监测、地质调查这类手握站点数据的从业者也适合做可视化大屏的前端开发工程师。这篇把从零到能跑的完整路径、参数设置和踩坑记录一次讲透你可以照着改成自己的数据。2. 克里金插值的前端路线先想清楚计算放在哪一侧2.1 为什么普通克里金能塞进浏览器克里金插值的核心不是画图而是求解。普通克里金Ordinary Kriging的计算可以拆成两个阶段第一步根据采样点对之间的距离和半方差拟合一个半变异函数模型第二步对每个待预测点用这个模型求一组权重再做加权平均。第二阶段本质上是解一个小规模线性方程组待预测点彼此独立天然适合循环。很多前端开发一听“插值”就觉得该上 Python实际上一百个站点、一百乘一百的网格浏览器里也只是毫秒级到秒级的事。浏览器端的克里金库虽然生态不大但核心算法是纯数学JavaScript 完全能承担。这里的关键是别把数据量和网格大小一起堆上去控制好规模纯前端方案完全可用。2.2 三条路线选型对比实时、预计算、混合预渲染做这个标题之前先回答一个问题克里金到底在哪一侧算完。我见过不少项目上来就写前端代码结果数据一到几千个站点就卡死最后灰溜溜地加后端。按数据规模选路线比选库更重要。路线计算位置典型规模优点缺点纯前端实时计算浏览器站点数 500网格 200×200一个 HTML 全搞定部署就是静态文件站点上千后计算时间明显拉长后端预计算Python / R几万站点网格任意地统计方法全支持交叉验证要维护服务前端只做展示混合预渲染后端算完存图片/GeoJSON海量数据 大屏前端零计算加载即显示更新数据要重跑后端任务我一般会这样判断如果是几十个站点做演示或内部工具直接走纯前端。如果数据量超过一千个站点或者要做时间序列逐小时更新就提前用 pykrige 或 gstat 把插值结果导出成网格图片前端只负责贴图。克里金本身是有代价的别为了“纯前端”三个字硬扛。2.3 kriging.js 的模型与必调参数纯前端路线里社区流传最广的库是 kriging.js一个单文件实现的普通克里金库。它提供的接口很精简train负责拟合半变异函数predict负责对单个点做预测grid可以生成网格面片但性能一般我习惯自己写循环。它的核心调用是kriging.train(values, xs, ys, model, sigma2, alpha)其中model是半变异函数模型支持exponential、gaussian、spherical三选一。sigma2是平滑噪声项默认 1/3值越大插值面越平滑但会牺牲局部细节。alpha是半变异函数的形状参数和数据的空间跨度强相关我通常先按数据跨度的 1/3 到 1/2 给初值再根据输出微调。这里有一个大多数教程不会说的关键点train里传入的x、y不能直接用经纬度。半变异函数计算的是点的距离经纬度差是“度”不是“米”一个站点跨度 10 度的数据直接算出来的变程是个奇怪的大数插值面完全失真。正确做法是先把经纬度投影成平面坐标算完再反投影回经纬度给 Cesium 用。第三节会给出具体的投影函数。3. 把站点数据装进 HTMLCesium 页面骨架与坐标的坑3.1 一个能直接跑的 Cesium 初始化骨架标题里带“html”三个字意味着这个项目应该能在一个 HTML 文件里跑起来。Cesium 的初始化并不复杂但默认的 Viewer 会在界面上铺满各种控件做数据展示时通常只保留必要的地球操作。!DOCTYPE html html langzh-cn head meta charsetutf-8 titlecesium克里金插值示例/title style html, body, #cesiumContainer { width: 100%; height: 100%; margin: 0; padding: 0; overflow: hidden; } /style !-- 替换成你项目里的 Cesium 构建路径这里用官网 latest 入口 -- link hrefhttps://cesium.com/downloads/cesiumjs/releases/latest/Build/Cesium/Widgets/widgets.css relstylesheet script srchttps://cesium.com/downloads/cesiumjs/releases/latest/Build/Cesium/Cesium.js/script /head body div idcesiumContainer/div script Cesium.Ion.defaultAccessToken 你的IonToken; // 可选没有就注释掉这行 const viewer new Cesium.Viewer(cesiumContainer, { baseLayerPicker: false, // 关掉底图切换器 geocoder: false, // 关掉搜索框 homeButton: false, // 关掉回到默认视角按钮 animation: false, // 关掉动画控件 timeline: false, // 关掉时间轴 infoBox: false, // 关掉点击弹窗 sceneMode: Cesium.SceneMode.SCENE3D, shouldAnimate: true }); // 去掉默认的版权信息不是必须的但做演示时可以减少干扰 viewer.cesiumWidget.creditContainer.style.display none; /script /body /html这段初始化里值得说的是shouldAnimate: true。Cesium 的时钟默认走的是仿真时间后面做时间序列轮播时这个开关必须打开。sceneMode固定为三维场景因为克里金插值面要贴在三维地形或地球表面上。如果你接的是静态页面没有 Ion Token 也能用默认的椭球体底图只是没有影像图层不影响插值结果展示。3.2 解析 CSV 站点数据并做平面投影站点数据最常见的格式是 CSV字段是站点名、经度、纬度、监测值。单文件 HTML 里直接用fetch读本地文件会撞上 file:// 协议的跨域限制所以我一般把数据直接内联成字符串既省去起服务的麻烦也方便拷贝给别人。// 内联 CSV 数据id, 经度, 纬度, 监测值 const csvText A01,116.391,39.907,12 A02,116.405,39.915,18 A03,116.380,39.890,9 A04,116.420,39.920,22 A05,116.360,39.930,15 ; // 解析 CSV 成对象数组 const stations []; csvText.trim().split(\n).forEach(line { const [id, lon, lat, value] line.split(,).map(s s.trim()); stations.push({ id, lon: parseFloat(lon), lat: parseFloat(lat), value: parseFloat(value) }); }); // 以数据中心纬度为基准做等距圆柱近似投影 // 目的是把经纬度“度”转成“米”让克里金的距离计算有意义 const centerLat stations.reduce((s, p) s p.lat, 0) / stations.length; const METER_PER_DEGREE 111319.9; function projectLonLat(lon, lat) { return [ lon * Math.cos(centerLat * Math.PI / 180) * METER_PER_DEGREE, lat * METER_PER_DEGREE ]; } // 给每个站点追加投影坐标 stations.forEach(s { const [x, y] projectLonLat(s.lon, s.lat); s.x x; s.y y; });投影这段是克里金能否正确出图的分水岭。METER_PER_DEGREE是赤道上 1 度纬度的米数乘以cos(centerLat)是修正纬线收缩。这个近似在城市级、省级范围内精度足够克里金的半变异函数要求的是相对距离而不是绝对精度。如果你处理的是全国范围的数据建议用正经的墨卡托投影或换成后端计算前端近似会很勉强。3.3 用 entity 画点验证数据位置开始插值之前先把站点画到地球上确认位置和数值范围这一步能筛掉大量“数据源本身有问题”的翻车现场。Cesium 里画点最简单的方式是用entity加point它适合少量标点数量超过几千个再考虑primitive。stations.forEach(s { viewer.entities.add({ id: station_${s.id}, position: Cesium.Cartesian3.fromDegrees(s.lon, s.lat), point: { pixelSize: 8, color: Cesium.Color.fromCssColorString(#ff5722), outlineColor: Cesium.Color.WHITE, outlineWidth: 1, disableDepthTestDistance: Number.POSITIVE_INFINITY }, description: b${s.id}/bbr数值${s.value} }); }); // 视角定位到数据范围方便快速检查 const positions stations.map(s Cesium.Cartesian3.fromDegrees(s.lon, s.lat)); viewer.camera.flyTo({ destination: Cesium.Rectangle.fromCartesianArray(positions), duration: 1.5 });disableDepthTestDistance设为正无穷是让点位不被地形遮挡否则点会陷进山体里看不到。flyTo定位到数据范围后旋转地球从侧面检查一次点位的空间分布如果发现某个点明显偏离预期区域先修数据再继续不要在脏数据上做插值。4. 克里金插值出图网格计算、canvas 上色、贴到三维地球4.1 用 kriging.js 训练半变异函数并预测网格kriging.js 的用法非常直接train一次训练模型然后对每个网格点做predict。但要看懂这里的性能账每个predict都要跟所有站点做一次距离计算和解小矩阵所以网格数乘站点数就是总计算量。一百个站点、一百乘一百网格是一百万次乘加浏览器能扛五百个站点、两百乘两百网格就跨过秒级门槛了。// 引入 kriging.js 后从站点数据里拆出数组 const xs stations.map(s s.x); const ys stations.map(s s.y); const values stations.map(s s.value); // 训练克里金模型模型选 exponential平滑参数用 0.5alpha 按数据跨度给初值 const variogram kriging.train(values, xs, ys, exponential, 0.5, getAlpha()); // 生成网格范围投影坐标下的边界 const xMin Math.min(...xs), xMax Math.max(...xs); const yMin Math.min(...ys), yMax Math.max(...ys); // 网格分辨率80x80 起步数据点多时降到 60演示足够 const GRID_W 80, GRID_H 80; const grid []; // grid[row][col] 存预测值 for (let row 0; row GRID_H; row) { grid[row] []; const lat yMin (row / GRID_H) * (yMax - yMin); for (let col 0; col GRID_W; col) { const lon xMin (col / GRID_W) * (xMax - xMin); const v kriging.predict(lon, lat, variogram); grid[row][col] v; } }kriging.train的sigma2我给了 0.5这比库默认的 1/3 略大面会更平滑适合站点只有十几个的场景。getAlpha()不是库内置函数而是根据数据跨度算出的初值通常取投影坐标下数据范围对角线长度的三分之一。这个参数没有银弹我一般写完第一版后打印几个已知站点的预测值和真实值做对比偏差大就调alpha这是最直接的调参路径。4.2 网格数据画成带透明度的 canvas 色块网格数据本身不是图要变成能让 Cesium 消费的纹理最省事的方式是把它画成一张 canvas。canvas 的每个像素对应网格的一个单元格颜色根据插值结果映射透明度留给 Cesium 的材质去控制。function drawGridToCanvas(grid, colsNumber, rowsNumber) { const canvas document.createElement(canvas); canvas.width colsNumber; canvas.height rowsNumber; const ctx canvas.getContext(2d); const imageData ctx.createImageData(colsNumber, rowsNumber); // 先扫一遍网格拿到数值范围 let vMin Infinity, vMax -Infinity; for (let y 0; y rowsNumber; y) { for (let x 0; x colsNumber; x) { const v grid[y][x]; if (v vMin) vMin v; if (v vMax) vMax v; } } // 预生成 256 级颜色映射表蓝色到红色渐变 const lut buildColorLUT(); for (let y 0; y rowsNumber; y) { for (let x 0; x colsNumber; x) { const t (grid[y][x] - vMin) / (vMax - vMin); // 0~1 归一化 const color lut[Math.min(255, Math.max(0, Math.round(t * 255)))]; const idx (y * colsNumber x) * 4; imageData.data[idx] color[0]; // R imageData.data[idx 1] color[1]; // G imageData.data[idx 2] color[2]; // B imageData.data[idx 3] 200; // A200 让底图透出来 } } ctx.putImageData(imageData, 0, 0); return canvas; } // 简单的蓝-青-黄-红渐变实际项目里可以换成自己的配色 function buildColorLUT() { const lut []; for (let i 0; i 256; i) { const t i / 255; lut.push([ Math.round(255 * Math.min(1, Math.max(0, (t - 0.33) * 3))), Math.round(255 * Math.min(1, Math.max(0, 1 - Math.abs(t - 0.33) * 3))), Math.round(255 * Math.min(1, Math.max(0, (0.66 - t) * 3))) ]); } return lut; }这里的颜色映射用的是线性归一化数据分布偏斜时容易一层色覆盖大面积区域第五节会讲怎么用分位数拉伸替代。网格和 canvas 像素一一对应putImageData效率远高于逐像素fillRect前者是直接写内存后者要跑绘制管线。地图方向在这里是“上北下南”Cesium 贴图时如果方向反了在后面贴图代码里通过缩放负值修正。4.3 贴到地球entity 的 Rectangle 方案与 primitive 的区别canvas 生成之后接下来的问题是“怎么贴到三维地球上”。常见做法是加一个entity用rectangle指定经纬度范围material用ImageMaterialProperty引用 canvas。这个方案声明式、代码少改颜色换图都方便适合插值面这种一张图带整个区域的数据。const krigingCanvas drawGridToCanvas(grid, GRID_W, GRID_H); // 把网格范围投影坐标反投影回经纬度 function unprojectX(x) { return x / (Math.cos(centerLat * Math.PI / 180) * METER_PER_DEGREE); } function unprojectY(y) { return y / METER_PER_DEGREE; } const rectangleEntity viewer.entities.add({ id: krigingSurface, rectangle: { coordinates: Cesium.Rectangle.fromDegrees( unprojectX(xMin), unprojectY(yMin), unprojectX(xMax), unprojectY(yMax) ), material: new Cesium.ImageMaterialProperty({ image: krigingCanvas, transparent: true }), clampToGround: true } });clampToGround: true让矩形贴到地形表面而不是悬在椭球面上如果你的场景加载了真实地形这行是必须的。如果发现贴图方向不对用 Cesium 的坐标系翻转可以解决但更常见的坑在 5.4 节。Cesium 里加载这种面数据entity 和 primitive 的区别值得说清楚。entity 是声明式的适合“我要一个贴地矩形”这种需求代码量少内部帮你处理状态变更primitive 是命令式的需要你手动构造 Geometry、Appearance、材质好处是渲染性能高、可控性强适合大量动态图元。克里金插值面通常只有一张图entity 完全足够如果你要做几十个时刻的插值面轮播且每个都是几百 KB 的 canvasprimitive 配合批量刷新会更稳。另一个区别是 entity 的description、id这些属性天然支持拾取和弹窗这在调试时非常方便。4.4 想要等值线用 turf.isolines 直接抽线插值面色块图之外业务上经常还要等值线比如等降水量线、等浓度线。网格数据已经算好了等值线就不需要再回克里金直接基于网格抽线即可。turf.js 的isolines方法是这里最顺手的工具它内部会基于点集生成三角网再按给定分级值抽取等值线。// 把网格转成 turf 点集每个点带 value 属性 const features []; for (let row 0; row GRID_H; row) { for (let col 0; col GRID_W; col) { features.push(turf.point( [unprojectX(xMin (col / GRID_W) * (xMax - xMin)), unprojectY(yMin (row / GRID_H) * (yMax - yMin))], { value: grid[row][col] } )); } } const pointFC turf.featureCollection(features); // 按值域分 6 级抽等值线 const breaks [5, 10, 15, 20, 25]; const isolines turf.isolines(pointFC, breaks, { zProperty: value }); // 加载到 Cesium每条线画成 polyline isolines.features.forEach(line { const positions line.geometry.coordinates[0].map(coord Cesium.Cartesian3.fromDegrees(coord[0], coord[1]) ); viewer.entities.add({ polyline: { positions: positions, width: 2, material: Cesium.Color.fromCssColorString(#333333).withAlpha(0.8) } }); });turf.isolines的breaks参数是分级边界数组有多少个边界就会抽出多少条线边界值要落在实测数据的值域里才有意义。另一个需要注意的点是zProperty必须和写入点的属性名一致写成value是常见做法。等值线抽出来之后配合 4.2 的色块图一起显示专业感一下就上来了这也是很多环境监测大屏的标准呈现方式。5. cesium克里金插值避坑五个翻车现场的现象、原因与解决5.1 插值面在数据边界外“翘边”现象插值面超出了站点分布范围在凸包外围出现一圈明显的高值或低值“帽子”看起来像面被风吹起来一块。原因普通克里金对预测点的外推能力很弱超出数据覆盖范围后半变异函数主导了权重分配预测值会逐渐回归到数据均值附近。这个回归过程不一定是单调的配合边缘处站点稀疏就会在凸包外形成翘边假象。解决用凸包把插值范围裁掉。常见做法是用 turf.js 的convex生成站点凸包多边形然后对网格点做 point-in-polygon 剔除。更平滑的效果用turf.concave但 concave 的参数敏感性高数据分布不规则时容易生成奇怪的形状。生产项目里用凸包已经足够把网格点过滤掉之后再画 canvas天然就没有翘边区域。5.2 经纬度当平面坐标直接算插值彻底失真现象train的参数直接传经纬度插值面要么是一块块色斑要么整体平滑到看不出梯度少数站点旁边出现不该有的高值孤岛。原因半变异函数内部按距离建模经纬度差到米之间存在非线性关系。北大荒一个站点跨度 10 度距离是“度数”而数值可能是几十两者量纲完全不匹配拟合出来的变程和基台值没有物理意义。解决按 3.2 的投影函数处理把经纬度转成以中心点为基准的米制平面坐标。算完预测值之后渲染 Cesium 时再通过反投影函数把坐标还原成经纬度。这里注意反投影函数要和正投影一一对应我在 4.3 里特意写了unprojectX/unprojectY目的就是保证投影前后能对得上。5.3 网格分辨率一上 150页面直接假死现象数据量和分辨率明明不大但一执行到插值循环页面白屏好几秒风扇狂转控制台还报阻塞警告。原因predict的循环是双重嵌套网格 150×150 是 22500 个预测点每个预测点要和全部站点做矩阵运算。如果站点数再是 300就是 675 万次距离计算和解小矩阵主线程被占满渲染帧率掉到 0。解决两件事。第一网格分辨率控制在 100×100 以内插值面展示的视觉信息量在 80×80 已经足够细节更多的地图操作应该交给图片缩放而不是增加计算点。第二如果确实需要高分辨率把预测计算挪到 Web Worker 里主线程只接收结果。这里有个血泪经验不要在viewer.scene.preUpdate里触发计算每一帧都会重算一次正确做法是初始化时算一次把结果缓存起来。5.4 贴地图片跟地形各算各的时不时报 renderError现象控制台出现viewer.scene.renderError贴地矩形在地形起伏区域闪边、花屏甚至整块消失换个视角又恢复。原因clampToGround的贴地方式是把几何体压到地形高度上对高程数据的分辨率非常敏感。默认的椭球体底图没有高程一切正常一旦加载了真实地形数据矩形边界的四个角按地形高度拉伸中间像素却按纹理 UV 线性映射跟实际地形起伏对不上就会出现撕裂感。另外canvas 尺寸过大也会触发纹理上限浏览器和显卡对单边像素上限普遍是 4096超出后纹理被截断。解决第一canvas 尺寸控制在 1024×1024 以下用的时候让它做双线性缩放视觉损失几乎看不出来。第二项目确实要用真实地形时先用 cesium terrain builder 之类工具把地形切好再决定贴图方式不要在运行时既开地形又挂 clampToGround 的图片面。我排查这类问题时的习惯是先用viewer.scene.globe.depthTestAgainstTerrain false做临时对比看是不是深度测试导致的闪边。5.5 色带断层严重低值区被一种颜色吞掉现象插值面大片区域是同一种颜色只有少数高值点周围有渐变色带明显分层像用了 8 位色一样。原因线性归一化遇到偏斜分布时的问题。监测数据大多服从对数正态或偏态分布低值数据密集、高值数据稀疏线性映射把 90% 的网格像素挤到色带前 10% 的区间里肉眼看起来就是一条色块。解决用分位数拉伸替代线性归一化。先对全部网格预测值排序取 5% 和 95% 分位值作为映射的上下界超出部分截断。这样中间 90% 的数据能铺满整个色带分布细节完全打开。这个技巧几乎适用所有克里金插值的可视化也是前端开发工程师面试题里“数据可视化如何避免颜色失真”的一个标准答案我每次带新人都会强调这一点。6. 从静态出图到动态序列让克里金插值动起来6.1 预渲染多时次图片用 Cesium clock 驱动轮播克里金插值最常见的业务形态不是单一时次而是逐小时、逐日更新。如果你在用户拖动时间轴时才现场算插值性能一定扛不住。我的做法是预渲染把所有时次的插值结果提前算好每个时次生成一张 canvas存成数组运行时只做图片替换。// 假设 timeSeries 是多个时次的预测结果数组 const canvasList timeSeries.map(frame { const c drawGridToCanvas(frame.grid, GRID_W, GRID_H); return c; }); // 用 viewer.clock 驱动轮播 const surfaceEntity viewer.entities.getById(krigingSurface); let currentIndex 0; viewer.clock.onTick.addEventListener(() { const second Math.floor(viewer.clock.currentTime.secondsOfDay); const index second % canvasList.length; if (index ! currentIndex) { currentIndex index; surfaceEntity.rectangle.material.image canvasList[index]; } });注意surfaceEntity.rectangle.material.image直接赋 canvas 对象Cesium 会把它当作纹理更新不需要重新创建 entity。这样切换的粒度是“一张图换另一张图”主线程只有一次纹理上传性能开销远小于重新做克里金预测。计算全部放在初始化阶段把阻塞留在一个可接受的加载时长里交互时保持流畅。6.2 验证插值结果的三个自检习惯插值面画出来之后我一般做三个自检都通过才会说这个面“能用”。第一站点值回验。把站点本身的点画在色块图上层看站点位置的颜色是否和它自身的数值颜色一致。如果偏差超过色带的一到两个等级说明sigma2或alpha参数需要回调。第二留一交叉验证。随机剔除 10% 的站点不参与训练用模型预测这些被剔除的点计算预测值和真实值的均方根误差。这个数值能直接告诉你插值结果的可靠程度比肉眼看图客观得多。第三侧面视角检查。把相机压到接近地面贴着插值面边缘看过去。翘边、闪边、颜色断层在这种视角下都会现原形。这也是我为什么习惯在项目里保留一个viewer.camera.flyTo的调试按钮一键切换俯视和侧视检查效率非常高。这套流程走下来静态单时次和动态序列两个方向都能落地。做这类三维开发最忌讳的就是把克里金当成黑匣子参数乱给、坐标不管、结果不看就往上贴。每次调参前先问自己这一步的数据是什么单位、在哪个坐标系、预期范围是多少。本文还有配套的精品资源点击获取 Cesium克里金插值实战：从散点到三维热面的前端实现
简介面向Web前端与三维可视化开发者的Cesium克里金插值示例包帮助理解如何在Cesium中结合geostat-js实现空间数据插值并生成动态3D点或数据表面。资源以HTML页面为入口搭配3个JavaScript库文件与1个GeoJSON数据文件另有1个RAR归档文件整体仅57KB结构精简适合已有Cesium基础、想深入Kriging算法应用的读者快速上手。资源包共6个文件浏览人数达1354人值得作为空间插值可视化的参考样板。通过阅读示例代码可掌握从准备经纬度与属性值数据、调用克里金模型预测、到在Cesium场景中渲染彩色点实体的完整链路还能了解不同变异性模型对插值精度和结果平滑程度的影响便于后续将这一能力扩展到地形拟合、环境监测、地质建模等真实场景。1. cesium克里金插值示例散点监测数据怎么变成三维热面拿到几十个离散站点的监测数据想做区域连片分析最常搜到的就是这个标题“cesium克里金插值示例html,三维开发实例 前端开发”。它背后其实是一整套前端三维开发写法在一个 HTML 页面里完成数据读取、空间插值、渲染上地球全程不依赖桌面 GIS。克里金是空间插值里最常用的地统计方法浏览器完全跑得动。适合气象、水文、环境监测、地质调查这类手握站点数据的从业者也适合做可视化大屏的前端开发工程师。这篇把从零到能跑的完整路径、参数设置和踩坑记录一次讲透你可以照着改成自己的数据。2. 克里金插值的前端路线先想清楚计算放在哪一侧2.1 为什么普通克里金能塞进浏览器克里金插值的核心不是画图而是求解。普通克里金Ordinary Kriging的计算可以拆成两个阶段第一步根据采样点对之间的距离和半方差拟合一个半变异函数模型第二步对每个待预测点用这个模型求一组权重再做加权平均。第二阶段本质上是解一个小规模线性方程组待预测点彼此独立天然适合循环。很多前端开发一听“插值”就觉得该上 Python实际上一百个站点、一百乘一百的网格浏览器里也只是毫秒级到秒级的事。浏览器端的克里金库虽然生态不大但核心算法是纯数学JavaScript 完全能承担。这里的关键是别把数据量和网格大小一起堆上去控制好规模纯前端方案完全可用。2.2 三条路线选型对比实时、预计算、混合预渲染做这个标题之前先回答一个问题克里金到底在哪一侧算完。我见过不少项目上来就写前端代码结果数据一到几千个站点就卡死最后灰溜溜地加后端。按数据规模选路线比选库更重要。路线计算位置典型规模优点缺点纯前端实时计算浏览器站点数 500网格 200×200一个 HTML 全搞定部署就是静态文件站点上千后计算时间明显拉长后端预计算Python / R几万站点网格任意地统计方法全支持交叉验证要维护服务前端只做展示混合预渲染后端算完存图片/GeoJSON海量数据 大屏前端零计算加载即显示更新数据要重跑后端任务我一般会这样判断如果是几十个站点做演示或内部工具直接走纯前端。如果数据量超过一千个站点或者要做时间序列逐小时更新就提前用 pykrige 或 gstat 把插值结果导出成网格图片前端只负责贴图。克里金本身是有代价的别为了“纯前端”三个字硬扛。2.3 kriging.js 的模型与必调参数纯前端路线里社区流传最广的库是 kriging.js一个单文件实现的普通克里金库。它提供的接口很精简train负责拟合半变异函数predict负责对单个点做预测grid可以生成网格面片但性能一般我习惯自己写循环。它的核心调用是kriging.train(values, xs, ys, model, sigma2, alpha)其中model是半变异函数模型支持exponential、gaussian、spherical三选一。sigma2是平滑噪声项默认 1/3值越大插值面越平滑但会牺牲局部细节。alpha是半变异函数的形状参数和数据的空间跨度强相关我通常先按数据跨度的 1/3 到 1/2 给初值再根据输出微调。这里有一个大多数教程不会说的关键点train里传入的x、y不能直接用经纬度。半变异函数计算的是点的距离经纬度差是“度”不是“米”一个站点跨度 10 度的数据直接算出来的变程是个奇怪的大数插值面完全失真。正确做法是先把经纬度投影成平面坐标算完再反投影回经纬度给 Cesium 用。第三节会给出具体的投影函数。3. 把站点数据装进 HTMLCesium 页面骨架与坐标的坑3.1 一个能直接跑的 Cesium 初始化骨架标题里带“html”三个字意味着这个项目应该能在一个 HTML 文件里跑起来。Cesium 的初始化并不复杂但默认的 Viewer 会在界面上铺满各种控件做数据展示时通常只保留必要的地球操作。!DOCTYPE html html langzh-cn head meta charsetutf-8 titlecesium克里金插值示例/title style html, body, #cesiumContainer { width: 100%; height: 100%; margin: 0; padding: 0; overflow: hidden; } /style !-- 替换成你项目里的 Cesium 构建路径这里用官网 latest 入口 -- link hrefhttps://cesium.com/downloads/cesiumjs/releases/latest/Build/Cesium/Widgets/widgets.css relstylesheet script srchttps://cesium.com/downloads/cesiumjs/releases/latest/Build/Cesium/Cesium.js/script /head body div idcesiumContainer/div script Cesium.Ion.defaultAccessToken 你的IonToken; // 可选没有就注释掉这行 const viewer new Cesium.Viewer(cesiumContainer, { baseLayerPicker: false, // 关掉底图切换器 geocoder: false, // 关掉搜索框 homeButton: false, // 关掉回到默认视角按钮 animation: false, // 关掉动画控件 timeline: false, // 关掉时间轴 infoBox: false, // 关掉点击弹窗 sceneMode: Cesium.SceneMode.SCENE3D, shouldAnimate: true }); // 去掉默认的版权信息不是必须的但做演示时可以减少干扰 viewer.cesiumWidget.creditContainer.style.display none; /script /body /html这段初始化里值得说的是shouldAnimate: true。Cesium 的时钟默认走的是仿真时间后面做时间序列轮播时这个开关必须打开。sceneMode固定为三维场景因为克里金插值面要贴在三维地形或地球表面上。如果你接的是静态页面没有 Ion Token 也能用默认的椭球体底图只是没有影像图层不影响插值结果展示。3.2 解析 CSV 站点数据并做平面投影站点数据最常见的格式是 CSV字段是站点名、经度、纬度、监测值。单文件 HTML 里直接用fetch读本地文件会撞上 file:// 协议的跨域限制所以我一般把数据直接内联成字符串既省去起服务的麻烦也方便拷贝给别人。// 内联 CSV 数据id, 经度, 纬度, 监测值 const csvText A01,116.391,39.907,12 A02,116.405,39.915,18 A03,116.380,39.890,9 A04,116.420,39.920,22 A05,116.360,39.930,15 ; // 解析 CSV 成对象数组 const stations []; csvText.trim().split(\n).forEach(line { const [id, lon, lat, value] line.split(,).map(s s.trim()); stations.push({ id, lon: parseFloat(lon), lat: parseFloat(lat), value: parseFloat(value) }); }); // 以数据中心纬度为基准做等距圆柱近似投影 // 目的是把经纬度“度”转成“米”让克里金的距离计算有意义 const centerLat stations.reduce((s, p) s p.lat, 0) / stations.length; const METER_PER_DEGREE 111319.9; function projectLonLat(lon, lat) { return [ lon * Math.cos(centerLat * Math.PI / 180) * METER_PER_DEGREE, lat * METER_PER_DEGREE ]; } // 给每个站点追加投影坐标 stations.forEach(s { const [x, y] projectLonLat(s.lon, s.lat); s.x x; s.y y; });投影这段是克里金能否正确出图的分水岭。METER_PER_DEGREE是赤道上 1 度纬度的米数乘以cos(centerLat)是修正纬线收缩。这个近似在城市级、省级范围内精度足够克里金的半变异函数要求的是相对距离而不是绝对精度。如果你处理的是全国范围的数据建议用正经的墨卡托投影或换成后端计算前端近似会很勉强。3.3 用 entity 画点验证数据位置开始插值之前先把站点画到地球上确认位置和数值范围这一步能筛掉大量“数据源本身有问题”的翻车现场。Cesium 里画点最简单的方式是用entity加point它适合少量标点数量超过几千个再考虑primitive。stations.forEach(s { viewer.entities.add({ id: station_${s.id}, position: Cesium.Cartesian3.fromDegrees(s.lon, s.lat), point: { pixelSize: 8, color: Cesium.Color.fromCssColorString(#ff5722), outlineColor: Cesium.Color.WHITE, outlineWidth: 1, disableDepthTestDistance: Number.POSITIVE_INFINITY }, description: b${s.id}/bbr数值${s.value} }); }); // 视角定位到数据范围方便快速检查 const positions stations.map(s Cesium.Cartesian3.fromDegrees(s.lon, s.lat)); viewer.camera.flyTo({ destination: Cesium.Rectangle.fromCartesianArray(positions), duration: 1.5 });disableDepthTestDistance设为正无穷是让点位不被地形遮挡否则点会陷进山体里看不到。flyTo定位到数据范围后旋转地球从侧面检查一次点位的空间分布如果发现某个点明显偏离预期区域先修数据再继续不要在脏数据上做插值。4. 克里金插值出图网格计算、canvas 上色、贴到三维地球4.1 用 kriging.js 训练半变异函数并预测网格kriging.js 的用法非常直接train一次训练模型然后对每个网格点做predict。但要看懂这里的性能账每个predict都要跟所有站点做一次距离计算和解小矩阵所以网格数乘站点数就是总计算量。一百个站点、一百乘一百网格是一百万次乘加浏览器能扛五百个站点、两百乘两百网格就跨过秒级门槛了。// 引入 kriging.js 后从站点数据里拆出数组 const xs stations.map(s s.x); const ys stations.map(s s.y); const values stations.map(s s.value); // 训练克里金模型模型选 exponential平滑参数用 0.5alpha 按数据跨度给初值 const variogram kriging.train(values, xs, ys, exponential, 0.5, getAlpha()); // 生成网格范围投影坐标下的边界 const xMin Math.min(...xs), xMax Math.max(...xs); const yMin Math.min(...ys), yMax Math.max(...ys); // 网格分辨率80x80 起步数据点多时降到 60演示足够 const GRID_W 80, GRID_H 80; const grid []; // grid[row][col] 存预测值 for (let row 0; row GRID_H; row) { grid[row] []; const lat yMin (row / GRID_H) * (yMax - yMin); for (let col 0; col GRID_W; col) { const lon xMin (col / GRID_W) * (xMax - xMin); const v kriging.predict(lon, lat, variogram); grid[row][col] v; } }kriging.train的sigma2我给了 0.5这比库默认的 1/3 略大面会更平滑适合站点只有十几个的场景。getAlpha()不是库内置函数而是根据数据跨度算出的初值通常取投影坐标下数据范围对角线长度的三分之一。这个参数没有银弹我一般写完第一版后打印几个已知站点的预测值和真实值做对比偏差大就调alpha这是最直接的调参路径。4.2 网格数据画成带透明度的 canvas 色块网格数据本身不是图要变成能让 Cesium 消费的纹理最省事的方式是把它画成一张 canvas。canvas 的每个像素对应网格的一个单元格颜色根据插值结果映射透明度留给 Cesium 的材质去控制。function drawGridToCanvas(grid, colsNumber, rowsNumber) { const canvas document.createElement(canvas); canvas.width colsNumber; canvas.height rowsNumber; const ctx canvas.getContext(2d); const imageData ctx.createImageData(colsNumber, rowsNumber); // 先扫一遍网格拿到数值范围 let vMin Infinity, vMax -Infinity; for (let y 0; y rowsNumber; y) { for (let x 0; x colsNumber; x) { const v grid[y][x]; if (v vMin) vMin v; if (v vMax) vMax v; } } // 预生成 256 级颜色映射表蓝色到红色渐变 const lut buildColorLUT(); for (let y 0; y rowsNumber; y) { for (let x 0; x colsNumber; x) { const t (grid[y][x] - vMin) / (vMax - vMin); // 0~1 归一化 const color lut[Math.min(255, Math.max(0, Math.round(t * 255)))]; const idx (y * colsNumber x) * 4; imageData.data[idx] color[0]; // R imageData.data[idx 1] color[1]; // G imageData.data[idx 2] color[2]; // B imageData.data[idx 3] 200; // A200 让底图透出来 } } ctx.putImageData(imageData, 0, 0); return canvas; } // 简单的蓝-青-黄-红渐变实际项目里可以换成自己的配色 function buildColorLUT() { const lut []; for (let i 0; i 256; i) { const t i / 255; lut.push([ Math.round(255 * Math.min(1, Math.max(0, (t - 0.33) * 3))), Math.round(255 * Math.min(1, Math.max(0, 1 - Math.abs(t - 0.33) * 3))), Math.round(255 * Math.min(1, Math.max(0, (0.66 - t) * 3))) ]); } return lut; }这里的颜色映射用的是线性归一化数据分布偏斜时容易一层色覆盖大面积区域第五节会讲怎么用分位数拉伸替代。网格和 canvas 像素一一对应putImageData效率远高于逐像素fillRect前者是直接写内存后者要跑绘制管线。地图方向在这里是“上北下南”Cesium 贴图时如果方向反了在后面贴图代码里通过缩放负值修正。4.3 贴到地球entity 的 Rectangle 方案与 primitive 的区别canvas 生成之后接下来的问题是“怎么贴到三维地球上”。常见做法是加一个entity用rectangle指定经纬度范围material用ImageMaterialProperty引用 canvas。这个方案声明式、代码少改颜色换图都方便适合插值面这种一张图带整个区域的数据。const krigingCanvas drawGridToCanvas(grid, GRID_W, GRID_H); // 把网格范围投影坐标反投影回经纬度 function unprojectX(x) { return x / (Math.cos(centerLat * Math.PI / 180) * METER_PER_DEGREE); } function unprojectY(y) { return y / METER_PER_DEGREE; } const rectangleEntity viewer.entities.add({ id: krigingSurface, rectangle: { coordinates: Cesium.Rectangle.fromDegrees( unprojectX(xMin), unprojectY(yMin), unprojectX(xMax), unprojectY(yMax) ), material: new Cesium.ImageMaterialProperty({ image: krigingCanvas, transparent: true }), clampToGround: true } });clampToGround: true让矩形贴到地形表面而不是悬在椭球面上如果你的场景加载了真实地形这行是必须的。如果发现贴图方向不对用 Cesium 的坐标系翻转可以解决但更常见的坑在 5.4 节。Cesium 里加载这种面数据entity 和 primitive 的区别值得说清楚。entity 是声明式的适合“我要一个贴地矩形”这种需求代码量少内部帮你处理状态变更primitive 是命令式的需要你手动构造 Geometry、Appearance、材质好处是渲染性能高、可控性强适合大量动态图元。克里金插值面通常只有一张图entity 完全足够如果你要做几十个时刻的插值面轮播且每个都是几百 KB 的 canvasprimitive 配合批量刷新会更稳。另一个区别是 entity 的description、id这些属性天然支持拾取和弹窗这在调试时非常方便。4.4 想要等值线用 turf.isolines 直接抽线插值面色块图之外业务上经常还要等值线比如等降水量线、等浓度线。网格数据已经算好了等值线就不需要再回克里金直接基于网格抽线即可。turf.js 的isolines方法是这里最顺手的工具它内部会基于点集生成三角网再按给定分级值抽取等值线。// 把网格转成 turf 点集每个点带 value 属性 const features []; for (let row 0; row GRID_H; row) { for (let col 0; col GRID_W; col) { features.push(turf.point( [unprojectX(xMin (col / GRID_W) * (xMax - xMin)), unprojectY(yMin (row / GRID_H) * (yMax - yMin))], { value: grid[row][col] } )); } } const pointFC turf.featureCollection(features); // 按值域分 6 级抽等值线 const breaks [5, 10, 15, 20, 25]; const isolines turf.isolines(pointFC, breaks, { zProperty: value }); // 加载到 Cesium每条线画成 polyline isolines.features.forEach(line { const positions line.geometry.coordinates[0].map(coord Cesium.Cartesian3.fromDegrees(coord[0], coord[1]) ); viewer.entities.add({ polyline: { positions: positions, width: 2, material: Cesium.Color.fromCssColorString(#333333).withAlpha(0.8) } }); });turf.isolines的breaks参数是分级边界数组有多少个边界就会抽出多少条线边界值要落在实测数据的值域里才有意义。另一个需要注意的点是zProperty必须和写入点的属性名一致写成value是常见做法。等值线抽出来之后配合 4.2 的色块图一起显示专业感一下就上来了这也是很多环境监测大屏的标准呈现方式。5. cesium克里金插值避坑五个翻车现场的现象、原因与解决5.1 插值面在数据边界外“翘边”现象插值面超出了站点分布范围在凸包外围出现一圈明显的高值或低值“帽子”看起来像面被风吹起来一块。原因普通克里金对预测点的外推能力很弱超出数据覆盖范围后半变异函数主导了权重分配预测值会逐渐回归到数据均值附近。这个回归过程不一定是单调的配合边缘处站点稀疏就会在凸包外形成翘边假象。解决用凸包把插值范围裁掉。常见做法是用 turf.js 的convex生成站点凸包多边形然后对网格点做 point-in-polygon 剔除。更平滑的效果用turf.concave但 concave 的参数敏感性高数据分布不规则时容易生成奇怪的形状。生产项目里用凸包已经足够把网格点过滤掉之后再画 canvas天然就没有翘边区域。5.2 经纬度当平面坐标直接算插值彻底失真现象train的参数直接传经纬度插值面要么是一块块色斑要么整体平滑到看不出梯度少数站点旁边出现不该有的高值孤岛。原因半变异函数内部按距离建模经纬度差到米之间存在非线性关系。北大荒一个站点跨度 10 度距离是“度数”而数值可能是几十两者量纲完全不匹配拟合出来的变程和基台值没有物理意义。解决按 3.2 的投影函数处理把经纬度转成以中心点为基准的米制平面坐标。算完预测值之后渲染 Cesium 时再通过反投影函数把坐标还原成经纬度。这里注意反投影函数要和正投影一一对应我在 4.3 里特意写了unprojectX/unprojectY目的就是保证投影前后能对得上。5.3 网格分辨率一上 150页面直接假死现象数据量和分辨率明明不大但一执行到插值循环页面白屏好几秒风扇狂转控制台还报阻塞警告。原因predict的循环是双重嵌套网格 150×150 是 22500 个预测点每个预测点要和全部站点做矩阵运算。如果站点数再是 300就是 675 万次距离计算和解小矩阵主线程被占满渲染帧率掉到 0。解决两件事。第一网格分辨率控制在 100×100 以内插值面展示的视觉信息量在 80×80 已经足够细节更多的地图操作应该交给图片缩放而不是增加计算点。第二如果确实需要高分辨率把预测计算挪到 Web Worker 里主线程只接收结果。这里有个血泪经验不要在viewer.scene.preUpdate里触发计算每一帧都会重算一次正确做法是初始化时算一次把结果缓存起来。5.4 贴地图片跟地形各算各的时不时报 renderError现象控制台出现viewer.scene.renderError贴地矩形在地形起伏区域闪边、花屏甚至整块消失换个视角又恢复。原因clampToGround的贴地方式是把几何体压到地形高度上对高程数据的分辨率非常敏感。默认的椭球体底图没有高程一切正常一旦加载了真实地形数据矩形边界的四个角按地形高度拉伸中间像素却按纹理 UV 线性映射跟实际地形起伏对不上就会出现撕裂感。另外canvas 尺寸过大也会触发纹理上限浏览器和显卡对单边像素上限普遍是 4096超出后纹理被截断。解决第一canvas 尺寸控制在 1024×1024 以下用的时候让它做双线性缩放视觉损失几乎看不出来。第二项目确实要用真实地形时先用 cesium terrain builder 之类工具把地形切好再决定贴图方式不要在运行时既开地形又挂 clampToGround 的图片面。我排查这类问题时的习惯是先用viewer.scene.globe.depthTestAgainstTerrain false做临时对比看是不是深度测试导致的闪边。5.5 色带断层严重低值区被一种颜色吞掉现象插值面大片区域是同一种颜色只有少数高值点周围有渐变色带明显分层像用了 8 位色一样。原因线性归一化遇到偏斜分布时的问题。监测数据大多服从对数正态或偏态分布低值数据密集、高值数据稀疏线性映射把 90% 的网格像素挤到色带前 10% 的区间里肉眼看起来就是一条色块。解决用分位数拉伸替代线性归一化。先对全部网格预测值排序取 5% 和 95% 分位值作为映射的上下界超出部分截断。这样中间 90% 的数据能铺满整个色带分布细节完全打开。这个技巧几乎适用所有克里金插值的可视化也是前端开发工程师面试题里“数据可视化如何避免颜色失真”的一个标准答案我每次带新人都会强调这一点。6. 从静态出图到动态序列让克里金插值动起来6.1 预渲染多时次图片用 Cesium clock 驱动轮播克里金插值最常见的业务形态不是单一时次而是逐小时、逐日更新。如果你在用户拖动时间轴时才现场算插值性能一定扛不住。我的做法是预渲染把所有时次的插值结果提前算好每个时次生成一张 canvas存成数组运行时只做图片替换。// 假设 timeSeries 是多个时次的预测结果数组 const canvasList timeSeries.map(frame { const c drawGridToCanvas(frame.grid, GRID_W, GRID_H); return c; }); // 用 viewer.clock 驱动轮播 const surfaceEntity viewer.entities.getById(krigingSurface); let currentIndex 0; viewer.clock.onTick.addEventListener(() { const second Math.floor(viewer.clock.currentTime.secondsOfDay); const index second % canvasList.length; if (index ! currentIndex) { currentIndex index; surfaceEntity.rectangle.material.image canvasList[index]; } });注意surfaceEntity.rectangle.material.image直接赋 canvas 对象Cesium 会把它当作纹理更新不需要重新创建 entity。这样切换的粒度是“一张图换另一张图”主线程只有一次纹理上传性能开销远小于重新做克里金预测。计算全部放在初始化阶段把阻塞留在一个可接受的加载时长里交互时保持流畅。6.2 验证插值结果的三个自检习惯插值面画出来之后我一般做三个自检都通过才会说这个面“能用”。第一站点值回验。把站点本身的点画在色块图上层看站点位置的颜色是否和它自身的数值颜色一致。如果偏差超过色带的一到两个等级说明sigma2或alpha参数需要回调。第二留一交叉验证。随机剔除 10% 的站点不参与训练用模型预测这些被剔除的点计算预测值和真实值的均方根误差。这个数值能直接告诉你插值结果的可靠程度比肉眼看图客观得多。第三侧面视角检查。把相机压到接近地面贴着插值面边缘看过去。翘边、闪边、颜色断层在这种视角下都会现原形。这也是我为什么习惯在项目里保留一个viewer.camera.flyTo的调试按钮一键切换俯视和侧视检查效率非常高。这套流程走下来静态单时次和动态序列两个方向都能落地。做这类三维开发最忌讳的就是把克里金当成黑匣子参数乱给、坐标不管、结果不看就往上贴。每次调参前先问自己这一步的数据是什么单位、在哪个坐标系、预期范围是多少。本文还有配套的精品资源点击获取