☰
Cesium + 克里金插值:把离散监测点变成连续三维场
2026/10/31 16:42:17 网站建设 项目流程

简介:Cesium 克里金插值示例包,面向前端开发与三维可视化学习者,演示如何在 HTML 页面中结合 geostat-js 与 Cesium 实现基于已知点数据的空间插值,并将预测结果渲染为 3D 实体。资源共 6 个文件,包含 3 个 JavaScript 脚本、1 个 HTML 页面、1 个 GeoJSON 数据文件及 1 个压缩包,整体不足 60KB,轻量易读。已有 1354 人学习下载。通过学习可以掌握克里金插值的完整流程:从 JSON/GeoJSON 数据准备、调用 kriging 函数构建插值模型,到依据经纬度与高度预测数值并动态生成 Cesium 点实体;同时可了解不同变异函数模型与最大插值距离等参数对插值精度的影响,适合作为三维 GIS 前端开发的入门参考与实验模板。

1. 克里金插值与Cesium:把离散监测点变成连续三维场的关键一步

在Cesium三维开发中,克里金插值是一个经常被点到名的空间估值方法。做前端开发的人第一次接触它,大多是在气象、污染或地质项目里被客户一句话问住:“我这里有几十个站的监测数据,你能给我在三维地球上画出一张连续的分布面吗?”如果没有克里金插值,你只能在Cesium里画出一堆点,既不成面也看不出趋势。克里金插值的价值,在于它基于空间自相关性,把离散点推算成平滑连续的栅格场,再交给Cesium渲染成三维表面。这份HTML前端开发示例包,正好打通了数据准备、插值计算、Cesium渲染的完整链路,可以直接拿来做原型。适合刚接触Cesium的Web开发,也适合已经被点数据逼疯、想找成熟套路的GIS从业者。

2. 克里金插值的原理与选型:变差函数不是玄学

2.1 为什么是克里金而不是反距离加权

做空间插值,第一反应往往是IDW——距离越近权重越大,代码十分钟就能写出来。但IDW的权重只由距离决定,完全不考虑样本点之间的空间关系,遇到监测点分布不均的时候,会在邻近区域形成明显的“牛眼”效应或扇形假象。克里金则不一样,它先通过半变异函数来刻画“两点之间数值差异随距离变化的统计规律”,再在这个规律上求解一组最优权重,让估值方差最小。换句话说,克里金不仅看“远不远”,更看“空间上像不像”。

半变异函数这个概念第一次接触会有点绕,但你只需要抓住一句话:它描述的是空间自相关性随距离的衰减曲线。横轴是两点之间的距离,纵轴是两点数值差异的统计量。选择的曲线形态不同,后续每个已知点分配的权重就完全不同。在geostat-js里,variogramModel这个参数控制的就是这条曲线的形态,linear、exponential、gaussian、spherical是四种常见形态,分别在平滑性、连续性、局部拟合能力上表现不同。第一次上手时不用深挖每个数学公式,先把“曲线形态→插值面长相”这个对应关系建立起来就够了。

2.2 geostat-js 是第三方算法库,加载顺序会影响生死

Cesium本身不提供克里金实现,开源社区里常用的配套库是geostat-js。它是纯粹的空间统计计算库,和Cesium没有耦合,负责把你给的采样点变成可复用的预测模型。在HTML页面里通过script标签引这三段脚本是最快的跑通方式:

<script src="https://cdn.jsdelivr.net/npm/cesium@1.95.0/Cesium.js"></script> <script src="https://cdn.jsdelivr.net/npm/geostat-js@1.5.3/dist/index.min.js"></script> <script src="https://cdn.jsdelivr.net/npm/geostat-js@1.5.3/dist/kriging.js"></script>

Cesium先加载能保证场景初始化时依赖完整,但这不强制;真正强制的是kriging.js必须在index.min.js之后加载,而且两者都要在你执行插值逻辑之前就位。如果顺序反了,调用geostat.kriging时只会拿到undefined。如果你把这段业务代码放在window.onload里执行,还需要确认脚本的加载时机早于onload触发,否则一样会扑空。

