Cesium克里金插值：前端三维温度场可视化渲染

简介：这是一份基于Cesium库与克里金插值算法实现温度数据三维渲染的工程资源，适合WebGIS开发者、三维可视化学习者以及从事气象、海洋等空间数据分析的技术人员。资源以离散温度测站数据为起点，演示了变异函数分析、半方差拟合、未知点位估值等克里金插值核心流程，并将结果集成到Cesium三维地球中进行连续温度场渲染，通过颜色渐变直观呈现温度高低。压缩包共计五百八十五个文件，以三百二十一个JavaScript脚本为核心，一百四十一个PNG与五十三个JPG图片作为图标纹理，二十九个JSON文件保存插值参数与场景配置，二十七个CSS和十个SVG负责界面样式，另有少量wasm等格式辅助，整体仅七点二七MB。借助该包，读者可以快速梳理Cesium场景搭建、克里金算法调用、色带映射与图层叠加的实现思路，理解从离散站点到连续曲面的数据管线，并直接移植到其他空间插值可视化项目中。该资源目前已有一百六十七人学习，对于想深入了解三维气象渲染、温度场构建以及Cesium与克里金结合的开发者，具有直接参考价值。

1. 一个能在地球上画温度场的项目：Cesium 与克里金插值怎么配合

做气象可视化时经常遇到一种尴尬：手里只有几十个离散监测站的数据，但地图上想看到一片连续的温度场。之前我的做法是把数据丢给后端跑 GIS，等出图、再切成瓦片，一来一回少说半天。某次临时要做一个交互 Demo，被迫把克里金插值直接搬到前端，结果发现用 JavaScript 实现普通克里金并不玄学，配合 Cesium 的多边形格网，几分钟就能在三维地球上渲染出一张完整的温度场。这份资源就是一个完整的 Cesium + Kriging 插值温度渲染示例，适合做气象大屏、环境监测、农业区划的朋友直接改。它把「插值计算」和「三维渲染」拆成两个独立模块，你可以单独替换算法，也可以单独换渲染方案。

2. 项目结构：从站点CSV到球面温度场的数据流水线

2.1 资源包目录：插值、渲染、配置三权分立

拿到资源包后，首先看目录结构。整个项目不是把所有逻辑塞进一个 HTML，而是按职责拆成三层：插值层、渲染层、配置层。这样做的直接好处是，你不需要为了换一个插值算法去翻渲染代码，也不需要为了改颜色方案去动数据读取逻辑。

kriging-cesium/ ├── index.html # 入口页面：初始化 Cesium Viewer ├── config.js # 全局配置：采样点路径、格网分辨率、颜色方案 ├── src/ │ ├── kriging.js # 克里金插值核心：变异函数拟合、矩阵求解、格网生成 │ ├── colorRamp.js # 温度值 -> 颜色的映射，内置蓝青黄红渐变 │ └── cesiumRenderer.js # Cesium 场景构建、格网坐标转换、Primitive 合并 ├── data/ │ ├── stations.csv # 模拟站点数据：id,lng,lat,temp │ └── grid.json # 可选的预生成格网缓存，避免每次刷新都重算 └── libs/ └── cesium/ # 推荐放本地 Cesium 静态包，避免 CDN 加载失败

入口文件只负责组装：读取 config → 加载 stations.csv → 调用 kriging.js 生成格网 → 调用 colorRamp.js 得到颜色 → 交给 cesiumRenderer.js 渲染。每一层的输入输出都是普通对象，没有全局变量污染。config.js 里我一般会留几个关键开关，比如gridCacheEnabled，改成true就直接读取 data/grid.json，做演示的时候刷新页面零等待。

2.2 数据格式约定：站点CSV与格网JSON

插值计算的第一步是读取站点数据。这里固定用 CSV 作为输入，因为气象站点的导出工具基本都支持 CSV，手工编辑也方便。字段顺序很重要，项目里统一为id,lng,lat,temp，如果有额外属性不影响解析，但前四列的顺序别乱改。

id,lng,lat,temp S001,112.52,27.81,21.3 S002,113.03,28.19,22.8 S003,111.97,28.55,19.6 S004,113.41,27.31,20.4 S005,114.05,28.02,23.1

