做导航定位相关的开发，十有八九会在某个晚上被坐标系问题折腾到怀疑人生。前段时间我帮一个做室内外一体化导航的朋友排查偏差，他的车辆轨迹上报用的是WGS84经纬度，后端想直接在电子地图上做轨迹匹配，结果发现直接拿经纬度当平面坐标用，几百米的误差都算小的。问题出在哪？就是没搞清楚 WGS84 基准椭球、地图投影方法、以及投影后的平面坐标这三者之间的关系。这篇就围绕"WGS84 基准推导 + 米勒投影算法 + C++/Matlab 实现"展开，把这套东西从头到尾捋一遍，顺便把我踩过的坑也写进去。适合刚接触地图开发、GNSS数据解析、或者正在做轨迹可视化/路径规划的开发者，看完至少能明白：经纬度到底怎么变成平面坐标，米勒投影的公式是怎么推出来的，以及代码里哪些细节不能省。

1. 坐标系与基准：WGS84 到底是什么

1.1 从经纬度到平面坐标的跨越

我们在导航设备、手机GPS里拿到的原始定位，绝大多数是 WGS84 坐标系下的经纬度对，形如(lat=31.2304, lon=121.4737)。这个经纬度描述的是一个位于椭球面上的点，椭球面是什么？简单理解，地球并非正球体，而是一个赤道略鼓、两极略扁的旋转椭球。WGS84 就是一套规定了"这个椭球多大、中心在哪、坐标轴怎么指"的全球基准。

但问题来了：人类的活动、地图的绘制、路径规划、渲染引擎，几乎全在"平面"上完成。你不可能拿一个曲面坐标直接去做两点间直线距离计算、去做网格索引、去做屏幕绘制。所以必须有一个数学手段，把椭球面上的经纬度"摊平"到二维平面上，这个过程就是地图投影。投影的本质是：把三维曲面映射到二维平面，同时尽量保持某些属性（角度、面积、距离）不变形。但任何投影都不可能完全不变形，区别只是"保什么、牺牲什么"。

1.2 WGS84 椭球的核心参数

WGS84 椭球的定义参数是固定的，理解这几个数字对后续投影计算至关重要：

在做米勒投影这类圆柱投影时，我们往往直接取长半轴 a 作为球体半径 R，这是一种工程近似。对于车载导航、WebGIS可视化来说，这种近似在绝大多数场景下误差可接受；但如果做高精度的工程测量、跨海大桥、长距离管线测绘，就不能这么粗放了，需要考虑椭球改正。

提示：很多人把"经纬度坐标"和"投影坐标"搞混，写代码时分不清谁是输入谁是输出。记住一条铁律：WGS84 是输入基准，平面坐标是输出结果，两者之间必须经过投影变换。

1.3 为什么不直接用等距圆柱投影

最简单的投影方式，就是直接把经纬度当成平面上的 x/y：x=lon，y=lat。这就是等距圆柱投影（Equirectangular），它计算最快，但变形巨大。在赤道附近凑合能用，一旦跑到中高纬度，纬线方向（东西）的拉伸和经线方向（南北）的拉伸严重不一致，一个圆形区域会被拉成椭圆，距离、面积、角度全都失真。导航里经常要计算行驶距离、判断转弯角度，这种投影根本扛不住。

所以行业内最常用的方案有两种：一种是保角的墨卡托投影及其变体（Web 墨卡托），另一种是保面积的各类等积投影。而米勒投影（Miller Cylindrical Projection）属于圆柱投影家族中偏向"视觉友好型"的一种，它比墨卡托在高纬度地区的变形小，又比等距圆柱投影的纬度分布更合理，因此在一些世界地图渲染、低精度导航态势展示里，米勒投影是个不错的中间选择。

2. 米勒投影：原理、公式与设计逻辑

2.1 为什么在墨卡托和等距圆柱之间折中

标准墨卡托投影的纵坐标公式是：

y = R × ln(tan(π/4 + φ/2))

