关于流传的 WGS-84 至 GCJ-02 转换算法


前言


在夏天去西藏之前,我曾仔细查询了前往西藏的路线、海拔等等资料。由于那些地方的地图稀缺,很多时候我不得不记录下关键路口、山口等地的经纬坐标。从 Google Earth 上查询卫星图和高程图也很有用。

在这个过程中,我发现 GCJ-02 坐标给我带来了很多麻烦。于是写此文与大家分享。

什么是 WGS-84、GCJ-02?


两个坐标都是以经度、纬度对来标识地球上任意一点的坐标系统。

WGS-84:

WGS 全称 World Geodetic System,于1984年建立。

通过 GPS 获得的经纬坐标,以及 Google Maps 上记录各地形、街道、建筑所使用的经纬坐标,都是 WGS-84 坐标系中的坐标。WGS-84 亦是最通用的地球坐标系。其起初就是因 GPS 的诞生而被设立。

GCJ-02

某毒百科声称:
GCJ-02是由中国国家测绘局制订的地理信息系统的坐标系统。
它是一种对经纬度数据的加密算法,即加入随机的偏差。
国内出版的各种地图系统(包括电子形式),必须至少采用GCJ-02对地理位置进行首次加密。

Wikipedia 亦有描述:
中国官方要求所有在中国运行的地图服务商要加装“国家保密插件”(亦称加密插件、加偏或SM模组),以“保障国家安全”。此插件会将真实的坐标加密成虚假的坐标,且此加偏并非线性加偏,所以各地的偏移情况都会有所不同。

事实上,GCJ-02 意味“国测局-2002”,也就是说,这是国家测绘局于2002年弄出的标准。

事实上,百毒等在线地图为了配合政府,以及保护自己的商业利益,在GCJ-02的基础上都使用了进一步的其它坐标加密算法。当然,那些不属于此篇文章的讨论范围之内。

Case study: Google 地图


上面说过,很多国内地图提供商都龌龊地进行了二次加密。但 Google 即使不得不配合中国政府,也只使用了一次加密,也就是 WGS-84 至 GCJ-02。

事实上,Google 的地图服务有两个版本:
一为 Google Maps,可以通过 maps.google.com 访问
二为 Google Ditu,可以通过 ditu.google.cnditu.google.commaps.google.cn 访问
显而易见,第二个版本是为了中国准备的。

可以注意到,在 Google Ditu 下,卫星图和公路图是重合的,没有任何问题;然而在 Google Maps 下,看中国内的区域,在同一个地方从公路图切换到卫星图时,可以发现发生了较大偏移。借助 Google 官方的经纬度标注工具,用下面这张图加以说明:


很容易看出,在 Google Ditu 上,公路图是在 GCJ-02 坐标系中的,为了让卫星图与之重合,Google 把卫星图也进行了同样的加偏。而在 Google Maps 中,虽然公路图依然是在 GCJ-02 坐标系中的(测绘阶段受中国政府限制),然而卫星图却是未经过加偏的。

通过取了很多点的两组不同坐标,进行偏移计算并与之前流传的算法进行比对,证明之前流传的算法是完全属实的

对于 GCJ-02 的评论


GCJ-02 与 WGS-84 间的转换实际上很容易做到,因为算法已经流传甚广。然而,这样的转换造成了大量的麻烦。

有人曾经评论过
国家的保密插件,是需要收费的,早期的时候,一个导航仪就需要10块钱的保密插件许可费。

至于这个插件的影响,使用过iPhone(iOS4.2以下版本)以及Android手机的用户应该体验深刻,用户在iPhone版的谷歌地图里打开“我的位置”,会发现自己的定位被偏移了几百米,在外出找路的时候尤其麻烦,让用户根本就无法通过谷歌地图来确定方位,以至于迷路。不明真相的用户纷纷谴责谷歌地图的定位准确度不好,其实真的是冤枉了谷歌。

好在这个所谓的“国家的保密插件”并不难破解,网友很容易就可以将其解密,把虚假的坐标转换成真实的坐标,这样,iPhone用户只需安装一个“中国地图校正”插件,就可以自动解决地图偏移的问题。对于Android用户来说,也有修改版的谷歌地图可以下载使用。

总之,国家测绘局的这个所谓的发明创新,最大的用处就是收钱和折腾用户,至于其实际保密效果就不敢恭维了。

WGS-84 与 GCJ-02间的转换


