转载:关于流传的 WGS-84 至 GCJ-02 转换算法
虎哥2014/07/23软件综合 IP:山东
本文是转载的,写于谷歌退出中国之前。

前言

在夏天去西藏之前,我曾仔细查询了前往西藏的路线、海拔等等资料。由于那些地方的地图稀缺,很多时候我不得不记录下关键路口、山口等地的经纬坐标。从 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,可以通过 XXXXXXXXXXXXXXX 访问
二为 Google Ditu,可以通过 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 访问
显而易见,第二个版本是为了中国准备的。

可以注意到,在 Google Ditu 下,卫星图和公路图是重合的,没有任何问题;然而在 Google Maps 下,看中国内的区域,在同一个地方从公路图切换到卫星图时,可以发现发生了较大偏移。借助 Google 官方的经纬度标注工具,用下面这张图加以说明:
compare.jpg
很容易看出,在 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 的转换算法出自XXXXXXXXXXXXXXXXXXXXXXXXXXX/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 米,位于青海省德令哈附近。
下面贴上两张图:第一张是进行计算的矩形内偏移量分布情况。颜色越深,偏移量越大,颜色越浅,偏移量越小。
origianal_result.jpg
第二张是贴上了政区图后的结果。红色为偏移量大,蓝色为偏移量小。
overall_result.jpg
两张图都使用 Web Mercator 投影法(和墨卡托投影法类似,但忽略了地球的椭率)。第一张图略有拉伸。

[修改于 9年2个月前 - 2015/11/03 03:28:47]

来自:计算机科学 / 软件综合
31
已屏蔽 原因:{{ notice.reason }}已屏蔽
{{notice.noticeContent}}
~~空空如也
虎哥 作者
10年5个月前 IP:未同步
704133
这个所谓的加密,对于老外来说一点保密效果都没有,国内爱好者都可以破解。

但对于国内用户来说,用“罄竹难书”这个成语来形容他的坏处,一点也不为过。

当年汶川地震救灾的时候,我就好几次因为前方测量的经纬度不准,把救援队调去了错误的地方。

当然,也有好处。那就是找直升飞机的时候,有人看到了一点可疑的东西,通报上来一个经纬度。然后,为了抢头功的某军队,竟然监听了救援队的通信,在我们核实现场的可信度之前,就火速下达命令派部队去找——显然,这个经纬度是“加偏”的。
+1
科创币
20!Dopaminor
2014-07-23
都要笑出泪了
引用
评论
2
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
xiaobai152
10年5个月前 IP:未同步
704134
现在的ios7里的和安卓4.2以上版本中各个地图app还是经过这个算法加偏了么?
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
rudolf
10年5个月前 IP:未同步
704135
和菜刀实名制有异曲同工之妙
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
虎哥作者
10年5个月前 IP:未同步
704145
本帖最后由 虎哥 于 2014-7-23 01:34 编辑

按国家政策,现在仍要加偏。向人民群众出售不加偏的地图是违法行为。

国外水货机内置的图,可能有没加偏的。但是涉及中国,原创者为了脱清干系,通常都说是第三方纠偏的。

比如下图,数字高程模型是用厂商发布的数据纠偏的(当然这是爱好者的“个人”行为,不关厂商的事)。等高线是依据NASA公开的数据算出来的,是无偏的。如果不对数字高程模型纠偏,那么等高线的山头和三维地形的山头就不会重合,非常难看。

1.jpg
+1
科创币
wuncle
2014-07-23
好漂亮的图!
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
20!Dopaminor
10年5个月前 IP:未同步
704270
本帖最后由 20!不曾存在 于 2014-7-23 12:25 编辑

有次用百度地图(也可能是谷歌地图 )电脑网页版,歪了几个屏幕,都要疯了,                     最基本的信任都没有了
只防君子不防贼..
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
89759
10年5个月前 IP:未同步
704274
赞一个!
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
89759
10年5个月前 IP:未同步
704275
国内连万分之一地形图都保密。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
89759
10年5个月前 IP:未同步
704276
怪不得用谷歌地图还走错路!
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
freeflying
10年5个月前 IP:未同步
704278
哈哈,2008年从国外带回一个GPS定位器,测得经纬度和地图有几百米误差,先怪GPS定位器不准,随后用谷歌卫星图坐标比较,吻合上。中国国内地图真实坐标害人,而且害的都是自己人。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
改装pcp
10年5个月前 IP:未同步
704290
汽车导航仪有这回事吗?用着好像不错。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
安忆lucky
10年5个月前 IP:未同步
704295
怪不得一次出去玩用谷歌地图还走错了路,万恶的天朝-_-||
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
安忆lucky
10年5个月前 IP:未同步
704297
害得我多走了10多里路~
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
89759
10年5个月前 IP:未同步
704343
安忆lucky 发表于 2014-7-23 12:30
怪不得一次出去玩用谷歌地图还走错了路,万恶的天朝-_-||