其中 φ 是纬度（弧度）。这个公式的特点很鲜明：赤道附近 y 随 φ 线性增长，越靠近两极 y 增长越快。当 φ 接近 ±90° 时，y 趋向无穷大，这就是为什么 Web 墨卡托要把纬度截断在 ±85.05112878° 左右。

墨卡托投影的优点是恒向线（等角航线）是直线，船在海上的航行规划特别方便；但它有个大问题：高纬度地区面积被严重夸大。格陵兰岛在墨卡托地图上看起来跟非洲差不多大，实际只有非洲的 1/14。

米勒投影的思路很巧妙：既然 φ 接近 ±90° 时 tan() 项会爆炸，那我能不能让"有效纬度"增长的慢一点？1942 年 Miller 提出的做法是：先对纬度做一个压扁变换，再套墨卡托公式。公式如下：

x = R × λ
y = R × 1.25 × ln(tan(π/4 + 0.4 × φ))

或者写成更常见的形式：

y = R × ln(tan(π/4 + 0.4×φ)) / 0.8

注意 1.25 和 0.8 互为倒数。这里的 0.4 就是压缩系数，把纬度乘 0.4，相当于喂给 tan() 的"输入角度"永远不超过 0.4×90° = 36°，也就是 π/4 + 0.4φ 永远不会超过 π/4 + 0.4×π/2 ≈ 1.727 弧度，约 98.96°，tan() 的结果不会爆炸。最后再乘 1.25 把纵向整体拉伸回来，使赤道附近的比例尺和等距圆柱投影保持一致。

2.2 数学推导与参数自洽性验证

我们来验证一下这个公式是否自洽。假设一个点在赤道上，φ=0：

y = R × 1.25 × ln(tan(π/4 + 0)) = R × 1.25 × ln(1) = 0

正确，赤道上的纵坐标是 0。

假设 φ = 45°（π/4 弧度），代入：

y = R × 1.25 × ln(tan(π/4 + 0.4×π/4)) = R × 1.25 × ln(tan(π/4 + π/10)) = R × 1.25 × ln(tan63°)

tan63° ≈ 1.9626，ln(1.9626) ≈ 0.674，乘以 1.25 ≈ 0.8425R，即约 5371 km。作为对照，标准墨卡托在 45° 纬度处的 y 值是 R×ln(tan(45°+22.5°)) ≈ R×ln(tan67.5°) ≈ 0.8814R ≈ 5618 km。可以看出米勒投影在 45° 纬度的纵向拉伸确实比墨卡托小一些。

再看反解公式。已知 y，求 φ：

y / (1.25R) = ln(tan(π/4 + 0.4φ)) exp(y / (1.25R)) = tan(π/4 + 0.4φ) 0.4φ = arctan(exp(y / (1.25R))) - π/4 φ = (arctan(exp(y / (1.25R))) - π/4) / 0.4

这是完全可逆的代数变换，正反解闭合成环。这一点在工程里很重要：如果你要把平面坐标再转回经纬度（比如从地图上点击获取坐标），反解公式必须简洁稳定。

2.3 米勒投影和 Web 墨卡托的适用边界

很多人问我："现在 WebGIS 不都用 Web 墨卡托吗？为什么还要折腾米勒投影？"我的理解是这样：

• Web 墨卡托是互联网地图的事实标准，瓦片切割、地图缩放层级全是按它设计的，如果做在线地图叠加，直接用 Web 墨卡托，省去一切兼容问题。
• 但 Web 墨卡托的纬度截断是硬伤。真要在极地附近部署传感器网络、做航线模拟，Web 墨卡托的 y 值会非常夸张，超出瓦片范围。
• 米勒投影在高纬度的 y 值增长比墨卡托平缓，视觉上更接近"真实形状"，所以在一些自建可视化和低精度导航系统里，米勒投影反而更合适。

实际项目中，如果只是做市区范围内的轨迹展示，米勒投影和 Web 墨卡托的差异很小。一旦跨度达到数百公里、涉及高纬度，就得掂量掂量选哪一种了。

3. C++ 实现：工程代码与边界处理

3.1 头文件与核心类设计

