three.js地图数学基础(二):Web Mercator

xiaominlaopodaren
122 阅读2分钟

Web Mercator：把经纬网压进单位正方形

矢量瓦片最常见的平面坐标基础是 Web Mercator（在 GIS 行业标准中对应 EPSG:3857，历史上曾使用过 EPSG:900913、EPSG:3785 等标识）。它把经度方向线性映射，把纬度方向按正割函数逐渐拉伸，使局部角度得以保持。代价是面积和距离在高纬度严重放大，而且南北极永远无法映射到有限坐标。

屏幕截图 2026-08-08 224811.png

1. 模型和符号

本书采用球面模型，而不是椭球 Mercator：

  • λ：经度弧度，范围通常为 [−π,π]；
  • φ：纬度弧度；
  • R：球体半径；
  • (X_m,Y_m)：以米或任意同尺度单位表示的未归一化投影；
  • (x,y)：归一化到单位世界的坐标；
  • 屏幕/瓦片约定让 y 向南增大。

角度输入来自用户时通常是度。必须先换为弧度：

θrad=θdegπ180\theta_{rad}=\theta_{deg}\frac{\pi}{180}

JavaScript 的 Math.sin、Math.tan、Math.atan 等函数全部接收弧度。

1.1 球面模型与椭球 Mercator 的区别

在地理信息系统中，Mercator 投影有两种主流实现：

  1. 椭球 Mercator（如 EPSG:3395）： 严格按照 WGS84 椭球体长短半轴公式计算。推导积分时需引入第一偏心率 ee，其纵坐标公式为： Ym=Rln⁡[tan⁡(π4+φ2)⋅(1−esin⁡φ1+esin⁡φ)e/2]Y_m = R \ln \left[ \tan\left(\frac{\pi}{4}+\frac{\varphi}{2}\right) \cdot \left( \frac{1 - e\sin\varphi}{1 + e\sin\varphi} \right)^{e/2} \right]
  2. 球面 Web Mercator（EPSG:3857）： 假定地球为完美正球体（半径取 WGS84 长半轴 R=6,378,137 mR = 6,378,137 \text{ m}），直接采用前述简化积分公式。

为什么 Web Mercator 放弃了更精确的椭球模型？

  • 计算性能极高： 球面公式省去了偏心率修正项，计算量显著降低，极其适合渲染引擎和客户端进行实时坐标变换；
  • 完美正方形切割： 球面模型下，南北边界截断后，整个投影世界是一个严格的 2πR×2πR2\pi R \times 2\pi R 正方形，非常适合用四叉树（Quadtree）进行逐级切片（XYZ 瓦片）；而椭球模型很难直接导出完美的正方形瓦片结构；
  • 定位与测绘的区别： Web 电子地图的核心目标是视觉渲染与交互而非测绘工程。虽然球面模型会导致最大约 20 公里左右的绝对纬度定位偏移（若不经转换直接对比），但由于所有 Web 瓦片与要素都统一使用这套球面规则，因此在地图上叠加矢量数据时，相对位置完全对齐，不存在视觉错位。

2. 经度为什么是线性的

Mercator 圆柱投影把赤道方向展开为直线。经度相差 2π 正好绕球一周，所以最自然的未归一化横坐标为：

Xm=RλX_m=R\lambda

将 λ=-π 平移到 0，再除以完整周长 2π，得到单位世界横坐标：

x=λ+π2π\boxed{x=\frac{\lambda+\pi}{2\pi}}

于是：

  • −180° → x=0；
  • 0° → x=0.5；
  • +180° → x=1。

公式本身不必截断经度。λ>π 或 λ<−π 会自然产生 [0,1] 外的 x，后续由 wrap 记录重复世界。这一点对跨反经线的连续平移很重要。

3. 纬度公式从哪里来

3.1 保角条件

球面上一小段纬线在东西方向的实际长度随纬度缩短。纬度为 φ 时，平行圈半径为 R cosφ，因此经度微分 dλ 对应的地表长度为：

dsx=Rcos⁡φ dλds_x=R\cos\varphi\,d\lambda

南北方向的纬度微分长度为：

