Android 火星坐标反向换算无逆公式？二分逼近 GCJ02 转 WGS84

原创于 2026-09-22 07:08:34 发布 · 154 阅读

0 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

收录于

Android 测绘开发实战

话题

#开学季·九月创作之星博客挑战赛

正向变换唾手可得，反向却查不到公式。把正向当黑盒，在区间上二分逼近即可。

前言

做地图、定位、轨迹类应用的人，几乎都绕不开"火星坐标"（GCJ02）与 WGS84 之间的换算。WGS84 是世界通用的经纬度基准，而国内很多地图平台在加密坐标下工作，两者之间不是简单的平移旋转，而是一条带三角级数的非线性映射。

正向换算（WGS84 → GCJ02）网上到处是现成公式，照着抄就行。可一旦你需要在反向（GCJ02 → WGS84）上工作——比如用定位 SDK 拿到火星坐标，再叠加到 WGS84 底图——麻烦就来了：大家都说"反向没有解析公式"。有的方案直接放弃精度，用固定偏移硬凑；有的照搬别人代码，却说不清为什么能收敛。

这篇文章想讲一个实测可行的反向思路：既然正向是现成的、连续的、可迭代的，那就别求逆公式了——把正向变换当成一个黑盒函数，在经纬度区间上做二分逼近。这套引擎在把火星坐标底图叠加回 WGS84 时，正是靠这一招完成反向换算。

一、问题长什么样

假设你在写一个工具方法，输入一个 GCJ02 经纬度，期望输出对应的 WGS84 经纬度。正向你已经有这么一段：

fun fromWgs84(lat: Double, lon: Double): Pair<Double, Double> {
    val a = 6378245.0
    val ee = 0.00669342162296594323
    val dLat = transformLat(lon - 105.0, lat - 35.0)
    val dLon = transformLon(lon - 105.0, lat - 35.0)
    val radLat = lat / 180.0 * PI
    val magic = sin(radLat)
    val sqrtMagic = sqrt(1 - ee * magic * magic)
    val glat = lat + dLat * 180.0 / ((a * (1 - ee)) / (magic * sqrtMagic) * PI)
    val glon = lon + dLon * 180.0 / (a / sqrtMagic * cos(radLat) * PI)
    return glat to glon
}

这里的 transformLat、transformLon 是一堆带 sin、sqrt 的三角级数叠加。问题在于：这个映射是"加密"出来的，设计上就刻意没有给你反函数。你只有一个正向函数，fromWgs84，它把 WGS84 送进去，吐出一对火星坐标。

现在给你一对火星坐标 (gcjLat, gcjLon)，你根本没法从表达式里直接解出 (lat, lon)。

二、为什么不能"直接反推"

首先排除朴素想法：把 fromWgs84 里的公式抄一遍、再硬编码一个相反的偏移。行不通，因为加密函数本身不是线性偏移，而是在经纬度上叠加了一组随位置变化的正弦波，不同位置偏移量不同。你把公式倒过来写，得到的还是一个错的映射，只是看起来像反向。

其次，"逆变换 = 正向变换的逆函数"在解析上不可行。正向函数由多个不可逆的初等函数复合而成，想化简出 lon = f⁻¹(glon) 这样的显式表达式，数学上不现实。

那剩下的路就清晰了：数值求解。把一个"求反函数"的问题，改造成"用二分法求方程根"的问题。

三、解法：把正向函数当黑盒，二分逼近

思路很简单：

给定目标火星坐标 (gcjLat, gcjLon)。
构造一个初始的 WGS84 候选区间：经纬度各自加一个初始步长 delta = 0.01，得到一个粗略范围 [min, max]。
取区间中点 (mLat, mLon)，用 fromWgs84 正向变换，得到它对应的火星坐标。
比较"变换出来的火星坐标"和"目标火星坐标"的偏差：
- 偏差为正，说明当前候选点偏大，把区间上界收缩到中点；
- 偏差为负，说明偏小，把区间下界抬升到中点。
不断把区间对半缩，直到偏差小于阈值 1e-9，或迭代超过 10000 次为止。

核心循环大致是这样：

fun toWgs84(gcjLat: Double, gcjLon: Double): DoubleArray {
    var dLat = 0.01
    var dLon = 0.01
    var mLat = gcjLat - dLat
    var mLon = gcjLon - dLon
    var pLat = gcjLat + dLat
    var pLon = gcjLon + dLon
    var i = 0
    while (true) {
        val wgs = doubleArrayOf((mLat + pLat) / 2, (mLon + pLon) / 2)
        val gcj = fromWgs84(wgs[0], wgs[1])
        dLat = gcj.first - gcjLat
        dLon = gcj.second - gcjLon
        if (abs(dLat) < 1e-9 && abs(dLon) < 1e-9) break
        if (dLat > 0) pLat = wgs[0] else mLat = wgs[0]
        if (dLon > 0) pLon = wgs[1] else mLon = wgs[1]
        if (++i > 10000) break
    }
    return doubleArrayOf((mLat + pLat) / 2, (mLon + pLon) / 2)
}

这个算法之所以能收敛，依赖正向函数的一个关键性质：它关于经纬度是单调递增的。加密映射只产生几十到几百米的偏移，不会"翻转"方向，所以二分搜索总能向目标收敛。这也是为什么阈值可以定到 1e-9（度，约 0.1 毫米级别）还能在有限步内到达。

正向映射单调递增，是二分法能够稳定收敛的数学保证。

四、升华：二分法的通用价值与代价

这套"没有逆函数就用数值逼近"的思路，远不止用于火星坐标：

任何单调可迭代的正向映射，都可以用同样的方法求反。比如某些高程修正、坐标投影的迭代反算，都符合这个模式。
代价是精度与速度的权衡。二分每次只得到半位精度，收敛到 1e-9 需要约 30 轮正向变换，好在每次变换只是几十次三角函数，毫秒级可接受。如果对实时性有更高要求，可改用牛顿迭代（需要导函数）或把结果做小规模缓存。
工程上更稳的取巧：既然偏移量级只有几百米，可以先算一次正向得到一个粗偏差，再直接减去该偏差做一次"预修正"，把初始区间缩得很小，让二分在更少的迭代内收敛。

反过来也提醒我们：正向变换一定要实现成可复用的纯函数（无状态、输入输出确定），否则任何基于它的数值算法都会变得不可预测。