我去办正事,让地图误了,干。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
89759
10年5个月前 IP:未同步
704354
真不知天朝这样做有什么意思!
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
freeflying
10年5个月前 IP:未同步
704393
改装pcp 发表于 2014-7-23 12:27
汽车导航仪有这回事吗?用着好像不错。


国内汽车导航仪座标都经过加密,地图也相应修正,行车没问题。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
23323
10年5个月前 IP:未同步
704551
万恶的天朝,我还以为Google maps出问题了
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
xiaobai152
10年5个月前 IP:未同步
704570
rudolf 发表于 2014-7-23 00:54
和菜刀实名制有异曲同工之妙


洋火都不让用了[s:12]
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
rudolf
10年5个月前 IP:未同步
704736
掩耳盗铃 劳民伤财
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
1176764177
9年11个月前 IP:北京
744230
现在淘宝出售的gps模块也是加偏的么?、
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
落叶兄
9年11个月前 IP:辽宁
744234
引用 1176764177:
现在淘宝出售的gps模块也是加偏的么?、
我刚买一个模块,用在火箭巡回上的,感觉平台就是百度地图。家里干的小买卖在地图上有显示,在家里平台上显示距某某12m,难道是我家当时位置上报错了?。。。我买的这个精度不错,实测应该偏差不到8m
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
1176764177
9年11个月前 IP:北京
744257
引用 落叶兄:
我刚买一个模块,用在火箭巡回上的,感觉平台就是百度地图。家里干的小买卖在地图上有显示,在家里平台上显示距某某12m,难道是我家当时位置上报错了?。。。我买的这个精度不错,实测应该偏差不到8m
今凌晨我和群里一哥们实验用谷歌地图对比表示误差咋30m左右
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
落叶兄
9年11个月前 IP:辽宁
744261
引用 1176764177:
今凌晨我和群里一哥们实验用谷歌地图对比表示误差咋30m左右
?
那一个?
虎哥不是说谷歌地图有俩吗
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
落叶兄
9年11个月前 IP:辽宁
744263
引用 1176764177:
今凌晨我和群里一哥们实验用谷歌地图对比表示误差咋30m左右
咱俩说的不是一个东西吧。。。我买的是这个,方便XXXXXXXXXXXXXXXXXXXXXXX/XXXXXXXm?spm=a230r.1.14.8.vYm5VX&id=19674569262&ad_id=&am_id=&cm_id=140105335569ed55e27b&pm_id=&abbucket=18
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
1176764177
9年11个月前 IP:北京
744429
引用 落叶兄:
咱俩说的不是一个东西吧。。。我买的是这个,方便XXXXXXXXXXXXXXXXXXXXXXX/XXXXXXXm?spm=a230r.1.14.8.vYm5VX&id=19674569262&ad_id=&am_id=&cm_id=14010...
谷歌现在除了中国就一个了吧
不太清楚
他那个是原子家的。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
虎哥作者
4年5个月前 IP:四川
882093

好多年没关心这事了,现在上梁不正下梁歪,学坏容易学好难,大恶棍百度也采用了类似办法,在GCJ02的基础上再套了一层壳,搞出了个BD09,其它地图公司也纷纷效仿,变成了商业壁垒。

这个东西糟糕之处在于很难求得反函数,常常只能依靠WGS84到GCJ02方向的转换函数来逼近。尽管可以无限精确,换句话说对保密没有什么卵用,但还是需要开销一些计算量,对于一般使用来说就很不方便。

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
phpskycn
4年5个月前 IP:云南
882189

Google中国的地图服务(XXXXXXXXXXXXXX)去年9月被下线了

一个临时的方法是访问XXXXXXXXXXXXXXXXXXXXX//maps,目前还可以获得修正过的中国卫星地图,但这并不是正式的方案

真的要解决这个问题还是得发展XXXXXXXXXXXXX之类的方案,大家多做标记


引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
warmonkey
4年5个月前 IP:广东
882205
引用落叶兄发表于20楼的内容
引用 1176764177:现在淘宝出售的gps模块也是加偏的么?、 我刚买一个模块,用在火箭巡回上...

查询GPS经纬度(硬件坐标)在地图的位置,可以用XXXXXXXXXX,已经修正。

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
westwarded
4年5个月前 IP:广东
882234

