做导航定位相关的开发十有八九会在某个晚上被坐标系问题折腾到怀疑人生。前段时间我帮一个做室内外一体化导航的朋友排查偏差他的车辆轨迹上报用的是WGS84经纬度后端想直接在电子地图上做轨迹匹配结果发现直接拿经纬度当平面坐标用几百米的误差都算小的。问题出在哪就是没搞清楚 WGS84 基准椭球、地图投影方法、以及投影后的平面坐标这三者之间的关系。这篇就围绕WGS84 基准推导 米勒投影算法 C/Matlab 实现展开把这套东西从头到尾捋一遍顺便把我踩过的坑也写进去。适合刚接触地图开发、GNSS数据解析、或者正在做轨迹可视化/路径规划的开发者看完至少能明白经纬度到底怎么变成平面坐标米勒投影的公式是怎么推出来的以及代码里哪些细节不能省。1. 坐标系与基准WGS84 到底是什么1.1 从经纬度到平面坐标的跨越我们在导航设备、手机GPS里拿到的原始定位绝大多数是 WGS84 坐标系下的经纬度对形如(lat31.2304, lon121.4737)。这个经纬度描述的是一个位于椭球面上的点椭球面是什么简单理解地球并非正球体而是一个赤道略鼓、两极略扁的旋转椭球。WGS84 就是一套规定了这个椭球多大、中心在哪、坐标轴怎么指的全球基准。但问题来了人类的活动、地图的绘制、路径规划、渲染引擎几乎全在平面上完成。你不可能拿一个曲面坐标直接去做两点间直线距离计算、去做网格索引、去做屏幕绘制。所以必须有一个数学手段把椭球面上的经纬度摊平到二维平面上这个过程就是地图投影。投影的本质是把三维曲面映射到二维平面同时尽量保持某些属性角度、面积、距离不变形。但任何投影都不可能完全不变形区别只是保什么、牺牲什么。1.2 WGS84 椭球的核心参数WGS84 椭球的定义参数是固定的理解这几个数字对后续投影计算至关重要参数数值说明长半轴 a6378137.0 m赤道半径扁率 f1 / 298.257223563椭球扁平程度短半轴 b约 6356752.3142 m通过 a×(1-f) 计算第一偏心率 e²约 6.69437999014e-3用于大地主题解算在做米勒投影这类圆柱投影时我们往往直接取长半轴 a 作为球体半径 R这是一种工程近似。对于车载导航、WebGIS可视化来说这种近似在绝大多数场景下误差可接受但如果做高精度的工程测量、跨海大桥、长距离管线测绘就不能这么粗放了需要考虑椭球改正。提示很多人把经纬度坐标和投影坐标搞混写代码时分不清谁是输入谁是输出。记住一条铁律WGS84 是输入基准平面坐标是输出结果两者之间必须经过投影变换。1.3 为什么不直接用等距圆柱投影最简单的投影方式就是直接把经纬度当成平面上的 x/yxlonylat。这就是等距圆柱投影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 数学推导与参数自洽性验证我们来验证一下这个公式是否自洽。假设一个点在赤道上φ0y 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.9626ln(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、有人传 lat910这些非法值如果不拦下来后面的计算会输出一堆不可信坐标。#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° 时这个参数就会超过 π/2tan 趋向无穷log 膨胀整个 y 值失控。所以压缩系数不是随便改的。第二反解时的除数是kMillerYScale我写的是除以 1.25。如果你用了y R * log(tan(...)) / 0.8这种写法反解时要乘以 0.8千万别搞混。我建议全项目统一用乘 1.25的写法减少心智负担。3.3 精度与边界测试写完了函数必须做测试这是工程习惯。我拿几个关键坐标点验证输入经纬度预期 x米预期 y米(0°, 0°)00(0°, 90°E)10018754.170(45°N, 0°)0约 5372543.58(80°N, 0°)0约 10111225.63这些预期值是用公式手算得到的。注意 y(45°) ≈ 0.8425Ry(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 坐标大得离谱。比如 lat45如果没转弧度直接代入公式tan(π/40.4×45)里的 45 被解释成弧度数值完全错乱。排查技巧很简单在转换函数入口打印原始值和转换后的弧度值确认弧度在 [-π/2, π/2] 范围内。或者写单元测试固定输入一个已知坐标比如北京的经纬度约 39.9°N, 116.4°E投影后 x 约 1297 万米、y 约 470 万米跟预期差太多就是弧度没转。常见错误表现解决度当作弧度直接代入y 值异常巨大或负值统一用 kD2R 转弧度经度纬度传反x/y 严重错位轨迹镜像接口命名注释必须清楚使用 3.14 近似 π坐标出现系统性偏差用 std::numbers::pi 或 constexpr反解时乘除法搞反闭合测试误差巨大对照正反公式手推一遍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 墨卡托公式后来做极地航线模拟发现纬度一高就绑不住换了米勒投影才算顺眼。回头想想坐标系这个东西最忌讳的就是拿来就用每个投影背后都有一套适用假设搞不清这些假设线上早晚要出幺蛾子。我的建议是代码实现之前先花半天把手推公式过一遍再把正反解、边界、误差都写成测试用例后面整个业务都会稳很多。米勒投影不算复杂但把它吃透了再去看兰伯特、高斯-克吕格这些投影思路都是通的。