提示:kriging.js依赖index.min.js先加载,两者都要在业务代码执行之前就位。如果加载顺序反了,调用geostat.kriging时只会拿到undefined。

这里补充一个工程化细节。示例包用script标签是为了方便打开即用,但如果你搭的是Vue或React项目,更推荐用npm装geostat-js然后在组件里按需import。全局script方式有个隐藏问题:geostat会挂到window对象上,一旦页面里还有别的库也定义了同名变量,会出现互相覆盖,排查起来很耗时间。另外,script方式把Cesium版本固定在URL里,如果你把1.95.0换成3.x版本,记得同步验证几个关键构造函数的兼容性,我见过有人换了新版Cesium之后,Primitive构造方式变了导致渲染黑屏。

2.3 四个核心参数:variogramModel、sigma2、alpha、maxDistance

拿到示例包后第一件事是把插值配置参数抄出来,搞清楚每个参数在干什么。这是运行前必做功课。下面是geostat-js里最常用的几个配置项,我按实战影响程度排了个序。

参数作用常见取值踩坑提醒
variogramModel变异函数模型linear / exponential / gaussian / sphericalexponential默认最稳,gaussian在小数据集上容易过拟合
sigma2样本方差的估计值默认由数据计算数据波动大时手动设大,否则预测值整体被压缩
alpha显著性水平默认0.1数据量小时可以放宽到0.2,但不建议硬调
maxDistance最大影响距离按站点间距的3~5倍设太小预测点易返回NaN

大多数翻车发生在maxDistance上。站点分布越稀疏,这个值越要放宽,否则模型训练出来之后,未知点一旦离所有样本都远,predict直接吐NaN。反过来,设太大又会让远处的异常点污染局部估值。我的调试套路是先量一下监测点之间的平均距离,乘3到5倍作为maxDistance初值,再用盲测误差迭代修正。

参数调优终归是要在“平滑”和“忠于原始数据”之间找平衡。exponential模型的插值面通常比较顺滑,能掩盖单点噪声;gaussian模型更贴近原始采样值,但对噪声也更敏感,局部会出现锯齿。做环境监测类展示时我一般先上exponential,显得专业且干净。你也可以把同一个数据集用两套参数各跑一遍,对比着看差异,比读论文理解得快。

3. 数据准备与插值计算:从JSON点到克里金模型

3.1 输入数据格式:经度、纬度、值一个都不能乱

geostat-js的kriging函数需要接收一个三维点的数组,每个点就是一组[经度, 纬度, 值]。在WGS84坐标下,经度和纬度使用十进制度。示例包里的数据结构是JSON数组,每个对象包含lon、lat、value三个字段:

[ { "lon": 120.15, "lat": 30.28, "value": 48.5 }, { "lon": 120.18, "lat": 30.31, "value": 62.3 }, { "lon": 120.12, "lat": 30.25, "value": 39.8 } ]

这里的关键是坐标必须统一为WGS84经纬度,不是投影坐标。如果源头数据是墨卡托投影的x/y,正确做法是用proj4转回经纬度,而不是直接当经纬度用。同一个城市的经纬度范围只有零点几度,投影坐标动辄几百万,输反了模型照样能跑,但画到Cesium上会直接飞到海里去,排查时你根本不会往坐标系方向想。

值字段的单位也要统一。现场采集的数据常出现同一字段混着ug/m³、mg/m³不同单位的情况,这在插值计算里不会报错,但会在单位断层处制造出假的峰值。我的习惯是在数据入口加一个标准化函数,把单位换算、缺失值过滤、重复点去重全部做掉,再喂给geostat。去重这一步尤其重要,经纬度完全相同的两个点会让样本矩阵奇异,predict结果变成NaN。

3.2 执行插值:kriging函数与predict函数的分工