dsy=R dφds_y=R\,d\varphi

投影横向微分为 dX_m=R dλ。若希望局部东西、南北使用同一放大比例，纵向投影必须满足：

dYmR dφ=dXmRcos⁡φ dλ=1cos⁡φ\frac{dY_m}{R\,d\varphi} = \frac{dX_m}{R\cos\varphi\,d\lambda} = \frac{1}{\cos\varphi}

因此：

dYm=Rsec⁡φ dφdY_m=R\sec\varphi\,d\varphi

对纬度积分：

Ym=R∫sec⁡φ dφ=Rln⁡∣sec⁡φ+tan⁡φ∣+CY_m=R\int\sec\varphi\,d\varphi =R\ln|\sec\varphi+\tan\varphi|+C

取赤道 φ=0 时 Y_m=0，可得 C=0。利用恒等式：

sec⁡φ+tan⁡φ=tan⁡(π4+φ2)\sec\varphi+\tan\varphi =\tan\left(\frac{\pi}{4}+\frac{\varphi}{2}\right)

得到经典球面 Mercator 纵坐标：

Ym=Rln⁡tan⁡(π4+φ2)Y_m=R\ln\tan\left(\frac{\pi}{4}+\frac{\varphi}{2}\right)

3.2 从向北为正改为向南为正

数学投影的 Y_m 向北增大；XYZ 瓦片和屏幕坐标却从北向南编号。完整世界在横向的跨度为 2πR，将赤道放在 y=1/2 并翻转轴：

y=12−ln⁡tan⁡(π4+φ2)2π\boxed{ y=\frac{1}{2} -\frac{\ln\tan\left(\frac{\pi}{4}+\frac{\varphi}{2}\right)}{2\pi} }

赤道上 tan(π/4)=1，ln1=0，所以 y=0.5。

在代码实现（如 JavaScript）中，利用双曲函数恒等式 ln⁡(sec⁡φ+tan⁡φ)=atanh(sin⁡φ)\ln(\sec\varphi + \tan\varphi) = \text{atanh}(\sin\varphi)，该公式可进一步表达为：

y=12−atanh(sin⁡φ)2π\boxed{ y = \frac{1}{2} - \frac{\text{atanh}(\sin\varphi)}{2\pi} }

工程优势： 使用 Math.atanh(Math.sin(phi)) 相比传统的 Math.log(Math.tan(...))，减少了一次加法、一次除法与多次三角函数嵌套，在计算机浮点运算中拥有更好的指令性能和数值稳定性。

4. 为什么必须截断纬度

当 φ→π/2 时，tan(π/4+φ/2)→∞，投影纵坐标也趋向无穷。南极同理。因此不可能把两个极点放进有限正方形。

XYZ 瓦片选择的截止纬度，使投影纵坐标刚好达到横向半周，即 normalized y=0 或 1。令北边界 y=0：

ln⁡tan⁡(π4+φmax2)=π\ln\tan\left(\frac{\pi}{4}+\frac{\varphi_{max}}{2}\right)=\pi

由反函数（下一节证明）可得：

φmax=arctan⁡(sinh⁡π)\varphi_{max}=\arctan(\sinh\pi)

换成角度：

φmax=85.0511287798066∘\boxed{\varphi_{max}=85.0511287798066^\circ}

任何更北或更南的输入都必须先截断到 ±φmax，否则会得到超出有限世界的数值，靠近极点时还可能出现 Infinity。

截断不是“修复错误坐标”，而是这套平面瓦片模型的定义域。若产品必须展示极点，应选择地球/极区投影，而不是继续扩大 Web Mercator。

工程截断时机： 在代码实现中，截断必须在传入任何三角函数计算之前（或紧随度转弧度之后）完成：

// 正确：先做范围约束，再进入三角函数，防止 ±90° 计算抛出 Infinity/NaN
const clampedLat = Math.max(-85.0511287798066, Math.min(85.0511287798066, lat));
const phi = (clampedLat * Math.PI) / 180;

5. 逆投影推导

横坐标直接变形：

λ=2πx−π\boxed{\lambda=2\pi x-\pi}

