简介面向前端三维开发者的Cesium克里金插值示例包聚焦在HTML环境中实现空间数据插值与三维可视化适合GIS、气象、地质建模等方向的开发人员学习参考。资源共6个文件包含1个可直接运行的HTML页面、3个JavaScript脚本、1个GeoJSON测试数据文件及1个RAR压缩包文件分工明确整体大小仅57KB轻量且结构清晰。该示例已有1354人学习代码按照数据准备、克里金模型构建、预测点计算、Cesium场景渲染的顺序完整展开演示了基于geostat-js实现插值并用点颜色映射数值强度的过程同时提供参数调整思路来改善插值平滑度。通过阅读和运行本项目可快速掌握在Web端集成空间插值能力的核心步骤并可将GeoJSON替换为自己的数据为环境监测、资源评估等实际业务构建动态3D数据表面提供直接参考。1. 克里金插值遇上 Cesium为什么前端要自己做空间插值手里只有十几个气象站温度读数或者几十个大气监测点的 PM2.5 数据却要在一张三维地球上铺出整个区域的连续分布色面这是气象、环保、地勘行业做前端可视化时绕不开的需求。克里金插值Kriging基于空间自相关统计用少量离散采样点估算任意位置的连续值在浏览器端就能完成整套计算再配合 Cesium 把结果渲染成贴地热力图。整个过程只需要一个 HTML 页面不依赖后端 GIS 服务。这篇笔记写给正在接触 Cesium 三维开发的前端工程师按「算法原理 → 纯 JS 实现 → 三维渲染 → 排错经验」的顺序把最小可运行方案和关键参数讲透。2. 克里金插值原理与选型为什么是它而不是反距离权重2.1 克里金插值的基本思想变异函数与空间自相关克里金插值最早出现在 20 世纪 50 年代的南非金矿储量估算中后来由法国数学家 Georges Matheron 整理成一套完整的地统计学理论。这门学科的基本假设是属性值在空间中不是独立分布的而是存在空间自相关——两个点离得越近它们的属性值越相似。克里金插值就是把这种“越近越像”的直觉量化通过变异函数来描述两点间属性差异随距离变化的规律。变异函数是克里金的核心。对距离为 h 的两个点它们属性值之差的平方的期望记作 γ(h)计算公式是 γ(h) [1/(2N(h))] × Σ[Z(xi) − Z(xih)]²。具体实现时算法会把点对按距离分箱每个箱子里算出一个平均 γ(h)然后用球状、指数或高斯模型拟合出一条平滑曲线。这个拟合结果就是后续所有预测的依据位置越近的采样点获得越高的权重但权重不是单纯按距离平方的倒数而是通过求解线性方程组得到的目标是让预测方差最小。前端开发最容易忽略的一点是克里金不是一把“给你一个插值公式”的锤子而是一套统计推断框架。它不仅输出每个位置的预测值还能输出预测方差。在业务可视化中预测方差场可以直接叠加成一层“可信度热力图”这一点是 IDW 做不到的。理解了变异函数和克里金方差后面调参的时候才不会靠肉眼反复试。2.2 普通克里金的三个关键参数块金值、基台值、变程在普通克里金Ordinary Kriging中变异函数曲线有三个特征参数它们直接决定插值结果的形态。调整这三个参数的优先级排在调整模型类型之前。参数含义对结果的影响块金值nugget距离趋近于 0 时仍存在的方差来自测量误差或小于采样间距的微观变异块金值偏大时插值面显碎细节不稳定容易出现颗粒状色斑基台值sill距离足够大后变异函数趋平的极限值接近样本方差基台值决定色带跨度。基台值太小色带层次会被压缩变程range变异函数到达基台值对应的距离超过变程后空间相关性消失变程小时插值受局部点支配色斑小而多变程大时整体更平滑三种理论模型的选择也直接影响结果。球状模型spherical在变程处平稳到达基台值适合土壤性质、地下水水位这类连续地学数据指数模型exponential渐近到达基台值曲线更长尾适合温度、风场等物理场高斯模型gaussian在原点附近更平缓整体过渡柔和用于高程、降水插值效果通常不错。前端库里一般只需要传入模型名称但要知道选错的后果——插值面不会报错只是层次要么过硬要么过软交叉验证的误差会暴露问题。2.3 与反距离权重、最近邻插值的对比什么时候选克里金做三维可视化时经常有人拿反距离权重IDW和克里金对比。两者的目标一致实现路径不同选型时需要结合数据量、场景和交付标准去看。方法原理优势劣势适合场景最近邻直接取最近点的值零计算成本面状边界不平滑分类数据、地块图IDW距离倒数加权平均算法直观参数少权重幂次敏感易出牛眼快速预览、原型验证克里金变异函数 最小方差估计有误差评估结果平滑参数多计算量大正式分析、交付图层我的选型经验是如果这个图层只是开发中的临时预览用 IDW 足够了计算快代码短如果图层要作为最终成果交付或者需要向别人解释“这个插值结果有多可靠”就用克里金。克里金在数据点相对稀疏、区域特征明显时比 IDW 更能还原真实的空间结构。但要留意克里金不是万能的当采样点少于 15 个时变异函数拟合非常不稳定出现牛眼色斑的几率比 IDW 还高。数据太少时老老实实增加采样点或改用 IDW。3. 纯 JS 在 HTML 页面里跑通克里金插值3.1 引入 kriging.js 并构造采样点数据浏览器端做克里金插值最常见的选择是 kriging.js 这个轻量库无依赖、API 也简单train、predict、grid、plot 四个方法覆盖了从拟合到渲染的完整流程。npm 安装或把源码放进本地 vendor 目录都行在 HTML 里直接用 script 标签引入即可。!DOCTYPE html html langzh-CN head meta charsetUTF-8 meta nameviewport contentwidthdevice-width, initial-scale1.0 title克里金插值最小示例/title script srcvendor/kriging.js/script /head body canvas idkrigingCanvas width800 height600/canvas script // 三个平行数组经度、纬度、观测值 const lngs [116.32, 116.48, 116.41, 116.57, 116.35, 116.52, 116.44]; const lats [39.90, 39.82, 39.95, 39.76, 39.88, 39.93, 39.79]; const values [22.4, 25.1, 24.6, 18.3, 20.8, 26.2, 19.5]; /script /body /html这里把采样点组织成三个平行数组一是 kriging.js 的 train 方法接受的就是这种格式二是方便后续把数据从 GeoJSON 或 CSV 转过来。如果数据源是 GeoJSON 的 FeatureCollection可以把 features 数组 map 成三个数组再传入。注意经纬度顺序很多 GeoJSON 里是 [lng, lat]但有些接口返回 [lat, lng]搞反了训练不会报错但插值结果会像被拧过的布。建议在入库前统一加一个坐标校验的函数。3.2 用 train 拟合变异函数模型与平滑参数怎么设数据准备好之后调用 train 方法拟合变异函数模型返回一个包含模型参数的对象后面所有 predict 都依赖这个对象。// 指数模型sigma20由算法估算alpha10平滑因子 const krigModel kriging.train(values, lngs, lats, exponential, 0, 10); // 训练成功后可以用 predict 在任意点插值 const testValue kriging.predict(116.40, 39.85, krigModel); console.log((116.40, 39.85) 处的预测值:, testValue);train 的六个参数依次是采样值数组 t、x 坐标数组 xs、y 坐标数组 ys、模型名 model、方差初始值 sigma2、平滑因子 alpha。sigma2 通常传 0让算法基于样本自己估计方差alpha 控制变异函数曲线的平滑程度值越大曲线越平缓插值面越柔和。实践时从 alpha10 起步如果结果出现明显的牛眼或条纹把 alpha 加大到 1520 再试。这里有点玄学alpha 和 sigma2 在没有交叉验证的情况下只能靠经验调跑完看一眼结果再决定下一个值。建议在页面上留一个参数面板实时改 alpha 触发重算比改代码刷新浏览器快得多。3.3 用 predict 生成规则格网包围盒和步长的折中单点 predict 没有意义Cesium 渲染需要一张覆盖目标区域的规则格网。做法是确定经纬度包围盒按固定步长逐行逐列预测。// 在采样点范围外扩一个边距避免插值面贴着数据边界 const padding 0.05; const minLng Math.min(...lngs) - padding; const maxLng Math.max(...lngs) padding; const minLat Math.min(...lats) - padding; const maxLat Math.max(...lats) padding; // 格网步长0.004 度约等于 400 米适合城市级范围 const step 0.004; const grid []; for (let lat minLat; lat maxLat; lat step) { const row []; for (let lng minLng; lng maxLng; lng step) { row.push(kriging.predict(lng, lat, krigModel)); } grid.push(row); }双层循环把预测值按“行优先”存成二维数组 grid行对应纬度、列对应经度。step 的选择是个权衡step 太细例如 0.001格网点数量会达到十几万predict 调用 15 万次以上浏览器明显卡顿step 太粗渲染出来是马赛克。我一般按目标区域宽度来定步长宽度不超过 0.5 度时用 0.004跨越 1 度以上时用 0.01 起步。另外predict 内部涉及矩阵求逆样本点数量在几百个级别时开销不小如果页面卡住优先在开发者工具里看主线程占用再决定是放大 step 还是减少训练集。3.4 色彩映射把格网变成 Canvas 图像格网矩阵中的数值必须映射成颜色才能显示。推荐自己实现归一化和色带插值而不是直接调 kriging.plot因为亲手掌控色带后在 Cesium 里做多图层颜色统一会方便很多。// 统计格网中的最小值和最大值 let minV Infinity, maxV -Infinity; for (const row of grid) { for (const v of row) { if (v minV) minV v; if (v maxV) maxV v; } } // 自定义四级色带蓝 - 青 - 橙 - 红 const palette [ [59, 76, 192], [102, 194, 164], [252, 141, 98], [217, 45, 45] ]; // 将数值 t0~1映射为调色板上的 RGBA 颜色 function valueToColor(t) { const seg t * (palette.length - 1); const idx Math.min(Math.floor(seg), palette.length - 2); const frac seg - idx; const c0 palette[idx]; const c1 palette[idx 1]; return [ Math.round(c0[0] (c1[0] - c0[0]) * frac), Math.round(c0[1] (c1[1] - c0[1]) * frac), Math.round(c0[2] (c1[2] - c0[2]) * frac) ]; } // 写入 Canvas每个格网点画成 2x2 像素避免贴图过细 const canvas document.getElementById(krigingCanvas); canvas.width grid[0].length * 2; canvas.height grid.length * 2; const ctx canvas.getContext(2d); const imgData ctx.createImageData(canvas.width, canvas.height); for (let j 0; j grid.length; j) { for (let i 0; i grid[0].length; i) { const t (grid[j][i] - minV) / (maxV - minV); const [r, g, b] valueToColor(Math.max(0, Math.min(1, t))); const base (j * 2 * canvas.width i * 2) * 4; for (let dy 0; dy 2; dy) { for (let dx 0; dx 2; dx) { const p base dy * canvas.width * 4 dx * 4; imgData.data[p] r; imgData.data[p 1] g; imgData.data[p 2] b; imgData.data[p 3] 255; } } } } ctx.putImageData(imgData, 0, 0);valueToColor 接收 01 的归一化值在调色板相邻两个颜色之间做线性插值。格网每个点绘制成 2x2 像素是为了避免 Canvas 尺寸太小导致后续贴地渲染时纹理发糊。生成的 canvas 变量可以直接传给 Cesium 作为贴图源。注意 imgData 的写入顺序是行优先、每行从左到右这个顺序和后面 Rectangle 的 UV 映射是对应关系如果写成列优先热力图会整体转 90 度。4. 把插值结果渲染到 Cesium 三维地球两种落地路线4.1 路线一Canvas 离屏绘制 RectangleGeometry 贴地拿到第三节生成的 canvas 后最直接的渲染方案是把它贴到一个覆盖插值范围的矩形地面上。Entity 的 rectangle 属性加上 ImageMaterialProperty可以在不到二十行代码内完成这项工作。// 先把 Canvas 转成 DataURL供材质使用 const imageUrl canvas.toDataURL(image/png); viewer.entities.add({ id: krigingLayer, rectangle: { coordinates: Cesium.Rectangle.fromDegrees(minLng, minLat, maxLng, maxLat), material: new Cesium.ImageMaterialProperty({ image: imageUrl, transparent: true }) } });这里用 toDataURL 把离屏画布转成材质可用的图片源。需要注意两件事一是 canvas 尺寸过大时 toDataURL 有明显的性能消耗建议把 canvas 宽度控制在 1000 像素以内否则首次加载会出现短暂的空白二是矩形区域四个角是直角而插值区域的边缘往往是圆滑或不规则的。不做处理的话贴图会带着一个明显的矩形边框观感很差。提示toDataURL 在 canvas 较大时会阻塞主线程几百毫秒建议先把 canvas 缩到 800px 宽再转或用离屏 canvas 直接作为 image 传入部分版本的 Cesium。解决办法是在离屏 canvas 上做透明边缘。计算采样点的凸包把凸包外部像素的 alpha 置为 0。如果不想引入 Turf.js可以做一个简化版本对每个像素计算它到最近采样点的距离超过某个阈值就认为该点不受数据支撑设为完全透明。这个简化方案在采样点分布均匀时效果不错但阈值需要试推荐用 Turf.js 的 convex 方法生成凸包再判断。4.2 路线二规则格网 Polygon 逐个叠加保留逐格交互Canvas 贴图的缺点是每个格网不可独立交互。业务方常常会提一个需求“我点击地图上任意位置能不能看到这个位置的插值是多少”用贴图方案你只能反查像素颜色再反算数值误差大且麻烦。如果对交互有要求直接在格网上创建逐个四边形 Entity。for (let j 0; j grid.length - 1; j) { for (let i 0; i grid[0].length - 1; i) { const lng0 minLng i * step; const lat0 minLat j * step; const v grid[j][i]; const t (v - minV) / (maxV - minV); const [r, g, b] valueToColor(Math.max(0, Math.min(1, t))); viewer.entities.add({ rectangle: { coordinates: Cesium.Rectangle.fromDegrees(lng0, lat0, lng0 step, lat0 step), material: Cesium.Color.fromBytes(r, g, b, 200) } }); } }每个四边形只有 0.004 度见方在常见缩放级别下看起来就是一整块连续色面。这个方案的优点很直接每个 cell 都是一个 entity可以挂 description、可以设置 pickable后面要做点击展示数值就很自然。缺点也明显entity 数量等于格网数量3000 个以内流畅超过 5000 个时场景旋转会明显掉帧。大批量数据下改成 Primitive GeometryInstance 批量渲染才压得住但代码复杂度会上升。我一般按最终需求来定路线只是出图展示用路线一要支持点击查看数值、要做联动分析用路线二但把 step 同步调大到 0.0080.01控制 entity 总数。4.3 颜色映射与色带设计从数值到视觉的最后一公里插值结果能不能让业务方信服颜色方案的影响甚至比算法参数更大。常见错误是拉满七彩虹带红橙黄绿青蓝紫全部上场看着唬人实际上把数据的中间层次全压扁了。热力图配色应该遵循“低值→中值→高值”的连续渐变且每段色相变化平缓。温度场推荐蓝→白→红三色渐变污染物浓度推荐蓝→青→黄→红四段地质品位推荐黄→橙→褐。如果数据有极值点线性归一化会让大多数格网的颜色挤在低端需要改用对数拉伸。// 对数拉伸把 0~1 的归一化值映射到对数区间增强低值部分色阶 function logNormalize(t) { if (t 0) return 0; return Math.log10(1 t * 9); // 结果范围约 0~1 }映射时先线性归一化再套 logNormalize再进入 valueToColor配色层次会明显变好。反过来如果数据本身差异不大、要求突出微小变化则避免对数拉伸保持线性。色带和拉伸方式也可以暴露在页面调试面板里让业务方现场切换这比事后改代码高效得多。5. 克里金插值三维可视化避坑指南数据、参数与渲染5.1 采样点太少出现牛眼现象插值结果里某个采样点周围出现一圈圈同心圆状的色斑像靶子一样明显不符合实际分布。原因采样点少于 15 个时变异函数拟合很不稳定权重几乎集中在最近的采样点身上稍远位置很快回到均值形成从高到低的环状过渡。解决优先增加采样点。项目数据无法扩充时把 alpha 平滑因子调大到 1520可以在一定程度上压平牛眼。另外把插值范围适当外扩也能减少边界处的环状效应。这个坑我在第一版示例里踩过七个点跑出的色斑图差点直接拿去给客户演示后来才发现是数据量的问题。5.2 经纬度距离变形导致插值失真现象高纬度地区同一经度差对应的实际距离比低纬度小很多插值结果到了高纬度区域明显沿纬向拉长色带变形。原因克里金内部用欧氏距离计算经纬度是角坐标直接用度数求距离在球面上会产生畸变。解决小范围数据做一次近似投影修正把经度乘以 cos(lat) 转成平面坐标再进入克里金运算。// 训练前对经纬度做等距圆柱投影近似 const factor Math.cos((39.85 * Math.PI) / 180); const xs lngs.map(lng lng * factor); const ys lats; // 用 xs/ys 训练和预测展示时把经度除以 factor 还原这段修正代码对城市级范围1 个经纬度以内足够用省级以上范围要换 UTM 或 Mercator 正规投影。坐标变换的思路虽然简单但在接大范围数据时经常被忽略等到色带变形才发现已经晚了。5.3 插值范围超出采样凸包边缘颜色失真现象插值区域是整个城市的矩形范围采样点只集中在市中心矩形边缘处出现大片不合理的怪异色块颜色向中心快速过渡。原因超出凸包的区域没有采样点支撑克里金的预测方差急剧增大预测值趋于均值视觉上就是一条条向边缘扩散的平淡色带。解决用 Turf.js 或手写射线法生成采样点的凸包多边形只对凸包内部的格网做渲染或者渲染时对凸包外像素做透明处理。先判断每个格网点是否在凸包内再进行下一轮处理即可。经过这一步热力图边缘会贴合数据实际覆盖范围不再是一整块直角矩形。5.4 Canvas 贴图在地表拉伸模糊现象热力图贴到 Cesium 地表面后色带边界糊成一团放大后纹理像隔着一层毛玻璃等高线边界完全看不清。原因离屏 canvas 分辨率偏低矩形几何体拉伸时纹理被浏览器插值平滑。这个问题在地形起伏区域更明显因为地形三角形网格会拉伸纹理。解决Canvas 像素尺寸至少设置为格网宽度的 23 倍写入像素时设置ctx.imageSmoothingEnabled false关闭浏览器对 image 的平滑缩放。另外贴图矩形要与地形贴合使用heightReference: Cesium.HeightReference.CLAMP_TO_GROUND避免漂浮在半空产生视觉错位。5.5 克里金参数过拟合插值面出现条纹现象插值结果出现平行的条纹或波纹状色带采样点越密集的区域越明显看起来不像是自然渐变。原因sigma2 设置过小、alpha 过低模型把采样点噪声当成了真实空间结构导致变异函数曲线过于陡峭权重矩阵病态。解决把 sigma2 提高到样本方差的 10% 以上alpha 加到 1015。用交叉验证算 RMSE在误差和视觉之间取平衡而不是一味追求模型和采样点的精确吻合。6. 进阶把插值做成可交互图层并验证精度6.1 点击查询任意位置的插值预测值Canvas 贴图方案配合事件机制可以在点击位置直接调 predict 算出预测值做到“点击即取值”。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 carto Cesium.Cartographic.fromCartesian(cartesian); const lng Cesium.Math.toDegrees(carto.longitude); const lat Cesium.Math.toDegrees(carto.latitude); const value kriging.predict(lng, lat, krigModel); document.getElementById(pickInfo).textContent 经度 ${lng.toFixed(4)}纬度 ${lat.toFixed(4)}预测值 ${value.toFixed(2)}; }, Cesium.ScreenSpaceEventType.LEFT_CLICK);这里的 pickEllipsoid 在相机视角倾斜时会有少量偏差如果图层贴合地形改用viewer.scene.pickPosition更准确。页面里加一个固定 div 来显示 pickInfo 内容即可。6.2 用交叉验证检查克里金插值精度克里金参数调了半天最终还是用数字说话。把采样点随机分成 80% 训练集和 20% 验证集训练完对验证集预测并计算 RMSE。const indices values.map((_, i) i).sort(() Math.random() - 0.5); const split Math.floor(indices.length * 0.8); const trainIdx indices.slice(0, split); const testIdx indices.slice(split); const trainModel kriging.train( trainIdx.map(i values[i]), trainIdx.map(i lngs[i]), trainIdx.map(i lats[i]), spherical, 0, 10 ); let rmse 0; testIdx.forEach(i { const p kriging.predict(lngs[i], lats[i], trainModel); rmse (p - values[i]) ** 2; }); rmse Math.sqrt(rmse / testIdx.length); console.log(交叉验证 RMSE:, rmse.toFixed(3));RMSE 再除以样本标准差得到相对误差。小于 20% 说明插值结果可用大于 30% 说明要么模型选错要么采样分布有硬伤。每次随机划分结果会有波动多跑几次取平均。我习惯把这个校验逻辑藏在 debug 开关后面项目现场让业务方输入不同 alpha 对比比口头解释“这个色块是误差大”有力气得多。克里金插值做三维可视化最危险的不是算法复杂而是没人校验结果就当真相出图。我现在每接一个项目第一件事就是把交叉验证脚本先写好参数怎么调都有数字兜底。这个习惯帮我避免过好多次交付后被打回改数据的场面。希望帮到你。本文还有配套的精品资源点击获取 纯JS+Cesium：克里金插值热力图三维可视化实战
简介面向前端三维开发者的Cesium克里金插值示例包聚焦在HTML环境中实现空间数据插值与三维可视化适合GIS、气象、地质建模等方向的开发人员学习参考。资源共6个文件包含1个可直接运行的HTML页面、3个JavaScript脚本、1个GeoJSON测试数据文件及1个RAR压缩包文件分工明确整体大小仅57KB轻量且结构清晰。该示例已有1354人学习代码按照数据准备、克里金模型构建、预测点计算、Cesium场景渲染的顺序完整展开演示了基于geostat-js实现插值并用点颜色映射数值强度的过程同时提供参数调整思路来改善插值平滑度。通过阅读和运行本项目可快速掌握在Web端集成空间插值能力的核心步骤并可将GeoJSON替换为自己的数据为环境监测、资源评估等实际业务构建动态3D数据表面提供直接参考。1. 克里金插值遇上 Cesium为什么前端要自己做空间插值手里只有十几个气象站温度读数或者几十个大气监测点的 PM2.5 数据却要在一张三维地球上铺出整个区域的连续分布色面这是气象、环保、地勘行业做前端可视化时绕不开的需求。克里金插值Kriging基于空间自相关统计用少量离散采样点估算任意位置的连续值在浏览器端就能完成整套计算再配合 Cesium 把结果渲染成贴地热力图。整个过程只需要一个 HTML 页面不依赖后端 GIS 服务。这篇笔记写给正在接触 Cesium 三维开发的前端工程师按「算法原理 → 纯 JS 实现 → 三维渲染 → 排错经验」的顺序把最小可运行方案和关键参数讲透。2. 克里金插值原理与选型为什么是它而不是反距离权重2.1 克里金插值的基本思想变异函数与空间自相关克里金插值最早出现在 20 世纪 50 年代的南非金矿储量估算中后来由法国数学家 Georges Matheron 整理成一套完整的地统计学理论。这门学科的基本假设是属性值在空间中不是独立分布的而是存在空间自相关——两个点离得越近它们的属性值越相似。克里金插值就是把这种“越近越像”的直觉量化通过变异函数来描述两点间属性差异随距离变化的规律。变异函数是克里金的核心。对距离为 h 的两个点它们属性值之差的平方的期望记作 γ(h)计算公式是 γ(h) [1/(2N(h))] × Σ[Z(xi) − Z(xih)]²。具体实现时算法会把点对按距离分箱每个箱子里算出一个平均 γ(h)然后用球状、指数或高斯模型拟合出一条平滑曲线。这个拟合结果就是后续所有预测的依据位置越近的采样点获得越高的权重但权重不是单纯按距离平方的倒数而是通过求解线性方程组得到的目标是让预测方差最小。前端开发最容易忽略的一点是克里金不是一把“给你一个插值公式”的锤子而是一套统计推断框架。它不仅输出每个位置的预测值还能输出预测方差。在业务可视化中预测方差场可以直接叠加成一层“可信度热力图”这一点是 IDW 做不到的。理解了变异函数和克里金方差后面调参的时候才不会靠肉眼反复试。2.2 普通克里金的三个关键参数块金值、基台值、变程在普通克里金Ordinary Kriging中变异函数曲线有三个特征参数它们直接决定插值结果的形态。调整这三个参数的优先级排在调整模型类型之前。参数含义对结果的影响块金值nugget距离趋近于 0 时仍存在的方差来自测量误差或小于采样间距的微观变异块金值偏大时插值面显碎细节不稳定容易出现颗粒状色斑基台值sill距离足够大后变异函数趋平的极限值接近样本方差基台值决定色带跨度。基台值太小色带层次会被压缩变程range变异函数到达基台值对应的距离超过变程后空间相关性消失变程小时插值受局部点支配色斑小而多变程大时整体更平滑三种理论模型的选择也直接影响结果。球状模型spherical在变程处平稳到达基台值适合土壤性质、地下水水位这类连续地学数据指数模型exponential渐近到达基台值曲线更长尾适合温度、风场等物理场高斯模型gaussian在原点附近更平缓整体过渡柔和用于高程、降水插值效果通常不错。前端库里一般只需要传入模型名称但要知道选错的后果——插值面不会报错只是层次要么过硬要么过软交叉验证的误差会暴露问题。2.3 与反距离权重、最近邻插值的对比什么时候选克里金做三维可视化时经常有人拿反距离权重IDW和克里金对比。两者的目标一致实现路径不同选型时需要结合数据量、场景和交付标准去看。方法原理优势劣势适合场景最近邻直接取最近点的值零计算成本面状边界不平滑分类数据、地块图IDW距离倒数加权平均算法直观参数少权重幂次敏感易出牛眼快速预览、原型验证克里金变异函数 最小方差估计有误差评估结果平滑参数多计算量大正式分析、交付图层我的选型经验是如果这个图层只是开发中的临时预览用 IDW 足够了计算快代码短如果图层要作为最终成果交付或者需要向别人解释“这个插值结果有多可靠”就用克里金。克里金在数据点相对稀疏、区域特征明显时比 IDW 更能还原真实的空间结构。但要留意克里金不是万能的当采样点少于 15 个时变异函数拟合非常不稳定出现牛眼色斑的几率比 IDW 还高。数据太少时老老实实增加采样点或改用 IDW。3. 纯 JS 在 HTML 页面里跑通克里金插值3.1 引入 kriging.js 并构造采样点数据浏览器端做克里金插值最常见的选择是 kriging.js 这个轻量库无依赖、API 也简单train、predict、grid、plot 四个方法覆盖了从拟合到渲染的完整流程。npm 安装或把源码放进本地 vendor 目录都行在 HTML 里直接用 script 标签引入即可。!DOCTYPE html html langzh-CN head meta charsetUTF-8 meta nameviewport contentwidthdevice-width, initial-scale1.0 title克里金插值最小示例/title script srcvendor/kriging.js/script /head body canvas idkrigingCanvas width800 height600/canvas script // 三个平行数组经度、纬度、观测值 const lngs [116.32, 116.48, 116.41, 116.57, 116.35, 116.52, 116.44]; const lats [39.90, 39.82, 39.95, 39.76, 39.88, 39.93, 39.79]; const values [22.4, 25.1, 24.6, 18.3, 20.8, 26.2, 19.5]; /script /body /html这里把采样点组织成三个平行数组一是 kriging.js 的 train 方法接受的就是这种格式二是方便后续把数据从 GeoJSON 或 CSV 转过来。如果数据源是 GeoJSON 的 FeatureCollection可以把 features 数组 map 成三个数组再传入。注意经纬度顺序很多 GeoJSON 里是 [lng, lat]但有些接口返回 [lat, lng]搞反了训练不会报错但插值结果会像被拧过的布。建议在入库前统一加一个坐标校验的函数。3.2 用 train 拟合变异函数模型与平滑参数怎么设数据准备好之后调用 train 方法拟合变异函数模型返回一个包含模型参数的对象后面所有 predict 都依赖这个对象。// 指数模型sigma20由算法估算alpha10平滑因子 const krigModel kriging.train(values, lngs, lats, exponential, 0, 10); // 训练成功后可以用 predict 在任意点插值 const testValue kriging.predict(116.40, 39.85, krigModel); console.log((116.40, 39.85) 处的预测值:, testValue);train 的六个参数依次是采样值数组 t、x 坐标数组 xs、y 坐标数组 ys、模型名 model、方差初始值 sigma2、平滑因子 alpha。sigma2 通常传 0让算法基于样本自己估计方差alpha 控制变异函数曲线的平滑程度值越大曲线越平缓插值面越柔和。实践时从 alpha10 起步如果结果出现明显的牛眼或条纹把 alpha 加大到 1520 再试。这里有点玄学alpha 和 sigma2 在没有交叉验证的情况下只能靠经验调跑完看一眼结果再决定下一个值。建议在页面上留一个参数面板实时改 alpha 触发重算比改代码刷新浏览器快得多。3.3 用 predict 生成规则格网包围盒和步长的折中单点 predict 没有意义Cesium 渲染需要一张覆盖目标区域的规则格网。做法是确定经纬度包围盒按固定步长逐行逐列预测。// 在采样点范围外扩一个边距避免插值面贴着数据边界 const padding 0.05; const minLng Math.min(...lngs) - padding; const maxLng Math.max(...lngs) padding; const minLat Math.min(...lats) - padding; const maxLat Math.max(...lats) padding; // 格网步长0.004 度约等于 400 米适合城市级范围 const step 0.004; const grid []; for (let lat minLat; lat maxLat; lat step) { const row []; for (let lng minLng; lng maxLng; lng step) { row.push(kriging.predict(lng, lat, krigModel)); } grid.push(row); }双层循环把预测值按“行优先”存成二维数组 grid行对应纬度、列对应经度。step 的选择是个权衡step 太细例如 0.001格网点数量会达到十几万predict 调用 15 万次以上浏览器明显卡顿step 太粗渲染出来是马赛克。我一般按目标区域宽度来定步长宽度不超过 0.5 度时用 0.004跨越 1 度以上时用 0.01 起步。另外predict 内部涉及矩阵求逆样本点数量在几百个级别时开销不小如果页面卡住优先在开发者工具里看主线程占用再决定是放大 step 还是减少训练集。3.4 色彩映射把格网变成 Canvas 图像格网矩阵中的数值必须映射成颜色才能显示。推荐自己实现归一化和色带插值而不是直接调 kriging.plot因为亲手掌控色带后在 Cesium 里做多图层颜色统一会方便很多。// 统计格网中的最小值和最大值 let minV Infinity, maxV -Infinity; for (const row of grid) { for (const v of row) { if (v minV) minV v; if (v maxV) maxV v; } } // 自定义四级色带蓝 - 青 - 橙 - 红 const palette [ [59, 76, 192], [102, 194, 164], [252, 141, 98], [217, 45, 45] ]; // 将数值 t0~1映射为调色板上的 RGBA 颜色 function valueToColor(t) { const seg t * (palette.length - 1); const idx Math.min(Math.floor(seg), palette.length - 2); const frac seg - idx; const c0 palette[idx]; const c1 palette[idx 1]; return [ Math.round(c0[0] (c1[0] - c0[0]) * frac), Math.round(c0[1] (c1[1] - c0[1]) * frac), Math.round(c0[2] (c1[2] - c0[2]) * frac) ]; } // 写入 Canvas每个格网点画成 2x2 像素避免贴图过细 const canvas document.getElementById(krigingCanvas); canvas.width grid[0].length * 2; canvas.height grid.length * 2; const ctx canvas.getContext(2d); const imgData ctx.createImageData(canvas.width, canvas.height); for (let j 0; j grid.length; j) { for (let i 0; i grid[0].length; i) { const t (grid[j][i] - minV) / (maxV - minV); const [r, g, b] valueToColor(Math.max(0, Math.min(1, t))); const base (j * 2 * canvas.width i * 2) * 4; for (let dy 0; dy 2; dy) { for (let dx 0; dx 2; dx) { const p base dy * canvas.width * 4 dx * 4; imgData.data[p] r; imgData.data[p 1] g; imgData.data[p 2] b; imgData.data[p 3] 255; } } } } ctx.putImageData(imgData, 0, 0);valueToColor 接收 01 的归一化值在调色板相邻两个颜色之间做线性插值。格网每个点绘制成 2x2 像素是为了避免 Canvas 尺寸太小导致后续贴地渲染时纹理发糊。生成的 canvas 变量可以直接传给 Cesium 作为贴图源。注意 imgData 的写入顺序是行优先、每行从左到右这个顺序和后面 Rectangle 的 UV 映射是对应关系如果写成列优先热力图会整体转 90 度。4. 把插值结果渲染到 Cesium 三维地球两种落地路线4.1 路线一Canvas 离屏绘制 RectangleGeometry 贴地拿到第三节生成的 canvas 后最直接的渲染方案是把它贴到一个覆盖插值范围的矩形地面上。Entity 的 rectangle 属性加上 ImageMaterialProperty可以在不到二十行代码内完成这项工作。// 先把 Canvas 转成 DataURL供材质使用 const imageUrl canvas.toDataURL(image/png); viewer.entities.add({ id: krigingLayer, rectangle: { coordinates: Cesium.Rectangle.fromDegrees(minLng, minLat, maxLng, maxLat), material: new Cesium.ImageMaterialProperty({ image: imageUrl, transparent: true }) } });这里用 toDataURL 把离屏画布转成材质可用的图片源。需要注意两件事一是 canvas 尺寸过大时 toDataURL 有明显的性能消耗建议把 canvas 宽度控制在 1000 像素以内否则首次加载会出现短暂的空白二是矩形区域四个角是直角而插值区域的边缘往往是圆滑或不规则的。不做处理的话贴图会带着一个明显的矩形边框观感很差。提示toDataURL 在 canvas 较大时会阻塞主线程几百毫秒建议先把 canvas 缩到 800px 宽再转或用离屏 canvas 直接作为 image 传入部分版本的 Cesium。解决办法是在离屏 canvas 上做透明边缘。计算采样点的凸包把凸包外部像素的 alpha 置为 0。如果不想引入 Turf.js可以做一个简化版本对每个像素计算它到最近采样点的距离超过某个阈值就认为该点不受数据支撑设为完全透明。这个简化方案在采样点分布均匀时效果不错但阈值需要试推荐用 Turf.js 的 convex 方法生成凸包再判断。4.2 路线二规则格网 Polygon 逐个叠加保留逐格交互Canvas 贴图的缺点是每个格网不可独立交互。业务方常常会提一个需求“我点击地图上任意位置能不能看到这个位置的插值是多少”用贴图方案你只能反查像素颜色再反算数值误差大且麻烦。如果对交互有要求直接在格网上创建逐个四边形 Entity。for (let j 0; j grid.length - 1; j) { for (let i 0; i grid[0].length - 1; i) { const lng0 minLng i * step; const lat0 minLat j * step; const v grid[j][i]; const t (v - minV) / (maxV - minV); const [r, g, b] valueToColor(Math.max(0, Math.min(1, t))); viewer.entities.add({ rectangle: { coordinates: Cesium.Rectangle.fromDegrees(lng0, lat0, lng0 step, lat0 step), material: Cesium.Color.fromBytes(r, g, b, 200) } }); } }每个四边形只有 0.004 度见方在常见缩放级别下看起来就是一整块连续色面。这个方案的优点很直接每个 cell 都是一个 entity可以挂 description、可以设置 pickable后面要做点击展示数值就很自然。缺点也明显entity 数量等于格网数量3000 个以内流畅超过 5000 个时场景旋转会明显掉帧。大批量数据下改成 Primitive GeometryInstance 批量渲染才压得住但代码复杂度会上升。我一般按最终需求来定路线只是出图展示用路线一要支持点击查看数值、要做联动分析用路线二但把 step 同步调大到 0.0080.01控制 entity 总数。4.3 颜色映射与色带设计从数值到视觉的最后一公里插值结果能不能让业务方信服颜色方案的影响甚至比算法参数更大。常见错误是拉满七彩虹带红橙黄绿青蓝紫全部上场看着唬人实际上把数据的中间层次全压扁了。热力图配色应该遵循“低值→中值→高值”的连续渐变且每段色相变化平缓。温度场推荐蓝→白→红三色渐变污染物浓度推荐蓝→青→黄→红四段地质品位推荐黄→橙→褐。如果数据有极值点线性归一化会让大多数格网的颜色挤在低端需要改用对数拉伸。// 对数拉伸把 0~1 的归一化值映射到对数区间增强低值部分色阶 function logNormalize(t) { if (t 0) return 0; return Math.log10(1 t * 9); // 结果范围约 0~1 }映射时先线性归一化再套 logNormalize再进入 valueToColor配色层次会明显变好。反过来如果数据本身差异不大、要求突出微小变化则避免对数拉伸保持线性。色带和拉伸方式也可以暴露在页面调试面板里让业务方现场切换这比事后改代码高效得多。5. 克里金插值三维可视化避坑指南数据、参数与渲染5.1 采样点太少出现牛眼现象插值结果里某个采样点周围出现一圈圈同心圆状的色斑像靶子一样明显不符合实际分布。原因采样点少于 15 个时变异函数拟合很不稳定权重几乎集中在最近的采样点身上稍远位置很快回到均值形成从高到低的环状过渡。解决优先增加采样点。项目数据无法扩充时把 alpha 平滑因子调大到 1520可以在一定程度上压平牛眼。另外把插值范围适当外扩也能减少边界处的环状效应。这个坑我在第一版示例里踩过七个点跑出的色斑图差点直接拿去给客户演示后来才发现是数据量的问题。5.2 经纬度距离变形导致插值失真现象高纬度地区同一经度差对应的实际距离比低纬度小很多插值结果到了高纬度区域明显沿纬向拉长色带变形。原因克里金内部用欧氏距离计算经纬度是角坐标直接用度数求距离在球面上会产生畸变。解决小范围数据做一次近似投影修正把经度乘以 cos(lat) 转成平面坐标再进入克里金运算。// 训练前对经纬度做等距圆柱投影近似 const factor Math.cos((39.85 * Math.PI) / 180); const xs lngs.map(lng lng * factor); const ys lats; // 用 xs/ys 训练和预测展示时把经度除以 factor 还原这段修正代码对城市级范围1 个经纬度以内足够用省级以上范围要换 UTM 或 Mercator 正规投影。坐标变换的思路虽然简单但在接大范围数据时经常被忽略等到色带变形才发现已经晚了。5.3 插值范围超出采样凸包边缘颜色失真现象插值区域是整个城市的矩形范围采样点只集中在市中心矩形边缘处出现大片不合理的怪异色块颜色向中心快速过渡。原因超出凸包的区域没有采样点支撑克里金的预测方差急剧增大预测值趋于均值视觉上就是一条条向边缘扩散的平淡色带。解决用 Turf.js 或手写射线法生成采样点的凸包多边形只对凸包内部的格网做渲染或者渲染时对凸包外像素做透明处理。先判断每个格网点是否在凸包内再进行下一轮处理即可。经过这一步热力图边缘会贴合数据实际覆盖范围不再是一整块直角矩形。5.4 Canvas 贴图在地表拉伸模糊现象热力图贴到 Cesium 地表面后色带边界糊成一团放大后纹理像隔着一层毛玻璃等高线边界完全看不清。原因离屏 canvas 分辨率偏低矩形几何体拉伸时纹理被浏览器插值平滑。这个问题在地形起伏区域更明显因为地形三角形网格会拉伸纹理。解决Canvas 像素尺寸至少设置为格网宽度的 23 倍写入像素时设置ctx.imageSmoothingEnabled false关闭浏览器对 image 的平滑缩放。另外贴图矩形要与地形贴合使用heightReference: Cesium.HeightReference.CLAMP_TO_GROUND避免漂浮在半空产生视觉错位。5.5 克里金参数过拟合插值面出现条纹现象插值结果出现平行的条纹或波纹状色带采样点越密集的区域越明显看起来不像是自然渐变。原因sigma2 设置过小、alpha 过低模型把采样点噪声当成了真实空间结构导致变异函数曲线过于陡峭权重矩阵病态。解决把 sigma2 提高到样本方差的 10% 以上alpha 加到 1015。用交叉验证算 RMSE在误差和视觉之间取平衡而不是一味追求模型和采样点的精确吻合。6. 进阶把插值做成可交互图层并验证精度6.1 点击查询任意位置的插值预测值Canvas 贴图方案配合事件机制可以在点击位置直接调 predict 算出预测值做到“点击即取值”。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 carto Cesium.Cartographic.fromCartesian(cartesian); const lng Cesium.Math.toDegrees(carto.longitude); const lat Cesium.Math.toDegrees(carto.latitude); const value kriging.predict(lng, lat, krigModel); document.getElementById(pickInfo).textContent 经度 ${lng.toFixed(4)}纬度 ${lat.toFixed(4)}预测值 ${value.toFixed(2)}; }, Cesium.ScreenSpaceEventType.LEFT_CLICK);这里的 pickEllipsoid 在相机视角倾斜时会有少量偏差如果图层贴合地形改用viewer.scene.pickPosition更准确。页面里加一个固定 div 来显示 pickInfo 内容即可。6.2 用交叉验证检查克里金插值精度克里金参数调了半天最终还是用数字说话。把采样点随机分成 80% 训练集和 20% 验证集训练完对验证集预测并计算 RMSE。const indices values.map((_, i) i).sort(() Math.random() - 0.5); const split Math.floor(indices.length * 0.8); const trainIdx indices.slice(0, split); const testIdx indices.slice(split); const trainModel kriging.train( trainIdx.map(i values[i]), trainIdx.map(i lngs[i]), trainIdx.map(i lats[i]), spherical, 0, 10 ); let rmse 0; testIdx.forEach(i { const p kriging.predict(lngs[i], lats[i], trainModel); rmse (p - values[i]) ** 2; }); rmse Math.sqrt(rmse / testIdx.length); console.log(交叉验证 RMSE:, rmse.toFixed(3));RMSE 再除以样本标准差得到相对误差。小于 20% 说明插值结果可用大于 30% 说明要么模型选错要么采样分布有硬伤。每次随机划分结果会有波动多跑几次取平均。我习惯把这个校验逻辑藏在 debug 开关后面项目现场让业务方输入不同 alpha 对比比口头解释“这个色块是误差大”有力气得多。克里金插值做三维可视化最危险的不是算法复杂而是没人校验结果就当真相出图。我现在每接一个项目第一件事就是把交叉验证脚本先写好参数怎么调都有数字兜底。这个习惯帮我避免过好多次交付后被打回改数据的场面。希望帮到你。本文还有配套的精品资源点击获取