我看见过军队测绘的部队进行测绘。他们回对范围里面的目标地形进行每年的更新。拿着去年的地图,然后再到同样位置,进行对目标从新测量。包括了地形的变化和建筑物的改编。甚至山被挖了个沟都要测绘出来。(反正我不懂:估计这种测绘估计是将地形变化进行修正,将卫星数据与现场测量进行修正。但是真实的地形数据对弹道导弹能否准确命中和巡航导弹地形匹配肯定是至关重要,卫星回来数据肯定是有畸变不能提升导弹精度)。所以从这个角度来说就能明白为何国家禁止未经许可的测量定义得那么高的罪。尤其是中国这种地形非常复杂的国家对巡航导弹来所匹配最好了。这就为啥海湾战争战斧要飞山区不飞沙漠这些地方。

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
虎哥作者
4年5个月前 修改于 4年5个月前 IP:四川
882236
引用westwarded发表于28楼的内容
我看见过军队测绘的部队进行测绘。他们回对范围里面的目标地形进行每年的更新。拿着去年的地图,然后再到同...

为何国家禁止未经许可的测量定义得那么高的罪”,不但语句不通,还一如既往的胡说八道。刑法上并无非法测绘罪,何来“那么高的罪”?

不论如何测绘,只有测绘成果能被归类到“国家秘密”、“军事秘密”范畴,才按照有关国家秘密、军事秘密的罪刑规定定罪处罚。非法测绘谋取利益,也可能按非法经营处理。

违法不等于犯罪,违法也可以行政处罚。

1.png

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
MayDay11neo73
4年5个月前 IP:浙江
882259

管的越来越严了啊。

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
WernerPleischner
6个月8天前 修改于 6个月8天前 IP:广东
932083

以前为了做带有网页地图的UI程序费劲找过某熊掌地图坐标的反变换方法,今日读了某乎上的老文才惊觉这种“加密”算法的脱裤子放屁之处,虽然专门精心设计得特别复杂来避免求解,但它为了不影响日常使用(比如日常导航方向和和距离总不能变得太离谱吧),小尺度下必须保持线性变换,还不能带缩放和旋转,于是这不就只剩下平移了嘛,这下已知变换函数时近似求逆就好办了(即使没有函数解析式用不多的点做个线性内插想必也是足够的)。

整理一下推论。当 x 与 y 临近时,███-██ 有以下性质:

局部线性:f(y)-y≈f(x)-x

相对位置近似:f(y)-f(x)≈y-x

初步反解:x≈f(x)-(f(f(x))-f(x))

精确反解:x≈y-(f(y)-f(x))

因为 ███-██  有局部线性性质,所以有 f(f(x))-f(x)≈f(x)-x,从而 x≈f(x)-(f(f(x))-f(x)),这就是  ███-██ 的反解原理。这样反解已经能达到米级精度。

转载源:XXXXXXXXXXXXXXXXXXXXXXXXXX/p/79027966

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论

想参与大家的讨论?现在就 登录 或者 注册

所属专业
上级专业
同级专业
虎哥
专家 进士 学者 机友 笔友
文章
1551
回复
13396
学术分
39
2005/08/24注册,46分51秒前活动

刘 虎

创新工程局主席

主体类型:个人
所属领域:无
认证方式:身份证号
IP归属地:未同步
文件下载
加载中...
{{errorInfo}}
{{downloadWarning}}
你在 {{downloadTime}} 下载过当前文件。
文件名称:{{resource.defaultFile.name}}
下载次数:{{resource.hits}}
上传用户:{{uploader.username}}
所需积分:{{costScores}},{{holdScores}}下载当前附件免费{{description}}
积分不足,去充值
文件已丢失

当前账号的附件下载数量限制如下:
时段 个数
{{f.startingTime}}点 - {{f.endTime}}点 {{f.fileCount}}
视频暂不能访问,请登录试试
仅供内部学术交流或培训使用,请先保存到本地。本内容不代表科创观点,未经原作者同意,请勿转载。
音频暂不能访问,请登录试试
支持的图片格式:jpg, jpeg, png
插入公式
评论控制
加载中...
文号:{{pid}}
投诉或举报
加载中...
{{tip}}
请选择违规类型:
{{reason.type}}

空空如也

加载中...
详情
详情
推送到专栏从专栏移除
设为匿名取消匿名
查看作者
回复
只看作者
加入收藏取消收藏
收藏
取消收藏
折叠回复
置顶取消置顶
评学术分
鼓励
设为精选取消精选
管理提醒
编辑
通过审核
评论控制
退修或删除
历史版本
违规记录
投诉或举报
加入黑名单移除黑名单
查看IP
{{format('YYYY/MM/DD HH:mm:ss', toc)}}