返回知识库

GIS

四叉树与 LOD

四叉树与 LOD 封面
GISengine-webgpu四叉树LOD

前情回顾02 坐标与投影 讲清了经纬度怎么变成 ECEF / Web 墨卡托米。本篇把这些”投影后的米”切成 256×256 的小方块,按相机位置动态调度——这就是你在 Web 地图上”无限缩放”看到的画面背后的算法。

直觉问题

打开 Google Maps,做两个动作:

  1. 滚轮放大到北京天安门,注意观察——画面会先模糊后清晰,中间有约 200 ms 的”瓦片加载”过程。这中间发生了什么?为什么不是一次性显示清晰画面?
  2. 拖动地图往任意方向滑动——新画面几乎是瞬间出现的,没有重新加载。但你的浏览器不可能预下载了整个地球的瓦片(z=18 全球有 680 亿张)。它是怎么知道接下来要画哪几张的?

这两个问题的答案合起来就是 Web GIS 引擎最核心的子系统:瓦片调度。它包含 4 个关键算法:

  • 瓦片编号(z/x/y 怎么编出来)
  • 四叉树空间索引(如何快速找当前视野的瓦片)
  • LOD 选择(用哪一层细节)
  • 剔除算法(视锥 + 地平线剔除,3D 地球专属)

读完这一篇,你能回答:“为什么瓦片是 256×256 不是 512×512?”、“为什么 TMS 和 OSM 编号 y 轴相反?”、“3D 地球为什么需要地平线剔除?“

核心概念白话讲

用”金字塔”理解瓦片层级

把整个地球的地图按缩放级别切成 23 层(z=0 到 z=22),就是瓦片金字塔

z 层全球瓦片数单瓦片覆盖典型用途
01整个地球世界地图
14一个半球大洲概览
364一个中型国家国家地图
716 384一个省省级导航
10约 100 万一座大城市城市导航
15约 10 亿几条街道街道地图
18约 680 亿一栋楼街景细节

关键公式:第 zz 层瓦片数为 4z4^z

为什么是 4 倍?因为每升一级,地图宽度和高度都翻倍,所以面积变 4 倍。这就是”四叉树”分裂的来源——每个瓦片分裂为 4 个子瓦片。

用”乐高积木”理解四叉树

每个瓦片是 256×256 像素的小图。升一级时,1 个父瓦片均匀分裂为 4 个子瓦片(2×2 网格):

父瓦片 z=N (1 张)              子瓦片 z=N+1 (4 张)
┌──────────────┐              ┌──────┬──────┐
│              │              │ 0,0  │ 1,0  │
│              │              ├──────┼──────┤
│              │              │ 0,1  │ 1,1  │
│              │              │      │      │
└──────────────┘              └──────┴──────┘

整个地球的瓦片集合构成一棵四叉树:根节点(z=0)有 4 个子节点(z=1),每个子节点又有 4 个子节点(z=2)……递归 22 层。

为什么不八叉树?地球表面是 2D 流形,不是 3D 体。3D 场景(如体素游戏 Minecraft)才用八叉树,2D 地图用四叉树就够。

为什么不用 KD-Tree / BVH?这些是动态对象空间索引(如游戏中几百个 NPC)。瓦片是静态空间划分——四叉树完美匹配瓦片的”均匀 4 分裂”几何。

用”舞台聚光灯”理解视锥剔除

3D 地球视图下,相机只能看到地球表面的 1/3 左右。视锥剔除(Frustum Culling)就是”聚光灯外的不画”:

  • 把相机视锥体定义为 6 个裁剪面(左/右/上/下/近/远)围成的截头金字塔
  • 每个瓦片计算自己的包围盒(AABB 或球)
  • 测试包围盒是否与视锥相交——不相交就跳过

效果:典型 3D 地球视图中,全球 10910^9 量级的瓦片,视锥剔除后只剩几百张要画。

用”地球仪的背面”理解地平线剔除

即使瓦片在视锥内(朝向相机那一面),它可能位于地球背面——被地球本体挡住了。地平线剔除(Horizon Culling)就是”看不到的不画”:

          相机

            \
             \  ← 视线
              \
   ─────────────●─── 地球切面(地平线)
              ╱ ╲
             ╱   ╲
            ╱  地  ╲
        隐藏│  球  │隐藏
           ╱       ╲
          ╱         ╲