插值计算分两步:第一步用kriging函数训练模型,第二步用模型.predict预测未知点。示例包里建议把模型定义为全局变量,页面初始化时训练一次,后续预测点时循环复用。下面是一段标准的调用方式:

// 从数据接口或本地JSON获取监测点 const data = await fetch('./data/station.json').then(res => res.json()); // 把对象数组映射成geostat-js需要的二维数组 const points = data.map(item => [item.lon, item.lat, item.value]); // valueIndex传2,表示数组索引2是值字段 const krigingModel = geostat.kriging(points, 2, { variogramModel: 'exponential', maxDistance: 50000, sigma2: 10 }); // 预测任意经纬度处的插值结果,注意传参顺序是经度、纬度 const prediction = krigingModel.predict([120.16, 30.29]); console.log('该点插值结果:', prediction.toFixed(2));

上面这几段代码是示例里最核心的计算流程。valueIndex这个参数容易被忽略,它告诉geostat-js二维数组中第几列是需要插值的值,这个位置不仅决定值字段,还决定了后续变异函数拟合的因变量。从接口拿数据时如果字段排序变了,这里要同步改,否则会拿经纬度当值去拟合,最终结果完全不可用。

predict参数顺序是[经度, 纬度],不是[纬度, 经度]。很多做GIS出身的人习惯lat/lon顺序,在Cesium的API里也常见lat在前,但这个库的predict明确是lon在前。顺序调反不会报错,所有结果点会在场景里对称偏移,属于最难排查的隐形错误。如果你发现预测值看起来合理但点位全部跑到另一个方向,先去检查这里的参数顺序。

3.3 模型验证:训练点回判与盲测点

插值模型建好了,先别急着上Cesium,用盲测做一次质量检查。把原始站点按五份拆分,四份训练一份验证,误差统计通过之后再去做可视化。这是避免“看着好看但全是噪声”的关键一步。

// 五分之一抽样做盲测 const trainSet = data.filter((_, index) => index % 5 !== 0); const testSet = data.filter((_, index) => index % 5 === 0); const model = geostat.kriging( trainSet.map(d => [d.lon, d.lat, d.value]), 2, { variogramModel: 'gaussian', maxDistance: 30000 } ); testSet.forEach(d => { const pred = model.predict([d.lon, d.lat]); console.log(`真实值 ${d.value}, 预测值 ${pred.toFixed(2)}`); });

留一法交叉验证能同时暴露两类问题:一是模型整体误差偏大,说明变异函数选型不当或数据本身缺乏空间相关性;二是局部点位误差突出,说明该点周围的数据形态和整体形态不一致。遇到第二种情况,不要删点,先检查原数据是否异常,或者考虑对该区域单独建立局部插值模型。抽样时如果数据本身是按地理区域排序的,直接按顺序每五个取一个会让验证集全部落在同一片区域,稳妥的做法是先把数组shuffle再切分。

4. 把插值结果画进Cesium:从单点到整面

4.1 先用Entity验证单点渲染链路

模型验证通过后,先不要一次性渲染整个插值面,而是先做单点验证。确认Cesium Viewer能正常新建,Cesium.Cartesian3.fromDegrees能把经纬度坐标正确转换成场景坐标,然后按值域映射一个渐变色:

const viewer = new Cesium.Viewer('cesiumContainer', { baseLayerPicker: false, animation: false, timeline: false }); // 预测一个具体经纬度的值 const q = krigingModel.predict([120.16, 30.29]); // 全局值域归一化,映射到蓝到红的色带 const ratio = Cesium.Math.clamp((q - globalMin) / (globalMax - globalMin), 0, 1); const color = Cesium.Color.fromHsl(0.7 - ratio * 0.7, 1.0, 0.5, 0.8); viewer.entities.add({ position: Cesium.Cartesian3.fromDegrees(120.16, 30.29), point: { color: new Cesium.ConstantProperty(color), pixelSize: 10 } });