C++ 实现的重点不是公式本身——公式就两行的事——重点是类的设计、参数校验和复用性。我习惯把坐标转换封装成命名空间，输入输出用结构体，避免魔法数字散落各处。

#ifndef GEO_CONV_H #define GEO_CONV_H #include <cmath> #include <stdexcept> namespace geo_conv { constexpr double kPi = 3.14159265358979323846; constexpr double kD2R = kPi / 180.0; constexpr double kR2D = 180.0 / kPi; constexpr double kA = 6378137.0; // WGS84 长半轴 constexpr double kF = 1.0 / 298.257223563; // WGS84 扁率 constexpr double kMillerLatScale = 0.4; // Miller 投影纬度压缩系数 constexpr double kMillerYScale = 1.25; // 纵向恢复系数 struct PointLL { double lat; // 纬度，单位：度 double lon; // 经度，单位：度 }; struct PointXY { double x; // 米 double y; // 米 }; PointXY wgs84ToMiller(const PointLL& ll); PointLL millerToWgs84(const PointXY& xy); } // namespace geo_conv #endif

这里把 kPi、kD2R 这些作为编译期常量，好处是编译器能内联展开，运行期零开销。你如果在项目里看到有人用#define PI 3.14，直接告诉他：double 精度不够，后面投影误差会让你哭。3.14 和真正的 π 在长距离计算里能差出好几公里。

3.2 转换函数实现

接下来是核心转换函数，我把正反解都写出来，并且加上了参数合法性校验。这一步很容易被忽略，但线上数据千奇百怪——有人给你传 NaN、有人传 lat=910，这些非法值如果不拦下来，后面的计算会输出一堆不可信坐标。

#include "geo_conv.h" namespace geo_conv { PointXY wgs84ToMiller(const PointLL& ll) { if (!std::isfinite(ll.lat) || !std::isfinite(ll.lon)) { throw std::invalid_argument("经纬度不是有限数值"); } if (ll.lat < -90.0 || ll.lat > 90.0 || ll.lon < -180.0 || ll.lon > 180.0) { throw std::invalid_argument("经纬度超出合法范围"); } const double phi = ll.lat * kD2R; const double lambda = ll.lon * kD2R; const double x = kA * lambda; const double y = kA * kMillerYScale * std::log(std::tan(kPi / 4.0 + kMillerLatScale * phi)); return {x, y}; } PointLL millerToWgs84(const PointXY& xy) { if (!std::isfinite(xy.x) || !std::isfinite(xy.y)) { throw std::invalid_argument("平面坐标不是有限数值"); } if (std::abs(xy.x) > kA * kPi) { throw std::invalid_argument("经度方向超出投影范围"); } const double phi = (std::atan(std::exp(xy.y / (kA * kMillerYScale))) - kPi / 4.0) / kMillerLatScale; const double lambda = xy.x / kA; return {phi * kR2D, lambda * kR2D}; } } // namespace geo_conv

实现的时候有两个细节容易出错：

第一，std::log(std::tan(...))中 tan 的参数kPi / 4.0 + kMillerLatScale * phi在正常经纬度范围内恒为正，开方和 log 不存在定义域问题，所以不需要额外判断。但要注意，如果哪一天你把 kMillerLatScale 改大了，比如改成 0.5，那么纬度 89° 时这个参数就会超过 π/2，tan 趋向无穷，log 膨胀，整个 y 值失控。所以压缩系数不是随便改的。

第二，反解时的除数是kMillerYScale，我写的是除以 1.25。如果你用了y = R * log(tan(...)) / 0.8这种写法，反解时要乘以 0.8，千万别搞混。我建议全项目统一用"乘 1.25"的写法，减少心智负担。

3.3 精度与边界测试

写完了函数必须做测试，这是工程习惯。我拿几个关键坐标点验证：

这些预期值是用公式手算得到的。注意 y(45°) ≈ 0.8425R，y(80°) 已经接近 10111 km。作为对比，Web 墨卡托在 80° 纬度的 y 值是 R×ln(tan(45°+40°)) ≈ R×ln(tan85°) ≈ 6378137×2.436 ≈ 15537 km，明显比米勒投影大得多。这说明米勒投影确实把高纬度的"膨胀"压下来了。