执行插值之后，会生成一个规则格网对象。为了排查问题，我建议在调试时把格网存成 JSON 看一眼结构，确认values长度等于nx * ny，不然渲染阶段经常会因为数组越界而报错。

{ "nx": 81, "ny": 61, "xmin": 111.5, "xmax": 114.2, "ymin": 27.1, "ymax": 29.0, "values": [21.3, 21.7, 22.1, 22.5, 22.9, ...] }

这里nx和ny是每个方向的插值点数，values按行优先存储，也就是先固定经度方向，遍历纬度方向。这个约定在后面渲染时决定坐标索引怎么算，写反了会出现温度场整体错位甚至跑到地球背面的情况。数据格式这一层看起来简单，却是整个项目最容易翻车的地方。

2.3 渲染策略对比：多边形格网、Canvas纹理与Entity Wall

在 Cesium 里画连续温度场，常见思路有三种：直接铺多边形小格网、用 Canvas 纹理贴到椭球体、用 Entity 的 Wall 拉伸。我最终选了多边形格网，而不是看起来更「高级」的纹理方案。

纹理方案最大的坑在于 Cesium 把地表纹理映射到球面时，会带一层经纬度到 UV 的间接转换，温度场的边缘容易拉伸变形，而且一旦要做鼠标悬浮查看某个格网的温度，纹理方案很难命中。多边形格网用经纬度直接建 Polygon，虽然 DrawCall 压力大，但通过 PerInstanceColor 合并后性能在几百个站点规模下完全能接受。如果遇到大范围高分辨率场景，后面章节会讲怎么用分片渲染补救。

3. 克里金插值实现：从空间自相关到普通克里金方程组

3.1 变异函数拟合：球状模型参数估算

克里金插值的核心不是「距离越近权重越大」，而是用变异函数描述距离与空间自相关之间的关系。生产中常用球状模型：如果两点距离小于变程a，半方差随距离增加；超过a后，半方差稳定在一个基台值附近。项目里我实现了一个最简版的参数拟合，适合温度这种空间连续变化的场。

// 球状模型半方差：h为距离，a为变程，c0为块金值，c1为基台值 function spherical(h, a, c0, c1) { if (h === 0) return 0; if (h >= a) return c0 + c1; const ratio = h / a; return c0 + c1 * (1.5 * ratio - 0.5 * Math.pow(ratio, 3)); } // 用样本点两两距离的经验值去估算模型参数 function fitVariogram(samples) { const maxDist = Math.max(...samples.map(s => s.maxDist)); const variance = samples.reduce((s, v) => s + v, 0) / samples.length; return { a: maxDist / 3, c0: 0.05 * variance, c1: variance }; }

以上代码是简化写法，更严谨的做法是统计距离区间、拟合曲线，但温度场在中小尺度空间里用固定比例也能得到不错效果。参数上重点看三个值：a控制影响半径，a太大会让远距离站点也参与插值，导致温度场过度平滑；c0是块金值，值越大，邻近样本点的权重会被拉低，能抑制“牛眼”现象；c1通常取样本方差。如果你发现插值结果出现一个个孤立的圆形色斑，优先调大c0，这是我反复踩出来的经验。

3.2 普通克里金方程组：高斯消元求权重

普通克里金除了要求预测值无偏，还要求估计方差最小，最终落到一个带约束条件的线性方程组。目标点周围的每个样本点都有一个权重，权重由样本点之间的半方差矩阵和目标点与样本点之间的半方差向量共同决定。这里我直接写了一个简单的高斯消元求解器，没有引入第三方矩阵库。

function solveLinear(A, b) { const n = b.length; for (let col = 0; col < n; col++) { let pivot = -1; for (let row = col; row < n; row++) { if (Math.abs(A[row][col]) > 1e-10) { pivot = row; break; } } if (pivot === -1) throw new Error('矩阵奇异'); [A[col], A[pivot]] = [A[pivot], A[col]]; [b[col], b[pivot]] = [b[pivot], b[col]]; for (let row = 0; row < n; row++) { if (row === col) continue; const factor = A[row][col] / A[col][col]; for (let j = col; j < n; j++) A[row][j] -= factor * A[col][j]; b[row] -= factor * b[col]; } } return b.map((x, i) => x / A[i][i]); } function ordinaryKriging(target, samples, variogram) { const n = samples.length; const A = []; const b = new Array(n + 1).fill(0); for (let i = 0; i < n; i++) { A[i] = new Array(n + 1).fill(0); for (let j = 0; j < n; j++) { A[i][j] = variogram(samples[i].dist, variogram); } A[i][n] = 1; A[i][i] += 1e-8; b[i] = variogram(target.dist, variogram); } A[n] = new Array(n + 1).fill(1); b[n] = 1; const weights = solveLinear(A, b).slice(0, n); return samples.reduce((sum, s, i) => sum + weights[i] * s.value, 0); }

