上篇结尾我留了这么一句话：
“下一篇我们先搁下地图控件这口气，去算一件要真动手的事——在一片地形上，挖掉多少、填上多少。空间分析-填挖方分析，我们下一篇见。”
前面几篇我们一直在做"看得见"的东西：遮罩、卷帘、比例尺，本质上都是把引擎已有的能力搬到一个能抓的地方。从今天起换一条线——空间分析：这一回不是把数据画出来，而是从数据里算出一个数，而且这个数要能拿去报价、拿去施工。
填挖方（cut & fill）是这类问题里最经典的一个：给你一片地形和一个设计标高，问"要挖掉多少立方米、要填上多少立方米"。听起来简单，可真动手你会发现，每一步都有一个可以悄悄算错的地方。
大家好，我是 Cesium 酱（也可以叫我"本猿"），一名在 WebGIS 领域摸爬滚打多年的前端开发者。
💡 本文涉及的全部代码均已开源，完整的仓库地址见文末，可自行取用、随意魔改。
一、场景与痛点：一片地上到底要动多少土
🎯 背景现状
任何一块场地动工前，甲方和施工方都要先对一件事达成共识：这片地要挖掉多少、填上多少。这个数字直接决定土方外运的车次、买土的成本、工期，甚至决定"这块地到底值不值得拿"。
传统做法是测量员跑现场、打方格网、用水准仪读高程，回来在 CAD 里算。WebGIS 把它搬到浏览器里，靠的是三样现成的东西：一个多边形区域 + 一份地形高程 + 一个设计标高。
常见的三种偷懒做法：
- 拿整个区域的平均高差乘面积——
V = S × Δh，一秒钟出结果，也一秒钟就错：地形起伏被抹平了，挖和填会互相抵消； - 只算挖不算填（或反过来）—— 报出来的数永远只有一半；
- 用"格网法"但格子粗到离谱—— 一个 100 m 的场地只切 4 格，填挖分界线根本对不上。
⚠️ 难点到底在哪
难在"填"和"挖"的分界不是画出来的，是算出来的：它是"地形面 = 设计面"那条等高线。这条线落在哪儿，决定了哪部分算挖、哪部分算填。而计算机只能按离散的采样点去逼近它——于是问题就从"算体积"变成了"你的网格有多细"。
二、原理拆解：把地形切成三棱柱，再一个一个加起来
2.1 一句话的原理
把多边形切成三角网，每个三角形配一个三棱柱：柱底是设计面、柱顶是地形面，柱子的体积就是这一块的填方或挖方。
这就是"三棱柱法"（也叫 TIN 体积法）。它比方格网法更贴地形，因为三角形的三个顶点直接踩在采样点上，顶面是一片斜平面，而不是"一格取一个代表高度"。
2.2 三角网从哪来：PolygonGeometry与 granularity
采样网格不是我们自己写的，是问 Cesium 要的：
constgranularity=Math.PI/Math.pow(2,11)/precision.valueconstgeometry=PolygonGeometry.fromPositions({positions,vertexFormat:PerInstanceColorAppearance.FLAT_VERTEX_FORMAT,granularity})constgeom=PolygonGeometry.createGeometry(geometry)granularity是 Cesium 里多边形细分的角度粒度（单位：弧度）。它默认取CesiumMath.RADIANS_PER_DEGREE，也就是1°——对洲际尺度的多边形刚好，对我们这种几百米的场地就等于"不切"。
本例的演示区域是写死在代码里的 4 个顶点，本猿把它们换算回经纬度：
| 项 | 实测值 |
|---|---|
| 经度范围 | 117.170153°E ~ 117.170879°E |
| 纬度范围 | 34.257988°N ~ 34.258899°N |
| 顶点椭球高 | 4.447 ~ 4.449 m |
| 周长 | 328.76 m |
| 面积 | 6480.2 ㎡（三角化后 6464.2 ㎡） |
把四档精度各跑一遍，网格是这样长的：
| 采样精度 | granularity | 顶点 | 三角形 | 平均三角形面积 | 平均边长 |
|---|---|---|---|---|---|
| 256（默认） | 5.9921e-6 rad | 19 | 24 | 269.34 ㎡ | ≈16.41 m |
| 512 | 2.9961e-6 rad | 61 | 96 | 67.34 ㎡ | ≈8.21 m |
| 1024 | 1.4980e-6 rad | 217 | 384 | 16.83 ㎡ | ≈4.10 m |
| 2048 | 7.4901e-7 rad | 817 | 1536 | 4.21 ㎡ | ≈2.05 m |
| 对照：默认 1° | 0.01745 rad | — | 2 | — | — |
三角形数是按精度的平方涨的（24 → 96 → 384 → 1536，每档整 4 倍），因为 granularity 跟1 / precision成正比，而细分是二维的。
这里有个值得单独说的认知："采样精度 256"听起来很高，实际只切出 24 个三角形。那个 256 不是"切 256 份"，它是分母上的一个系数。要判断网格够不够细，只有"平均三角形面积"或"平均边长"这两个数说了算——本例默认档是16.4 m 一格，量一个 60×100 m 的场地，只有 24 块。
2.3 高度从哪来：两条完全不同的路
网格建好之后要逐顶点取高程，代码里有个很扎眼的判断：
constterrainProvider=viewer.terrainProviderasunknownas{_layers?:unknown[]}|undefinedif(terrainProvider?._layers){constsampled=awaitsampleTerrainMostDetailed(viewer.terrainProvider,cartos)heights=sampled.map((c)=>c.height)}else{heights=cartos.map((c)=>viewer!.scene.globe.getHeight(c)??0)}_layers是CesiumTerrainProvider的私有字段（CesiumTerrainProvider.jsL110provider._layers = this.layers）。这里靠嗅探私有字段来判断"是不是量化网格地形"，能跑但很脆：换一个没这个字段的地形源，就悄悄走另一条路。- 有
_layers→sampleTerrainMostDetailed：联网去取最高精度的高程，代价是异步、要等、可能失败。 - 没有 →
globe.getHeight(c)：读当前已经加载进内存的瓦片的高程，瓦片没加载就返回undefined，再被?? 0兜成0 米。
两条路的精度差着好几个数量级，而切换是静默的。这一点在第五节会变成一个真问题。
2.4 体积怎么算：三顶点均值的那一次"精确"
constmean=(ha+hb+hc)/3if(mean<design){fillArea+=area;fillVolume+=area*(design-mean)}elseif(mean>design){cutArea+=area;cutVolume+=area*(mean-design)}else{noArea+=area}mean = (ha + hb + hc) / 3这一步是精确的，不是近似。三个顶点撑起一个平面，平面在三角形重心处的高度恰好等于三个顶点高度的算术平均；而斜顶棱柱的体积恰好等于底面积乘重心高。
本猿用一个纯平面地形（顶面 = 一片斜平面）验证了一遍：
| 项 | 实测值 |
|---|---|
| 三棱柱累加 | 172020.162631 m³ |
| 解析真值 | 172020.162631 m³ |
| 相对误差 | 0.000e+0 |
所以"三棱柱法"本身没问题，误差全在别的地方。
2.5 反直觉：填挖分界线是被"咬"出来的
注意那三个分支——一个三角形只看均值，要么整块算填、要么整块算挖。可真实世界里，一个三角形常常一半在填区、一半在挖区。
于是那条"地形 = 设计面"的等高线，被强行吸附到了三角形的边上：网格越粗，锯齿越大。本猿用一片合成地形（斜坡 + 正弦脊）做了一组对照，真值用 1200×1200 的细网格积出来：
| 采样精度 | 三角形 | 平均边长 | 挖方体积偏差 | 填方体积偏差 |
|---|---|---|---|---|
| 256 | 24 | 16.41 m | −2.59% | −2.59% |
| 512 | 96 | 8.21 m | −0.96% | −0.96% |
| 1024 | 384 | 4.10 m | −0.25% | −0.25% |
| 2048 | 1536 | 2.05 m | −0.08% | −0.08% |
网格每加密一倍，误差大约降到原来的三分之一——这是一阶收敛，因为误差主要来自那条边界带（面积 ≈ 边界长度 × 格子边长）。换句话说：填挖方量的精度，是被"分界线有多长"和"格子有多大"这两个数锁死的，跟"总面积多大"反而关系不大。
顺带说一句那个else：只有mean === design才会走它。高程是浮点数，精确相等几乎不可能发生——所以面板上的"无须填挖面积"平时就是 0，它真正有意义的场景见第五节。
三、代码实战：从多边形到六个数字
3.1 主流程
asyncfunctionrunAnalysis():Promise<void>{constpositions=manualPositions.length>=3?[...manualPositions]:ORIGINAL_POSITIONSisAnalyzing.value=truestatusMessage.value=''awaitnewPromise((resolve)=>setTimeout(resolve,50))// 让出一帧，把"分析中…"画出去try{// ① 切网格 ② 取高程 ③ 定设计标高 ④ 逐三角形累加}catch(error){statusMessage.value=errorinstanceofError?error.message:String(error)}finally{isAnalyzing.value=false}}那句await new Promise(r => setTimeout(r, 50))值得单独说：它不是在等什么数据，而是故意让出 50 毫秒，好让isAnalyzing变成true之后能被 Vue 渲染成"分析中…"。不加这一句，同步的密集计算会把状态变更一起憋住，用户点下去毫无反应。这是所有"同步重活"都要补的一脚刹车。
3.2 设计标高：填 0 的时候发生了什么
constmeanHeight=heights.reduce((sum,h)=>sum+h,0)/heights.lengthconstdesign=designHeight.value>0?designHeight.value:meanHeight effectiveDesignHeight.value=Number(design.toFixed(2))面板上写着"填 0 时自动采用采样区域的平均地形高度"，这看起来是个贴心功能：平均高度就是填挖平衡的那个高度（挖出来的正好填回去，净方为零）。
但注意它是顶点高度的算术平均——顶点在三角网里被重复计数：内部顶点被 6 个三角形共用，边界顶点只被 2~3 个用。所以它不是面积加权平均，也就不是真正的平衡高度。
本猿用合成地形实测（precision = 1024）：
| 设计标高取法 | 数值 | 挖方 | 填方 | 净方 |
|---|---|---|---|---|
| 顶点算术平均（demo） | 6.1006 m | 9966.7 m³ | 10132.3 m³ | −165.7 m³ |
| 面积加权平均（真平衡） | 6.0750 m | 10049.7 m³ | 10049.7 m³ | 0.0 m³ |
两者只差0.0256 m，却留下了165.7 m³ 的净方（约占单侧方量的 1.65%）。工地上的土方是按车算钱的，1.65% 不是小数。
另外那个designHeight.value > 0也值得留意：填 0 会被当成"没填"，所以你没法把设计标高设成 0 米（想设成 0，请用 0.01）。负数同理被忽略。
3.3 体积怎么变成看得见的三棱柱
constinstances=selected.map((solid)=>newGeometryInstance({geometry:newPolygonGeometry({polygonHierarchy:newPolygonHierarchy(solid.top),perPositionHeight:true,extrudedHeight:effectiveDesignHeight.value,vertexFormat:PerInstanceColorAppearance.VERTEX_FORMAT}),attributes:{color:ColorGeometryInstanceAttribute.fromColor(colorForType(solid.type))}}))volumePrimitive=viewer.scene.primitives.add(newPrimitive({geometryInstances:instances,appearance:newPerInstanceColorAppearance({flat:true,translucent:true})}))perPositionHeight: true+extrudedHeight：顶面按三个顶点的真实高程起伏，底面统一压在设计标高上——正好就是那个三棱柱。- 挖方是橙色
#ff8c2e、填方是绿色#37c870、无须填挖是灰色#8a9bb0。 - 一个
Primitive装下全部实例，靠PerInstanceColorAppearance逐实例上色（这一套和 #59~#61 那几篇的批量实例是同一条路）。 - 实测单个三棱柱：18 个顶点 / 24 条索引 / 8 个三角面。于是 2048 档一次要建1536 个实例、27648 个顶点，顶点 + 索引约1.09 MB。
⚠️ 这里有个小浪费：外观写的是flat: true，而PerInstanceColorAppearance在flat: true时自己的顶点格式就是FLAT_VERTEX_FORMAT（只要 position）（PerInstanceColorAppearance.jsL82-84）。可上面给几何传的是VERTEX_FORMAT（position + normal，36 字节/顶点）——法线算了、传了、着色器一眼都不看。2048 档下这是27648 × 12 = 331,776 字节，约 0.32 MB 的死重量。有意思的是同一份代码里采样网格用的FLAT_VERTEX_FORMAT是对的，三棱柱这一处用错了。
3.4 交互：三次点击定一个区域
handler.setInputAction((event)=>onLeftClick(event),ScreenSpaceEventType.LEFT_CLICK)handler.setInputAction((event)=>onMouseMove(event),ScreenSpaceEventType.MOUSE_MOVE)handler.setInputAction(()=>finishDrawing(),ScreenSpaceEventType.RIGHT_CLICK)handler.setInputAction(()=>finishDrawing(),ScreenSpaceEventType.LEFT_DOUBLE_CLICK)左键一个点一个点地采，右键或双击闭合。采点用的是pickPosition（measure-lib/pick.ts），它有个很实用的三级回退：scene.pickPosition（深度拾取，最准）→globe.pick（地形射线）→camera.pickEllipsoid（椭球解析）。深度拾取要 GPU 支持（scene.pickPositionSupported），不支持时自动降级——这套三级回退值得抄到任何需要"点哪落哪"的功能里。
预览也有讲究：>= 3个点时画半透明多边形，不足 3 个点时退化成clampToGround的折线，两种状态切换得很自然。
四、运行效果
初始视角是这么设的：
viewer.camera.setView({destination:ORIGINAL_POSITIONS[0].clone()})只给了 destination，没给 orientation——按Camera.jsL1528-1530，默认姿态是heading = 0, pitch = -π/2, roll = 0，也就是镜头笔直朝下。而destination是相机自己的位置，不是看向的目标（Camera.jsL1493-1501：只有传Rectangle才会被当成"看向的范围"）。
于是相机被放在第 1 个顶点上，椭球高只有4.447 m，正对着脚下的地皮俯视。这里本猿要提醒一句：这个视角相当危险——相机离椭球面才 4.4 米，只要当地地形高于这个数（绝大多数地方都会），相机一上来就在地面以下。想要一个能看的初始视角，请给orientation，或者干脆用viewer.zoomTo(entity)/ 一个Rectangle。
面板上那六个数字是这么来的：
| 面板项 | 含义 | 累加方式 |
|---|---|---|
| 总分析面积 | cut + fill + no | 逐三角形面积求和 |
| 填方面积 / 填方体积 | 地形低于设计面 | area、area × (design − mean) |
| 挖方面积 / 挖方体积 | 地形高于设计面 | area、area × (mean − design) |
| 无须填挖面积 | 恰好相等 | area（平时为 0） |
五、避坑与优化
5.1 精度不是越高越好，但默认档一定不够
三档数字摆一起看：
| 采样精度 | 三角形 | 体积误差 | 可视化顶点 | 缓冲 |
|---|---|---|---|---|
| 256（默认） | 24 | −2.59% | 432 | 0.02 MB |
| 512 | 96 | −0.96% | 1728 | 0.07 MB |
| 1024 | 384 | −0.25% | 6912 | 0.27 MB |
| 2048 | 1536 | −0.08% | 27648 | 1.09 MB |
工程上一般要求土方量误差≤ 3%，所以 256 档在"合规线"上勉强够用，但你要是拿它去对账，最好先说明分界线被切成了 16 米一粒的锯齿。真要较真，1024 档（0.25%）是个不错的平衡点。别把"精度"当成一个可以无脑拉满的旋钮——它是平方级的代价。
5.2 采样失败时，界面上什么都看不出来
这是本篇最该记住的一条。两条高程来源，各有各的"静默失败"：
①sampleTerrainMostDetailed失败 →undefined
它的第三个参数rejectOnTileFail默认false（sampleTerrainMostDetailed.jsL42-44），也就是瓦片请求失败时不抛错，只把 height 留成undefined。接着mean = (ha+hb+hc)/3变成NaN，而NaN < design与NaN > design都是 false——于是这个三角形走进else分支，被记成**“无须填挖面积”**。
本猿模拟了一遍（全部高程 undefined）：挖方 0 ㎡、填方 0 ㎡、“无须填挖” 6464.2 ㎡（占总面积 100%），体积显示 NaN。界面上没有任何报错，只是多了一个看起来很合理的"无须填挖面积"。
顺带一提，它还会在terrainProvider.availability未定义时直接 reject（源码 L57-61）——那种情况反而好，会被catch住显示出来。
②globe.getHeight未加载 →?? 0
这条路更狠：瓦片还没加载完时返回undefined，被?? 0兜成0 米。于是"没有数据"变成了"这里是海平面"，如果该地实际高程是 300 米，每个三角形都会多算 300 米的挖方。
修法：把"有没有拿到高程"当成一等公民。
constsampled=awaitsampleTerrainMostDetailed(viewer.terrainProvider,cartos,true)// 失败就抛// 回退分支里也别用 ?? 0consth=viewer.scene.globe.getHeight(c)if(h===undefined){missCount++;continue}至少要做到：统计缺失点数，超过阈值就别出结果，或者先await viewer.scene.globe.tilesLoaded/ 等地形加载完再放行（本例的loadWorldTerrain已经有 await，但用户在加载完成前点按钮仍会踩到）。
5.3 设计标高：要么让用户输入，要么用面积加权平均
2.2 节那 0.0256 m 的差别值 165.7 m³。改法很简单——加权时用面积而不是顶点：
letwSum=0,aSum=0for(consttoftris){constA=area(t);constm=mean(t);wSum+=A*m;aSum+=A}constdesign=aSum>0?wSum/aSum:meanHeight// 这才是真正的平衡标高顺便把designHeight.value > 0的判断改成"是否填了值"（比如用??或加一个useAutoHeight开关），否则用户永远设不了 0。
5.4 顺手能做的几件小事
vertexFormat改回FLAT_VERTEX_FORMAT（外观已经是flat: true），省掉 0.32 MB 的死法线。- 别嗅探
_layers：改用typeof (viewer.terrainProvider as any).sampleTerrain === 'function'或者直接判断availability是否存在，私有字段随时会变。 PerInstanceColorAppearance的translucent默认就是 true，半透明体之间不写深度、也不排序，密密麻麻的柱体会互相穿。想要扎实的效果可以把填挖体分开两个不透明 Primitive（按类型筛选已经有了，只差一个translucent: false）。- 单位与量纲要写清楚：面板写的是 ㎡ 和 m³，输入的设计标高是米；跨单位项目（英尺/码）请在进口处统一换算，别在体积公式里改。
- 大面积区域要分块：本例只有 6480 ㎡；真来一块几平方公里的场地，网格顶点数会平方级爆炸（参考 #59~#62 那几篇的规模算账），记得做分块 + 进度条。
5.5 上线前 checklist
- ✅ 网格密度用平均三角形面积 / 平均边长说话，别信"精度 256"这种数字
- ✅ 高程采样失败要显式报错，不要让
undefined流进mean，也不要用?? 0兜底 - ✅ 自动设计标高用面积加权平均，否则"填挖平衡"只是看起来平衡
- ✅ 同步重活前让出一帧（
setTimeout 50），否则状态变更渲染不出来 - ✅ 外观
flat: true时几何就用FLAT_VERTEX_FORMAT - ✅ 别依赖私有字段（
_layers）做分支 - ⚠️ 记得在结果旁写一句：填挖分界线被网格咬成了锯齿，精度档位是多少
六、写在最后
Cesium 是头优雅的大象。
我们曾经以为，算体积的难点是"体积公式"；
今天才发现，公式本身精确到了小数点最后一位——难的是承认那些没量到的地方：分界线被咬成了锯齿，没取到的高程变成了 0 米，而"自动平衡"其实差了 165 立方。
这一篇最让本猿坐不住的，是同一个算法里的两种"精确"：
第一，(ha + hb + hc) / 3是精确的。三棱柱体积公式对斜顶面一点误差都没有（实测相对误差 0.000e+0）。数学这一半，干净利落。
第二，"这块算填还是算挖"是糙的。一个三角形只凭均值二选一，于是那条等高线被强行掰到网格线上。精度 256 档下体积差2.59%——误差不来自公式，来自离散化。
第三，最危险的不是算错，是算错了还长得对。采样失败时undefined静默变成"无须填挖面积"，未加载的瓦片静默变成 0 米高程；两种情况下界面都端出一个完整、自洽、可以用来汇报的数字。一个会报错的错误只是麻烦，一个不报错的错误才是事故。
这个案例还不完美：初始视角贴在地皮上、设计标高填不了 0、精度档位的名字有误导性、法线白算了 0.32 MB。但它把"一片地要动多少土"这件事从公式到网格到数字，整条链路都摊开讲了一遍，这就够了。🎉
下一篇我们回到地图控件那条线上——这次不是量距离，而是告诉用户"哪边是北"。罗盘控件，我们下一篇见。
如果你也在项目里做过土方量计算，欢迎留言说说你是怎么处理"填挖分界线被网格咬碎"这件事的——是拉高网格精度硬扛，还是先勾出等高线再分区算？
🌟 项目已开源，欢迎 Star！（你的 Star ⭐ 就是对我最大的支持！）
仓库地址：https://github.com/jianlei-wang/CESIUM-TREASURE-BOX
代码在src/cases/cut-fill-analysis/目录下，核心就一个CutFillAnalysisDemo.vue，把ORIGINAL_POSITIONS换成你自己的场地经纬度就能跑。
Cesium酱の百宝箱 · 第 67 篇
三棱柱量尽高低土，一条界线只认网格那一格
2026 年的秋夜，与你共探"没量到的地方该怎么被说出来"