另外建议做一次"往返闭合测试"：随机生成几万个经纬度点，正转再反转，看误差是否在毫米级别。理论上由于双精度浮点运算，往返误差应在 1e-6 米量级。如果发现误差达到米级，多半是单位转换写错了（比如漏了 kD2R）。

4. Matlab 实现：快速验证与可视化

4.1 向量化实现

做原型验证的时候，我一般先用 Matlab 把算法跑通，确认公式和边界行为没问题，再写 C++ 工程代码。Matlab 的优势是向量化计算和绘图方便，几百行 C++ 才能做完的可视化，Matlab 几行就出来了。

function [x, y] = wgs84_to_miller(lat, lon) % WGS84 经纬度转米勒投影平面坐标 % 输入 lat, lon 可为标量或等尺寸矩阵，单位：度 % 输出 x, y 单位：米 R = 6378137.0; phi = lat * pi / 180; lambda = lon * pi / 180; x = R * lambda; y = R * 1.25 * log(tan(pi/4 + 0.4 * phi)); end

反向函数：

function [lat, lon] = miller_to_wgs84(x, y) % 米勒投影平面坐标转 WGS84 经纬度 % 输入 x, y 可为标量或等尺寸矩阵，单位：米 % 输出 lat, lon 单位：度 R = 6378137.0; phi = (atan(exp(y / (R * 1.25))) - pi/4) / 0.4; lambda = x / R; lat = phi * 180 / pi; lon = lambda * 180 / pi; end

向量化实现的要点是：lat、lon 如果是矩阵，整个运算过程矩阵化，不要 for 循环。上述代码天然支持矩阵输入，因为 log、tan、exp 都是逐元素运算。你在脚本里传一个 100×100 的网格坐标，一次调用就能得出全部平面坐标。

4.2 用真实轨迹数据做可视化验证

为了验证投影效果，我建议用真实 GNSS 轨迹做一次可视化。可以从自己的设备导出一条轨迹，或者用公开数据集。假设你有一段轨迹数据存在trajectory.csv里，包含 lat 和 lon 两列，读取并投影的代码：

data = readmatrix('trajectory.csv'); lat = data(:, 1); lon = data(:, 2); [x, y] = wgs84_to_miller(lat, lon); figure; subplot(1, 2, 1); plot(lon, lat, '.-'); xlabel('lon (deg)'); ylabel('lat (deg)'); title('原始经纬度轨迹'); axis equal; subplot(1, 2, 2); plot(x, y, '.-'); xlabel('x (m)'); ylabel('y (m)'); title('米勒投影轨迹'); axis equal;

这里有个细节很多人不注意：第一个图里axis equal其实是"假的等比例"，因为经纬度单位不同也没法真等比例；第二个图里axis equal才有真实意义——x、y 单位都是米，等比例之后你才能直观看出轨迹的实际形状。如果发现轨迹在投影图上形状和地图上明显不一致，通常是你把经度当成了 x、纬度当成了 y，或者单位没转弧度。

4.3 生成全球网格做变形分析

还有一种测试方法，就是生成全球经纬度网格，投影到平面上看网格形态。这个能帮你建立对投影变形的直觉：

[lat_grid, lon_grid] = meshgrid(-80:10:80, -180:20:180); [x_grid, y_grid] = wgs84_to_miller(lat_grid, lon_grid); figure; plot(x_grid, y_grid, 'k'); hold on; plot(x_grid', y_grid', 'k'); axis equal; grid on; xlabel('x (m)'); ylabel('y (m)'); title('Miller 投影全球网格');

运行之后你会看到：经线是一组竖直等距直线，纬线是一组水平曲线，越往两极间距越大，但不会像墨卡托那样在 85° 附近突然炸开。这个网格图拿去做汇报或者写技术文档，比纯文字说明直观得多。

5. 常见问题排查与工程技巧

5.1 度/弧度混用是最隐蔽的 bug