代码里A[i][i] += 1e-8是为了防止样本点几乎重合时矩阵奇异，这一点能用高斯消元顺利解出来的关键。普通克里金通过拉格朗日乘子把约束条件融入矩阵，所以最后一行和最后一列都是 1，求解得到的最后一维就是拉格朗日乘子，没有实际用途但必须参与求解。如果你的站点数量超过 300，建议先做空间索引，只取目标点周围最近的 20 个样本点参与求解，否则矩阵规模会急剧膨胀。

3.3 生成规则格网：遍历经纬度写入值

有了单点插值函数，生成格网就变得很直接。你要先确定覆盖范围，这块范围最好略大于站点分布的外包矩形，避免边缘温度场缺少数据兜底。然后按分辨率拆分成长方形网格，逐格计算温度值。

function generateGrid(stations, options) { const { nx, ny, xmin, xmax, ymin, ymax } = options; const dx = (xmax - xmin) / nx; const dy = (ymax - ymin) / ny; const values = []; for (let i = 0; i <= nx; i++) { for (let j = 0; j <= ny; j++) { const x = xmin + i * dx; const y = ymin + j * dy; const v = ordinaryKriging( { x, y, dist: 0 }, stations.map(s => ({ ...s, dist: haversine(x, y, s.lng, s.lat) })), variogram ); values.push(v); } } return { nx: nx + 1, ny: ny + 1, xmin, xmax, ymin, ymax, dx, dy, values }; }

这个循环是 O(nx * ny * m²) 的复杂度，m是站点数量。如果nx=80, ny=60, m=50，意味着大约 40 万次距离计算，浏览器里耗时约 2 到 4 秒。想提速可以先把站点两两距离的半方差缓存起来，每次调用时直接查表，而不是重新算距离。这个优化在调试时能明显减少等待感，跑完之后把格网存成 JSON，后面几章渲染阶段就不需要重复计算了。

4. Cesium 渲染：把插值结果变成球面上的色块

4.1 坐标转换与格网构建

经纬度格网和 Cesium 球面的坐标体系不一样，Cesium 内部用三维笛卡尔坐标，而且角度必须转成弧度。在渲染之前，要把每个小格网的四个角从经纬度转换成Cartesian3，并按逆时针顺序排列成 PolygonHierarchy，否则出现多边形重叠或者渲染方向错误。

function gridToCartesianCells(grid) { const cells = []; for (let i = 0; i < grid.nx - 1; i++) { for (let j = 0; j < grid.ny - 1; j++) { const west = Cesium.Math.toRadians(grid.xmin + i * grid.dx); const south = Cesium.Math.toRadians(grid.ymin + j * grid.dy); const east = Cesium.Math.toRadians(grid.xmin + (i + 1) * grid.dx); const north = Cesium.Math.toRadians(grid.ymin + (j + 1) * grid.dy); const height = 0; cells.push({ positions: [ Cesium.Cartesian3.fromRadians(west, south, height), Cesium.Cartesian3.fromRadians(east, south, height), Cesium.Cartesian3.fromRadians(east, north, height), Cesium.Cartesian3.fromRadians(west, north, height) ], value: grid.values[i * grid.ny + j] }); } } return cells; }

这里有个容易踩坑的地方：grid.values的索引是i * grid.ny + j，但我在generateGrid里返回的nx已经是加一之后的值，渲染时循环用的是grid.nx - 1，所以索引必须和生成时保持一致。如果你在生成和渲染环节使用了不同的循环顺序，温度场就会变成一条条斜线。高度设置成 0 是指贴在地表，如果你的温度场需要放在某个高度上，可以统一调整height。

4.2 颜色映射与实例属性