这段代码里有两个细节要注意。globalMin和globalMax必须在预测之前由全部原始站点统计出来,并作为全局常量使用。如果定义在循环里用局部最小值做归一化,每个点都会把自己的值映射成红色,最终画面会变成整片红,看不出任何梯度。另外,Entity的position属性可以直接传位置对象,但明确写出ConstantProperty更利于后续动态更新。

单点测试的意义在于快速验证三件事:坐标转换、颜色映射逻辑、渲染层的可见性。任何一环有问题,在这个最小例子里几行代码就能定位,比直接刷整个网格再回头查要快得多。这一个步骤我省过,结果网格一起渲染的时候发现所有点全都叠在地心,排查了半天才发现是Cesium.Cartesian3.fromDegrees的调用参数顺序出了问题。

4.2 批量预测:用网格把离散点推成连续表面

插值的本质是逐点预测。把目标区域按经纬度步长切成规则网格,每个网格点调一次predict,就能得到一组连续变化的值。再用面片渲染出来,就是视觉上的连续表面。步长的选择直接决定视觉效果和计算成本。

// 定义研究区域范围与步长 const west = 120.10, east = 120.22; const south = 30.22, north = 30.36; const step = 0.002; // 约200米一个点 const grid = []; for (let lon = west; lon <= east; lon += step) { for (let lat = south; lat <= north; lat += step) { const value = krigingModel.predict([lon, lat]); grid.push({ lon, lat, value }); } } console.log('网格点数:', grid.length);

step参数每缩小一半,网格点数变成原来的四倍。当站点数量多、数据范围大时,predict耗时随之增加,页面可能卡好几秒。示例包里的数据规模不大,可以直接跑;在真实项目里建议把这段逻辑放进Web Worker,主线程只接收一次性计算完的网格数组,避免阻塞Cesium的渲染线程。这个优化属于老生常谈,但真的能救回一倍的帧率。

生成的grid数组就是后面渲染的原材料。每个网格点都有经纬度和插值结果,接下来要做的就是根据值域给它上色,然后以点或面的形式加入场景。

4.3 渲染方案:Entity点云与Primitive批量实例化

网格点从几千到几万,渲染方式的选择会直接影响交互流畅度。Entity虽然定位直观、调试方便,但每个Entity在Cesium里都是一个独立对象,几千个Entity一起出现,场景管理开销非常明显。Cesium Primitive则是批量实例化绘制,同样的数据量帧率要稳得多:

// Primitive实例化点云渲染 const instances = grid.map(p => { const ratio = (p.value - globalMin) / (globalMax - globalMin); const color = Cesium.Color.fromHsl(0.7 - ratio * 0.7, 1.0, 0.5, 0.8); return new Cesium.GeometryInstance({ geometry: new Cesium.PointGeometry({ position: Cesium.Cartesian3.fromDegrees(p.lon, p.lat, 0), vertexFormat: Cesium.VertexFormat.POSITION_ONLY }), attributes: { color: Cesium.ColorGeometryInstanceAttribute.fromColor(color) } }); }); viewer.scene.primitives.add(new Cesium.Primitive({ geometryInstances: instances }));

从Entity迁移到Primitive的收益在千点以上最为明显:图元数量不再影响场景树遍历,浏览器端只需要维护一份buffer。代价是实例管理的代码变多了,每个点都要显式创建GeometryInstance。示例包里两种方案都有实现,我建议先跑Entity版本验证算法,再切Primitive版本解决性能,这样可以隔离调试变量。

如果希望插值面贴合地形,还可以在生成网格点坐标时传入高度。Cesium提供terrain采样接口,常见做法是先生成经纬度网格,再用地形服务取高程,这样插值面虽然代表的是value分布,但表面高低起伏呈现的是真实地形轮廓,视觉上更像一张贴在地球表面的数据膜。需要注意的是,viewer没有启用地形时,任何高度值都会被当作零高或忽略,接口调用也不会报错,只会看到平面,这块容易产生“我设了高度怎么没反应”的困惑。