我见过不下十次这种问题：代码里公式写对了，但调用时忘了把度转弧度，结果 y 坐标大得离谱。比如 lat=45，如果没转弧度直接代入公式，tan(π/4+0.4×45)里的 45 被解释成弧度，数值完全错乱。

排查技巧很简单：在转换函数入口打印原始值和转换后的弧度值，确认弧度在 [-π/2, π/2] 范围内。或者写单元测试，固定输入一个已知坐标（比如北京的经纬度约 39.9°N, 116.4°E），投影后 x 约 1297 万米、y 约 470 万米，跟预期差太多就是弧度没转。

5.2 投影变形不是 bug，是特性

很多新手拿到投影坐标，第一反应是"怎么距离和地图上量的对不上？"我举个例子：在北纬 60° 附近，同一经度间隔 1° 对应的地面距离大约是 111.32km × cos(60°) ≈ 55.7km，但在米勒投影的 x 方向上，1° 经度对应的平面距离始终是 R×(1°×π/180) ≈ 111.3km。所以你在投影图上的"东西距离"被拉长了约一倍。这不是公式错误，而是等角圆柱投影的固有特性。

如果你需要在投影平面上做距离计算，必须根据当前纬度乘以 cos(lat) 修正。我自己在项目里常用的做法是：小范围（半径 10km 以内）直接用米勒平面坐标做欧氏距离，误差在可接受范围；大范围或者需要精确距离的场合，改用测地线距离算法（如 Vincenty 公式），不要跟投影变形硬刚。

5.3 椭球近似引起的误差边界

之前提到米勒投影默认把地球当作球体，半径取 WGS84 长半轴 6378137m。那这种近似到底误差多大？我做过一个估算：在纬度 60° 处，WGS84 椭球的纬线半径约为 a×cos(60°)/sqrt(1-e²sin²60°)，和球体半径 a 相比，差异主要体现在 y 方向的子午线弧长上。横轴 x 本来就只与经度差和纬度有关，直接用球体近似误差很小；纵轴 y 在纬度 60° 处，球面近似和椭球面的子午线弧长差异大约 0.2%~0.3%，换算到 1000km 的跨度上就是 2~3km，这个量级做海洋导航绝对不行，做城市道路可视化则可以忽略。

如果你的业务场景对纵轴精度有硬指标，需要先把大地纬度转换成等量纬度（authalic latitude）或归化纬度（reduced latitude），再喂给米勒公式。

5.4 工程落地时的性能考量

最后聊点工程上的性能问题。导航设备通常每秒上报 1~10 个坐标点，单次转换计算量微乎其微，纯 CPU 计算完全不是瓶颈。真正要关注的是：

• 如果每秒处理百万级点（比如做历史轨迹回放、网格聚合统计），可以用查表法：把纬度分成 18000 份（每 0.01° 一个采样点），预计算 y 值表，运行时线性插值。实测可以把单次转换时间压到纳秒级。
• 如果用 GPU 做批量投影（比如 WebGL 可视化），尽量把转换放在着色器里做，把 R、kPi 这些常量作为 uniform 传入，避免 CPU-GPU 往返拷贝。
• 多线程环境下，转换函数是纯函数，无共享状态，可以放心并行。唯一要注意的是异常处理：如果你在 millions 级循环里抛异常，内存开销会非常难看。我在生产代码里通常会提供一个不抛异常的版本，用返回值或输出参数表示失败。

结尾：一点个人体会

这套投影转换我前前后后用过好几种方案：最先是项目里直接抄 Web 墨卡托公式，后来做极地航线模拟发现纬度一高就绑不住，换了米勒投影才算顺眼。回头想想，坐标系这个东西最忌讳的就是"拿来就用"，每个投影背后都有一套适用假设，搞不清这些假设，线上早晚要出幺蛾子。我的建议是：代码实现之前，先花半天把手推公式过一遍，再把正反解、边界、误差都写成测试用例，后面整个业务都会稳很多。米勒投影不算复杂，但把它吃透了，再去看兰伯特、高斯-克吕格这些投影，思路都是通的。