温度值本身是浮点数，要渲染成色块就得做归一化和颜色插值。项目里预设了一套适合温度展示的色带：低温蓝色、中温青色、高温黄色、极高温红色，这套颜色对应到colorRamp.js里。

function createColorRamp() { const stops = [ { t: 0.0, color: [0, 0, 255] }, { t: 0.5, color: [0, 255, 255] }, { t: 0.75, color: [255, 255, 0] }, { t: 1.0, color: [255, 0, 0] } ]; return function (value, min, max) { const t = Math.max(0, Math.min(1, (value - min) / (max - min))); let lower = stops[0], upper = stops[stops.length - 1]; for (let i = 0; i < stops.length - 1; i++) { if (t >= stops[i].t && t <= stops[i + 1].t) { lower = stops[i]; upper = stops[i + 1]; break; } } const local = (t - lower.t) / (upper.t - lower.t); const r = Math.round(lower.color[0] + local * (upper.color[0] - lower.color[0])); const g = Math.round(lower.color[1] + local * (upper.color[1] - lower.color[1])); const b = Math.round(lower.color[2] + local * (upper.color[2] - lower.color[2])); return `rgb(${r},${g},${b})`; }; }

min和max一般取数据里的实际最小/最大温度，而不是固定 -20 和 40。如果你发现渲染出来的色块几乎全是一个颜色，十有八九是归一化分母写死了，导致温度区间被压得极窄。颜色插值里还有一个细节：stops的t不一定要等分，比如高温区想更敏感，可以把 0.8 处加一个颜色节点，色带会自动在区间内插值。

4.3 性能优化：用 GeometryInstance 合并 DrawCall

直接用viewer.entities.add添加几千个 Polygon 会卡到没法看。更好的做法是使用Cesium.Primitive+GeometryInstance，把每个小格网作为独立的 Instance 塞进同一个 Primitive，颜色用逐实例属性传给 GPU。

const instances = cells.map(cell => { return new Cesium.GeometryInstance({ geometry: new Cesium.PolygonGeometry({ polygonHierarchy: new Cesium.PolygonHierarchy(cell.positions), perPositionHeight: true, vertexFormat: Cesium.PerInstanceColorAppearance.VERTEX_FORMAT }), attributes: { color: Cesium.ColorGeometryInstanceAttribute.fromColor( Cesium.Color.fromCssColorString(colorRamp(cell.value, min, max)) ) } }); }); const primitive = new Cesium.Primitive({ geometryInstances: instances, appearance: new Cesium.PerInstanceColorAppearance({ flat: true, translucent: false }) }); viewer.scene.primitives.add(primitive);

这里最关键的是perPositionHeight: true，它告诉 Cesium 使用每个点自带的弧度坐标作为高度基准，而不是拍扁到椭球表面。flat: true表示去掉光照和渐变，让每个色块看起来平整清晰。如果格网数量超过 10000，建议把温度场水平切成几块，分片创建多个 Primitive，每片控制在 5000 个小格以内，否则首次渲染还是会卡顿。这个问题我调了很久才反应过来，Cesium 合并 DrawCall 的核心就是尽量少创建 Primitive，而不是少创建格子。

5. 运行手册与避坑排查

5.1 环境准备：静态服务器与 Cesium 本地化

这个项目是纯前端，没有构建步骤，但绝对不能直接用file://打开 HTML，因为浏览器会拦截 CSV 和 JSON 的本地加载请求。最容易跑通的方式是起一个静态文件服务器，任选 Python 或 Node 都行。

# 在项目根目录执行 python -m http.server 8080

打开浏览器访问http://localhost:8080。如果浏览器控制台里 Cesium 加载报错，多半是config.js里的CesiumBaseUrl指向了 CDN，而当前网络访问不到。强烈建议把 Cesium 的 Build 目录放到项目libs/cesium/下，并设置CesiumBaseUrl = 'libs/cesium/'，这样离线也能跑。

5.2 从命令到看到温度场

启动静态服务器后，页面会自动加载data/stations.csv，并在页面右上角显示当前站点数量、格网分辨率和插值耗时。如果之前生成过data/grid.json且config.js里gridCacheEnabled为true，页面会直接渲染缓存，不会看到计算等待过程，这对演示场景非常友好。