4.4 网格面的上色策略:色带映射而不是逐点随机

有些新手做出来的插值面花花绿绿每一块都不连通,问题往往出在色带策略上。插值面的价值在于让人一眼看出来“哪里高、哪里低、哪里过渡”,这要求整张图使用同一套颜色映射函数。用HSL沿着色相轴从蓝渐变到红,是最直观的方案:

function getValueColor(value, min, max) { const ratio = Cesium.Math.clamp((value - min) / (max - min), 0, 1); const hue = (1 - ratio) * 0.7; // 0.7=蓝, 0=红 return Cesium.Color.fromHsl(hue, 1.0, 0.5, 0.8); }

这个函数建议单独放到公共工具里,所有渲染路径都调用同一个getValueColor,保证Entity版本、Primitive版本、后期换成Polygon版本时颜色完全一致。min和max不要用数组的Math.min/Math.max同步算,用一次全量遍历的结果缓存下来即可。在高值区想突出重点时,可以把饱和度在ratio大于0.8的区间做二次插值,让亮点更干净,但不建议把对比度拉太高,否则插值面噪点感会加重。

5. 避坑与排查:Cesium + 克里金常见问题记录

5.1 现象:插值结果整体漂移,站点落到海里

原因:数据源坐标实际是投影坐标或GCJ-02坐标,被当作WGS84经纬度直接喂给了模型。投影坐标数值范围大,偏差会被Cesium放大到几十公里;GCJ-02的偏移则更隐蔽,看起来量级正常但整体平移了几百米。

解决:先确认数据坐标系。如果是投影,用proj4转回WGS84后再进模型;如果是GCJ-02,先做坐标纠偏。判断方法很简单:在Cesium里随机选一个原始站点经纬度,看看落点和真实监测站是否重合,偏差超过几十米就要先解决坐标系问题再谈插值。

5.2 现象:调用geostat.kriging报undefined,或window上找不到geostat

原因:加载顺序不对。kriging.js在index.min.js之前执行,或者执行业务代码时外部脚本还没同步加载完成。有些情况下,业务代码写在DOMContentLoaded回调里,但浏览器对多个外部脚本的解析顺序和你理解的script排列顺序不一致,也会出现这种问题。

解决:把两个script标签放在业务代码之前,保持index.min.js在前、kriging.js在后。如果是异步加载,不要在DOMContentLoaded里直接调用,改成检查window.geostat存在后再执行。用npm方式引入一般不会遇到这个问题,反而要留意打包时geostat-js的ESM入口被tree-shaking意外剪掉。我建议在调用前加一行console.log(window.geostat),确认对象存在再继续。

5.3 现象:网格点一多页面卡死,CPU占用率飙升

原因:predict是纯CPU计算,双重循环里做大量predict把主线程占满,渲染线程无法及时绘制。Cesium的绘制本身压力不大,瓶颈在插值计算。

解决:把批量预测放进Web Worker,完成后把结果数组postMessage回主线程再渲染;或者在循环里做分帧调度。最直接的优化是把step调大降网格密度,先跑通再谈精细度。经验数据是两千个网格点以上建议走Worker,一万个以上必须走,否则用户看到的就是一张白屏加一个转圈的鼠标指针。

5.4 现象:插值面出现明显色带断层,过渡突兀

原因:网格数据里混入了NaN或超出值域范围的异常点;也可能是每个点单独用局部min/max做归一化,导致同一数值在不同区域渲染出不同颜色,画面看起来像打了补丁。

解决:全局统一次min和max,所有点共用一份映射参数。对NaN做显式检测——predict异常时有时会返回-9999之类的占位值,处理时过滤掉或置为null。在getValueColor入口加一个Number.isFinite判断,兜住所有异常值,色带过渡会自然得多。

5.5 现象:predict结果全是NaN,插值面一片空白

原因:预测点距离所有已知点都超过maxDistance,样本矩阵在求解时没有有效邻居;或者训练数据中存在经纬度重复的点,导致样本矩阵奇异。数据量过少(比如只有三四个站点)时,变异函数模型拟合本身就会失败。

解决:把maxDistance调大覆盖整个研究区域。对输入点先做去重,经纬度相同的数据取平均值或保留最新值。站点数太少时克里金没有意义,至少保留15个以上的点再谈插值。如果必须用少量点出图,退而求其次用IDW反而更稳。

这几条是我拿这个示例包跑业务数据时真实遇到过的,尤其是坐标系漂移和NaN这两条,排查起来非常费时间。建议在接入业务数据之前,把上述检查做成工具函数放在数据管道的最前端,能省掉后续大量debug时间。

6. 进阶:把插值面做成色带热力图并贴到地球上验证

6.1 在Cesium里构建色带热力图

前面讲的是点云渲染,真实交付时更常用的是热力图形态。思路是先把网格点之间的空白填起来——用Polygon把每个网格单元画出来,颜色用统一的色带映射,拼接起来就是一张贴合地球表面的热力图。做这件事的关键是构造几何体时把每个网格单元的四角经纬度都算出来,按规则四边形的顶点顺序提交:

grid.forEach((p, i) => { // 以格点为中心向外偏移半个步长做面片 const lon1 = p.lon - step / 2, lon2 = p.lon + step / 2; const lat1 = p.lat - step / 2, lat2 = p.lat + step / 2; viewer.entities.add({ polygon: { hierarchy: Cesium.Cartesian3.fromDegreesArray( [lon1, lat1, lon2, lat1, lon2, lat2, lon1, lat2] ), material: new Cesium.ColorMaterialProperty( getValueColor(p.value, globalMin, globalMax) ), clampToGround: true } }); });

clampToGround设为true后,多边形会贴到当前地形表面,不会悬在半空,这是三维地球展示插值面最常用的呈现方式。如果大量Polygon影响性能,Cesium也提供GroundPrimitive来合并地面图元,本质上和Primitive思路一致,只是针对贴地多边形做了优化。

这里有个容易被忽略的交付细节:Polygon的hierarchy顶点顺序必须沿顺时针或逆时针闭合,否则部分版本会出现面片扭曲或渲染反面。建议把经纬度数组按左上、右上、右下、左下依次填写,简单且不出错。

6.2 验证方法:站点叠图与十字盲测

交付前我固定执行一套验证流程。第一步,把原始的监测站点用白色点画在地图上,和插值面对照,原站点处插值面的颜色应该和它的实测值量级基本对应;第二步,随机抽几个点位做盲测,在Cesium里用明显颜色的标记点标注,同时console输出预测值和实测值的误差;第三步,切换地形服务,看clampToGround后的热力图是否依然贴合。

// 十字盲测点位定义 const testPoints = [ { lon: 120.16, lat: 30.28, actual: 51.2 }, { lon: 120.17, lat: 30.30, actual: 60.5 } ]; testPoints.forEach(t => { console.log(`实测${t.actual}, 插值${krigingModel.predict([t.lon, t.lat]).toFixed(2)}`); });

这一步做下来,数值是否可靠、贴地是否准确一目了然,比任何视觉效果都更有说服力。做三维可视化项目最怕客户问“你这颜色代表什么值,可信吗”,有了盲测误差列表,至少能给出一份量化边界,哪些区域可信、哪些区域是外推都说得清楚。

从那以后我每次接Cesium插值项目,都强制自己按这套流程走:数据清洗 → 坐标系确认 → 模型盲测 → 单点渲染 → 批量网格 → Primitive优化 → 贴地交付。每一步留下可复现的日志和截图,任何环节出现异常都能快速定位到是数据问题还是渲染问题。这套流程用久了就会发现,克里金在自己的业务场景里其实就是“空间相关性合理刻画”的统计工具,难点从来不在算法本身,而在数据与场景的衔接。这份示例包把最关键的几步已经拆好了,希望这篇拆解能帮你在自己的项目里少走一段弯路。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询