地平线之上的瓦片(朝相机那半球)可见;地平线之下的瓦片(背半球)不可见。这一步可减 30-50% 瓦片。

NOTE

关键认知:视锥剔除 + 地平线剔除是串联的——先视锥筛掉视锥外的,再地平线筛掉背面的。两者一起把”全球瓦片”筛到”当前可见瓦片”。

原理与数学/机制

1. z/x/y 编号数学

瓦片编号 z/x/y 是 Web 地图通用约定:

  • zz = 缩放级别(0-22)
  • xx = 该层的列号(从左到右),范围 [0,2z1][0, 2^z - 1]
  • yy = 该层的行号(从上到下),范围 [0,2z1][0, 2^z - 1]

Web 墨卡托瓦片(Google/OSM/Bing 标准)

从经纬度 (λ,φ)(\lambda, \varphi)z/x/yz/x/yzz 已选定):

{n=2zx=nλ+180360y=n1ln(tanφ+secφ)/π2\begin{cases} n = 2^z \\ x = \left\lfloor n \cdot \dfrac{\lambda + 180}{360} \right\rfloor \\ y = \left\lfloor n \cdot \dfrac{1 - \ln(\tan\varphi + \sec\varphi)/\pi}{2} \right\rfloor \end{cases}

z/x/yz/x/y 反求经纬度(瓦片左上角):

{λ=x2z360180φ=arctan(sinh(π(12y/2z)))180π\begin{cases} \lambda = \dfrac{x}{2^z} \cdot 360 - 180 \\ \varphi = \arctan(\sinh(\pi \cdot (1 - 2y/2^z))) \cdot \dfrac{180}{\pi} \end{cases}

EPSG:4326 瓦片方案(geodetic / WMTS 标准)

EPSG:4326 不用 Web 墨卡托,直接用经纬度切瓦片。与 Web 墨卡托瓦片的差异:

维度Web 墨卡托瓦片EPSG:4326 瓦片
投影EPSG:3857(球面墨卡托)EPSG:4326(直接经纬度)
单瓦片形状正方形(赤道与极地相同)矩形(纬度越高越窄)
y 编号方向从上到下(北极 = 0)从下到上(南极 = 0,TMS 风格)
z=0 瓦片数1(整个地球)2(东西半球各一)
主要使用者Google/OSM/Bing/百度WMTS OGC 标准、NASA GIBS

EPSG:4326 瓦片正变换(与 Web 墨卡托不同):

{x=2zλ+180360yTMS=2z+1φ+90180\begin{cases} x = \left\lfloor 2^z \cdot \dfrac{\lambda + 180}{360} \right\rfloor \\ y_{\text{TMS}} = \left\lfloor 2^{z+1} \cdot \dfrac{\varphi + 90}{180} \right\rfloor \end{cases}

注意 yTMSy_{\text{TMS}}南极开始计数(与 OSM 相反)。换算关系:

yOSM=2z1yTMSy_{\text{OSM}} = 2^z - 1 - y_{\text{TMS}}

2. 四叉树空间索引

四叉树调度伪代码:

输入: 当前相机位置 + 视锥 + LOD 阈值
输出: 本帧需要绘制的瓦片集合

函数 traverse(node):
  if 视锥剔除(node) == 完全在外面:
    return  # 跳过整个子树
  if 地平线剔除(node) == 完全在背面(仅 3D):
    return
  if LOD 选择(node) == 足够细:
    加入绘制集合
  else:
    for child in node.children:  # 递归 4 个子瓦片
      traverse(child)

从根节点 z=0 开始 traverse

关键优化

  • 早期剔除:父节点视锥/地平线剔除失败时,整个子树都跳过
  • 逐层细化:每个节点单独判断 LOD,不必全局统一层级
  • 包围盒层级继承:父节点包围盒 = 4 个子节点包围盒的并集

3. LOD 选择策略

LOD(Level of Detail)决定一个瓦片”够不够细”。两种主流策略:

3.1 距离阈值(2D 简单版)

最简单的 LOD:基于相机到瓦片的距离 dd 与瓦片层级 zz 的固定映射。

z=clamp ⁣(log2 ⁣(R2πcosφdtilePixels),0,zmax)z = \text{clamp}\!\left(\left\lfloor \log_2\!\left(\dfrac{R \cdot 2\pi \cdot \cos\varphi}{d \cdot \text{tilePixels}}\right) \right\rfloor, 0, z_{max}\right)

其中 RR 是地球半径,φ\varphi 是瓦片纬度,tilePixels=256\text{tilePixels} = 256

优点:计算简单,2D 地图完全够用。 缺点:精度差,“画面变化”和”瓦片切换”不严格对齐。

3.2 屏幕空间误差 SSE(3D 精确版)

3D 引擎(如 Cesium)用屏幕空间误差(Screen Space Error):

SSE=ϵHfd2tan(θ/2)\text{SSE} = \dfrac{\epsilon \cdot H \cdot f}{d \cdot 2\tan(\theta/2)}

符号说明

  • ϵ\epsilon — 瓦片的几何误差(geometric error,该层瓦片与”完美细节”的差距,米)
  • HH — 屏幕高度(像素)
  • ff — 视锥缩放因子(透视投影 near/far 比值)
  • dd — 相机到瓦片距离(米)
  • θ\theta — 相机垂直 FOV(弧度)

判定规则:如果 SSE>SSEthreshold\text{SSE} > \text{SSE}_{\text{threshold}}(典型 16 像素),瓦片”太粗”,递归到子瓦片;否则当前瓦片够细,停止分裂。

直觉ϵ\epsilon 是”如果用粗瓦片代替细瓦片,会差多少米”;乘以”屏幕每米对应的像素数”,就是”屏幕上看起来差多少像素”。差超过 16 像素,肉眼能看见,必须细化。

IMPORTANT

关键不变量:SSE 是 3D Tiles 标准Phase 9 详讲)的 LOD 核心。2D Web 瓦片(如 OSM)通常不需要 SSE,因为相机始终正交俯视,距离阈值足够准确。

4. 视锥剔除的数学

相机视锥体是 6 个裁剪面围成的截头金字塔。每个面对应一个不等式:

insidei(p)=nip+di0,i{L,R,B,T,N,F}\text{inside}_i(\vec{p}) = \vec{n}_i \cdot \vec{p} + d_i \geq 0, \quad i \in \{L, R, B, T, N, F\}

其中 ni\vec{n}_i 是第 ii 个面的法向量(指向视锥内部),did_i 是常数项。

瓦片测试:计算瓦片包围盒的 8 个顶点,全部满足 6 个不等式 → 完全在视锥内;全部不满足任一不等式 → 完全在外面;否则部分相交。

优化:用球测试代替包围盒测试——只测瓦片中心点和半径,计算量从 8 个顶点 × 6 个面 = 48 次点积,降到 6 次点积。

5. 地平线剔除的数学

3D 地球特有的剔除。原理:把瓦片顶点变换到以地球中心为原点的坐标系,判断它是否被地球遮挡。

设瓦片顶点向量 v\vec{v}(从地心出发),相机向量 c\vec{c}(从地心出发),则瓦片被遮挡当且仅当:

vcv2c2R2c2<0\vec{v} \cdot \vec{c} - |\vec{v}|^2 \cdot \dfrac{|\vec{c}|^2 - R^2}{|\vec{c}|^2} < 0

其中 RR 是地球半径。

直觉推导:从相机往地球切线方向画线,切线以下的部分(地平线下)被遮挡。点积 vc\vec{v} \cdot \vec{c} 衡量”瓦片在相机方向上的投影”,配合阈值即可判断。

TIP

Cesium 的实现技巧:用** cone-culling**(圆锥剔除)变种——把整个背面半球用一个圆锥近似,瓦片顶点点积与圆锥半角比较。计算量更小,精度略损但视觉无差异。

6. 瓦片调度流程

图 1:瓦片调度主线流程(虚线未画,A→B→...→K 是同步主循环,H→J 异步)

关键时序

  1. 同步(主循环每帧):A → E,得到本帧 draw list
  2. 异步(Web Worker):F → J,新瓦片下载/解码
  3. 过渡:在 J 完成前,用父瓦片(粗 LOD)兜底,子瓦片就绪后淡入

可视化对比与动手实验

对比表一:TMS vs Google/OSM 编号

维度TMS(Tile Map Service)Google/OSM/Slippy
y 轴方向从下到上(南极 = 0)从上到下(北极 = 0)
z=0 起点左下角左上角
URL 形式/tile/{z}/{x}/{y}.png/{z}/{x}/{y}.png/{z}/{x}/{y}.jpg
标准来源OSGeo OSGeo-WMSOpenStreetMap Wiki
主要使用者老 OGC 服务、NASA GIBSGoogle/OSM/Bing/Mapbox

转换公式yOSM=2z1yTMSy_{\text{OSM}} = 2^z - 1 - y_{\text{TMS}}

对比表二:Web 墨卡托瓦片 vs EPSG:4326 瓦片

维度Web 墨卡托瓦片(EPSG:3857)EPSG:4326 瓦片(geodetic)
投影球面墨卡托直接经纬度(无投影)
瓦片形状正方形(256×256)矩形(宽 256,高随纬度变化)
z=0 瓦片数1 张(全球)2 张(东西半球)
yy 编号方向OSM 风格(北极=0)TMS 风格(南极=0)
两极处理截断在 ±85.05°完整覆盖到 ±90°
主要使用者Google / OSM / Bing / 百度 / 高德WMTS OGC 标准 / NASA GIBS / 部分政府服务
适用场景Web 浏览器地图(视觉一致)科学数据交换(坐标精确)

对比表三:剔除算法对比

算法2D / 3D数学复杂度平均剔除率必需性
视锥剔除两者都用6 个面点积60-80%必需
地平线剔除仅 3D 地球向量点积 + 阈值30-50%3D 必需,2D 不用
距离剔除仅 2D 地图简单阈值5-10%可选优化
遮挡剔除两者都用复杂(Hi-Z depth buffer)10-20%可选优化

对比表四:LOD 策略对比

策略计算复杂度精度适用场景典型引擎
距离阈值低(log 计算)2D Web 地图Leaflet / OpenLayers
SSE 屏幕空间误差中(含几何误差)3D 地球 / 3D TilesCesium / engine-webgpu
混合策略最高高端 3D 引擎Unreal Engine Lumen

动手实验一:z/x/y 手算(北京天安门)

北京天安门 (116.397°E, 39.908°N),求 z=10z=10 的 Web 墨卡托瓦片编号:

  • n=210=1024n = 2^{10} = 1024
  • x=1024(116.397+180)/360=842.4=842x = \lfloor 1024 \cdot (116.397 + 180) / 360 \rfloor = \lfloor 842.4 \rfloor = 842
  • yytan(39.908°)+sec(39.908°)0.835+1.300=2.135\tan(39.908°) + \sec(39.908°) \approx 0.835 + 1.300 = 2.135ln(2.135)0.758\ln(2.135) \approx 0.758(10.758/π)/20.379(1 - 0.758/\pi)/2 \approx 0.379y=10240.379=388y = \lfloor 1024 \cdot 0.379 \rfloor = 388

结果:北京天安门在 z=10z=10 时瓦片编号为 (10, 842, 388)

验证:浏览器打开 https://tile.openstreetmap.org/10/842/388.png 应看到北京区域。

动手实验二:MapTiler 在线工具

打开 MapTiler - tile bounds

  1. 搜索任意地点
  2. 调整 z 滑块,看瓦片边界如何变化
  3. 注意 zz 每升一级,瓦片数 ×4,单瓦片地理范围 ÷4

动手实验三:在纸上画四叉树

在纸上画 z=0/1/2 的瓦片网格与编号(OSM 风格,y 从上到下):

z=0 (1 张):                  z=1 (4 张):                  z=2 (16 张):
┌────────────┐               ┌─────┬─────┐               ┌──┬──┬──┬──┐
│            │               │ 0,0 │ 1,0 │               │00│10│20│30│
│   (0,0)    │               ├─────┼─────┤               ├──┼──┼──┼──┤
│            │               │ 0,1 │ 1,1 │               │01│11│21│31│
└────────────┘               └─────┴─────┘               ├──┼──┼──┼──┤
                                                          │02│12│22│32│
                                                          ├──┼──┼──┼──┤
                                                          │03│13│23│33│
                                                          └──┴──┴──┴──┘

注意 z=2 的瓦片 (1, 2) 在 z=1 对应父瓦片 (0, 1)(右下子瓦片)。这就是四叉树父子关系的可视化。

常见误区

WARNING

误区 1:以为瓦片是按屏幕像素切的。 ✅ 正确:瓦片按地理范围切(256 像素是输出大小,不是切分依据)。同一瓦片在不同缩放级别下地理范围不同。

WARNING

误区 2:以为 LOD 只看缩放级别(滚轮滚了几次)。 ✅ 正确:3D 地球的 LOD 看相机到地球的距离——离得远用粗 LOD,离得近用细 LOD。同一缩放级别下,旋转视角会让某些瓦片切到不同 LOD。

WARNING

误区 3:以为子瓦片没加载时屏幕是空白。 ✅ 正确:引擎用父瓦片兜底——子瓦片下载时,先拉伸父瓦片填充;子瓦片就绪后淡入(alpha 渐变)。这就是 Google Maps 放大时”先模糊后清晰”的视觉原理。

WARNING

误区 4:以为 TMS 和 OSM 是同一套编号。 ✅ 正确yy 轴方向相反——TMS 从南极开始(下到上),OSM 从北极开始(上到下)。混用会让地图上下颠倒

WARNING

误区 5:以为 Web 墨卡托瓦片和 EPSG:4326 瓦片只是 y 轴方向不同。 ✅ 正确:根本差异是投影方式不同——Web 墨卡托瓦片是投影后的米切分(瓦片是正方形),EPSG:4326 瓦片是直接经纬度切分(瓦片是矩形,纬度越高越窄)。

WARNING

误区 6:以为视锥剔除就够了,不需要地平线剔除。 ✅ 正确:3D 地球视图中,视锥内仍有 30-50% 的瓦片在背面(被地球本体遮挡)。地平线剔除可显著减负。2D 地图不需要(无背面概念)。

WARNING

误区 7:以为 SSE 公式中的 ϵ\epsilon(几何误差)是写死的常量。 ✅ 正确ϵ\epsilon 是每个瓦片元数据自带的字段(3D Tiles 标准的 geometricError),描述”该瓦片离完美细节还差多少米”。城市建模 ϵ10\epsilon \approx 10 m,地形 ϵ100\epsilon \approx 100 m,全球底图 ϵ1000\epsilon \approx 1000 m。

延伸阅读与自测

权威参考

自测题

  1. z/x/y 手算:给定上海东方明珠 (121.50°E, 31.24°N),手算 z=12z=12 时的 Web 墨卡托瓦片编号。验证应得到约 (12, 3369, 1546)
  2. TMS/OSM 转换:某瓦片在 OSM 编号是 (z=10,x=842,y=388)(z=10, x=842, y=388),TMS 编号下 yy 是多少?(应得 2101388=6352^{10} - 1 - 388 = 635
  3. SSE 阈值推导:给定 ϵ=100\epsilon = 100 m,H=1080H = 1080 px,θ=60°\theta = 60°f=1.0f = 1.0d=106d = 10^6 m(1000 km 高空),算 SSE。如果阈值是 16 px,要不要细化到下一层?(应得 SSE ≈ 9.4 px,不需要细化)
  4. 3D vs 2D 设计取舍:2D 地图为什么不用 SSE 而用距离阈值?(提示:2D 相机正交投影,没有 perspective foreshortening)
  5. 瓦片过渡:放大地图时,为什么”先模糊后清晰”而不是”先空白后清晰”?引擎如何实现父子过渡?给出 2 个理由。

下一篇导引04 影像图层 将打开瓦片的内容黑盒——瓦片里装的是什么数据?XYZ / WMS / WMTS 协议差异是什么?矢量图层与栅格图层的本质区别?多图层纹理如何在 GPU 上混合?