# 首次运行时可以手动生成格网缓存 # 在浏览器控制台执行： generateGridFromCurrentData().then(grid => { console.log(JSON.stringify(grid)); });

把控制台输出的 JSON 复制到data/grid.json，之后刷新页面速度会快很多。注意这里的generateGridFromCurrentData是项目在 index.html 里暴露的调试函数，不是全局 API。如果你看到页面底部有 Cesium 的 ion Logo，那是没配 Token，不影响功能，但生产环境建议替换成自己的 Token。

5.3 高频雷区清单

以下是这个项目最容易遇到的问题，每一条我都亲自踩过，写出来帮你省调试时间。

现象一：插值结果出现一个个孤立的圆形色斑，看起来像“牛眼”。 原因：块金值c0太小，邻近样本点权重过大，离目标点稍远一点的点几乎不参与贡献。 解决：把config.js里的nugget参数从 0 调到0.1 * variance，再观察效果。

现象二：温度场整体偏移，或者部分色块跑到地球背面。 原因：经纬度字段顺序写反，lng当成lat用。 解决：检查 CSV 表头，统一为id,lng,lat,temp，并在loadStations解析时打印第一条记录核对。

现象三：页面白屏，控制台报Cesium is not defined。 原因：Cesium 静态资源路径错误，或 CDN 加载被拦截。 解决：确认libs/cesium/下存在Cesium.js，把CesiumBaseUrl指向本地相对路径。

现象四：温度场颜色几乎全部相同，分不出高低温。 原因：归一化时候min和max写成了固定值，或者数值数组中包含 NaN。 解决：在计算颜色前先过滤非法值，再用Math.min.apply(null, values)和Math.max求实际范围。

现象五：格网很密集时，拖动地图明显掉帧。 原因：每个小格网都创建了独立 Primitive 或 Entity，DrawCall 爆炸。 解决：改用GeometryInstance合并到单个 Primitive，同时把nx、ny控制在 100 以内，必要时开启gridCacheEnabled缓存。

6. 二次开发思路：换数据源、换算法、加时间轴

6.1 换数据：从温度到 PM2.5 或高程

这个项目的核心是插值和渲染，数据源可以非常轻松替换。比如你想做空气质量可视化，只需要把 CSV 里的temp列改成pm25，然后改colorRamp.js里的色带，把蓝青黄红换成适合空气质量的紫黄橙。插值过程完全不关心字段名，只读取value属性，所以你要做的只是调整loadStations里取值的那一行。高程数据也能用，但要注意单位，温度是摄氏度，高程是米，渲染时如果数值范围差异大，需要先做一次对数归一化。

6.2 换算法：用 IDW 做对照

克里金插值不是唯一选择，做实验时经常需要和反距离加权（IDW）对比。你可以在src下新增一个idw.js，实现核心公式value = Σ(wi * vi) / Σ(wi)，其中wi = 1 / distance²。然后在config.js里加一个algorithm字段，kriging或idw。由于generateGrid依赖策略模式，你只需要在插值函数入口做一个分支即可，不用动渲染代码。这个对照实验非常适合验证：克里金到底比 IDW 好在哪里，是不是值得多付出的计算时间。

6.3 动态时间序列与交互

温度场最常见的二次需求是按时间轮播，展示一天内不同小时的变化。这里不推荐 setInterval 直接修改 Primitive 的颜色，因为 Cesium 会非常频繁地重绘。我一般会预生成多个格网 JSON，每个时间点一份，然后用时间轴作为驱动，在clock的回调里替换 Primitive 的geometryInstances。如果不要求平滑过渡，简单做法是先移除旧 Primitive，再添加新 Primitive，每次切换控制在几百毫秒内。记得给每个时间帧做一次预计算，不要在播放时实时插值，否则 CPU 会瞬间打满。

从那次临时 Demo 开始，我养成一个习惯：任何插值类可视化项目，第一次跑通后都会强制检查一次「最小分辨率 + 最大采样点」下的性能下限，再往上加复杂度。特别是克里金这种算法，看似只是矩阵求解，真正卡住你的往往不是算法本身，而是离站点最近的几个点因为距离重合导致矩阵奇异。希望这份解析能帮你在自己的地图可视化里少踩几个坑。

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