两个坐标系都要能够唯一地标识地球上的每一个点。因而 WGS-84 和 GCJ-02 间的转换必然是一一对应的。一个 WGS-84 坐标便有一个对应的 GCJ-02 坐标,反之亦然。也就是说,他们相互之间的转换讲的也就是同一个点在两个坐标系里坐标的转换而已,没有什么特别的。

也就是说,转换的输入为一对 WGS-84 下的经纬度对,输出为 GCJ-02 下同一点的经纬度对。

我所能找到的最早流出的 WGS-84 至 GCJ-02 的转换算法出自 https://on4wp7.codeplex.com/SourceControl/changeset/view/21483。我们看一下它的代码中的一小部分,处理纬度加偏的前半部分:
-100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.Sqrt(Math.Abs(x)) + (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0 + (20.0 * Math.Sin(y * pi) + 40.0 * Math.Sin(y / 3.0 * pi)) * 2.0 / 3.0 + (160.0 * Math.Sin(y / 12.0 * pi) + 320 * Math.Sin(y * pi / 30.0)) * 2.0 / 3.0

这一部分已经体现了 GCJ-02 设计精巧的几点:

这是一个很长的多项式,并且是很复杂的多项超越式。于是它难以进行反向运算、难以根据已知数据点反推出计算公式。且就算我们拿到了这个已知的转换算法的公式,要想求出它的反函数是非常困难的,于是如果我们要完成 GCJ-02 到 WGS-84 的逆转换,只能通过二分法来逼近。

这个转换算法需要达到的要求是显而易见的:它必须是连续的,并且是单调的。连续才能确保地图不出现断点且覆盖所有区域;单调才能保证原先的位置相互关系转换后依旧成立。

进一步分析


既然转换太复杂,我实在是连求个导的心情都没有,我们就进行一些数值计算看看吧。

若以墨卡托投影来看,进行计算的范围在一个矩形内:
最西为东经 72°,最东为东经 135°。最男为北纬 18°,最北为北纬 54°。这个矩形恰好能包含除了南海诸岛外的中国领土。

通过以很小的步长进行穷算,并且再进一步进行逼近,我们发现:

在这个矩形内,纬度偏差最大处在 26.71°N, 72.60°E。纬度偏差为 -0.0040°。该地位于印度。
经度偏差最大处在 54.02°N, 131.40°E。经度偏差为 0.0102°。该地位于俄罗斯。

这体现了 GCJ-02 的精妙:误差最大地方的都在国外。也就是说,GCJ-02 把在国内的误差控制在了较小的水平。

借助计算地球表面任意两点距离的算法,我们可以计算误差的具体距离。这个过程中考虑到了地球是个椭球体。不过由于偏移量较小,我们可以部分地忽略两点是球面而不是平面这个事实。

计算结果是:
在矩形中,便宜最大的点偏移量为 718.3 米,坐标为 53.89°N, 131.40°E。该地位于俄罗斯。
在中国实际领土上:
偏移量最大值为 700 米,位于黑龙江省。
偏移量最小值为 18 米,位于青海省德令哈附近。

下面贴上两张图:第一张是进行计算的矩形内偏移量分布情况。颜色越深,偏移量越大,颜色越浅,偏移量越小。


第二张是贴上了政区图后的结果。红色为偏移量大,蓝色为偏移量小。


两张图都使用 Web Mercator 投影法(和墨卡托投影法类似,但忽略了地球的椭率)。第一张图略有拉伸。

结论



  • 流传的转换算法完全正确属实。可以很容易实现自动将 .kml、.gpx 等等坐标和轨迹文件自动转换为另一个坐标系下适用的。

  • 所谓“GCJ-02 偏移有两三公里”的说法是不靠谱的。事实证明 WGS-84 到 GCJ-02 所带来的偏移不超过 700 米。


评论

Devin 2013-12-21 09:41:05
高人啊
fail 2014-01-10 16:56:18
你好,感谢你的文章!依据网上流传的WGS->GCJ算法,我尝试了对覆盖全国的一千多万的POI进行加偏运算,由于该一千多万量的点,本身存在WGS与GCJ两种坐标,所以通过算法运算,再与原始的GCJ进行距离分析,最大误差也只有11m,但这也是比较郁闷的,有1千万点量距离小于1m,但发现云南、西藏、新疆等,距离就比较大了,,,,请问您知道这个是什么回事吗?
草陌博客 2014-01-11 21:51:19
博主很久没更新了
foo 2014-11-03 17:17:11
学习了,顺便吐槽一下那个类名:EvilTransform
。。。
LeadersFirst 2015-09-22 16:51:01
你好请问可否全文转载?
kmxz 2015-10-03 21:25:38
@LeadersFirst

没问题。