对于纬度，令：

a=π(1−2y)a=\pi(1-2y)

由正向公式可知：

a=ln⁡tan⁡(π4+φ2)a=\ln\tan\left(\frac{\pi}{4}+\frac{\varphi}{2}\right)

指数化后：

ea=tan⁡(π4+φ2)e^a=\tan\left(\frac{\pi}{4}+\frac{\varphi}{2}\right)

可写成常见反式：

φ=2arctan⁡(ea)−π2\varphi=2\arctan(e^a)-\frac{\pi}{2}

利用 atan(sinh a) 与上式等价，得到数值实现更直接的形式：

φ=arctan⁡(sinh⁡(π(1−2y)))\boxed{\varphi=\arctan\left(\sinh\left(\pi(1-2y)\right)\right)}

最后把 λ/φ 从弧度乘 180/π 转回度。

6. 两组可手算的检查

6.1 原点检查

输入 (0°,0°)：

x=0+π2π=0.5x=\frac{0+\pi}{2\pi}=0.5
y=12−ln⁡tan⁡(π/4)2π=12−ln⁡12π=0.5y=\frac12-\frac{\ln\tan(\pi/4)}{2\pi} =\frac12-\frac{\ln1}{2\pi}=0.5

所以结果必须精确为 (0.5,0.5)。若不是，最常见原因是把角度直接传给三角函数，或忘记 y 轴翻转。

6.2 上海检查

输入：

(λdeg,φdeg)=(121.4737,31.2304)(\lambda_{deg},\varphi_{deg})=(121.4737,31.2304)

投影结果约为：

(x,y)=(0.8374269444444444,  0.40860372471489265)(x,y)=(0.8374269444444444,\;0.40860372471489265)

将这两个值送入逆投影，应恢复：

(121.4737∘,31.2304∘)(121.4737^\circ,31.2304^\circ)

双精度往返误差应远小于 10⁻⁹°。这里检验的是公式与实现的一致性，不代表地面测量精度达到纳米级。

7. 失真如何影响瓦片引擎

Mercator 的局部线性放大比例为：

k=sec⁡φ=1cos⁡φk=\sec\varphi=\frac{1}{\cos\varphi}

赤道 k=1；纬度 60° 时 k=2。这意味着投影平面上一段固定长度对应的真实地表距离随纬度变化。它会影响：

  • “每像素多少米”的计算；
  • 以米定义的建筑高度或线宽；
  • 相机高度与缩放级别之间的换算；
  • 距离标尺。

但瓦片寻址仍在 normalized Mercator 中进行。不要为了修正米制失真而改变 x/y 或 XYZ 规则；应该在需要物理尺度的模块中显式加入纬度比例。

7.1 归一化坐标向像素/瓦片坐标的映射

归一化坐标 (x,y)∈[0,1](x,y) \in [0, 1] 是跨缩放级别的统一表达。在特定缩放级别（Zoom Level ZZ）及瓦片像素尺寸 SS（通常为 256px 或 512px）下，坐标换算公式为：

  • 世界像素坐标： px=x⋅S⋅2Zpx = x \cdot S \cdot 2^Z py=y⋅S⋅2Zpy = y \cdot S \cdot 2^Z

  • 瓦片索引 (Tile X / Tile Y)： TileX=⌊x⋅2Z⌋TileX = \lfloor x \cdot 2^Z \rfloor TileY=⌊y⋅2Z⌋TileY = \lfloor y \cdot 2^Z \rfloor

8. 实现契约

数学模块应满足：

  1. projectLngLat 接收度，输出无量纲 normalized Mercator；
  2. 纬度在三角函数之前截断；
  3. 经度不自动折返，以保留连续世界信息；
  4. unprojectMercator 接收 normalized 数值，输出度；
  5. 模块不依赖 DOM、Three.js、瓦片缓存或当前相机；
  6. 返回值不可变，避免检查器或渲染代码意外改写共享结果。

下一篇教程将把这些公式写成一个依赖为零的 ES 模块，并用原点、上海和边界输入逐项核对。

avatar
文章
阅读
粉丝