AIS和雷达信息融合(1)坐标转换

AIS和雷达信息融合（1）坐标转换

虽然很多⽂章资料的坐标转换是说的⾼斯-克吕格投影等，但最后的实现仿真步骤很多是把AIS的经纬度转换为和雷达⼀样的极坐标

⾼斯-克吕格投影我⼜按照公式写个程序，但是貌似错了，后来查找资料，分带和中央经线把我搞得越来越头晕。

现在做的就是将AIS的WGS-84⼤地坐标系中的经纬度转换为⽬标船相对本船的距离和⽅位，与雷达数据格式统⼀。

防⽌被删除，我复制过来。

最近⾃⼰做的⼀个⼩东西要⽤到经纬度⽅⾯的计算，查遍中⽂⽹页见到的要么基本上是⼀帮惜字如⾦装⼤爷的“砖家”，要么就是像贴膏药⼀样，啪，⼀⼤堆代码

往上⼀贴，⼀点说明都没有，让⼈看不懂，有的看了半天看懂了，结果他⽤的公式要么有使⽤局限（但没有半点声明）要么根本就是个错的。所以现在将⾃⼰⼏

天学习来的在这⾥总结⼀下，⽅便后来⼈少⾛弯路。

这⾥主要解决四个问题：

1、已知两点经纬度，求⼀点相对于另⼀点⽅位⾓；

2、已知两点经纬度，求两点间距离；

3、已知⼀点经纬度及与另⼀点距离和⽅位⾓，求另⼀点经纬度；

4、问题1与问题2的简化算法。

注：简化算法的运算量和对系统的运算精度要求都⼤⼤降低，但只在短距离内（⾼纬地区建议10km以下）可以保证精度，除简化算法之外

的算法可适⽤于地球上任意两点。这⾥只是出于便于理解的⽬的来解释“原理”，具体到不同的编程环境还要⾃⼰做化简和注意单位。

在求算前我们先对符号及单位进⾏约定：

此处设定求B相对于A的⽅位⾓，即A为当前位置，B为⽬标位置

Aj：A点经度

Aw：A点纬度

Bj：B点经度

Bw：B点纬度

北纬为正，南纬为负；东经为正，西经为负

经纬度使⽤度，°，⾮度分或度分秒。

度数未加说明均采⽤⾓度制

R：地球平均半径

说到这⾥，⼈们或许会⾸先想因为地球是个球体，如果AB两点⾜够近（如相距1Km）可以当做平⾯三⾓形已知两临边求夹⾓，把两点的经

度、纬度各⾃做差，差值作为两边的长度，再⽤反正切函数⼀算就得到了⾓度，似乎很有道理，但是如果将计算结果与实际测量值做⽐较，

就会发现⽐较⼤的误差，⽽这种误差在某些地区甚⾄能⽤普通的量⾓器测量出来，其实这已经不叫误差，⽽叫错误了。这种近似利⽤平⾯⼏

何知识解决问题的算法只适合于低纬度地区（例如南北纬10度），如果在⾼纬度使⽤这种⽅法，偏差会⽐较⼤，并且这种偏差会随着纬度

的升⾼⽽⼤幅变⼤。例如，在北纬10度，AB两点经纬差0.001时，近似算法与真实值之间的偏差为0.44度。当纬度改为北纬40度其他不

变，这时的偏差就达到了7.55度。

为什么会这样呢？其实原因就在于经线、纬线划定不同。如果把地球简化成⼀个球体，每条经线的长度都等于球体周长的⼀半，每条经线均

在两极相交。但是，纬线之间互不相交，纬线所围成的圆均为“同轴”圆，所以每条纬线的长度会因纬度的不同⽽不同，也就是，纬度不

同，1经度差所对应的球⾯距离是不⼀样的，例如，在⾚道处，1纬度的跨度约为111.3千⽶，1经度的跨度也是约为111.3千⽶；在北

纬40附近，1纬度的跨度没有变，但1经度的跨度却变为85.3千⽶。（转不过弯来的可以去看地球仪，看看就了然了~）

——————————————————卖关⼦环节结束，正⽂开启———————————————————

那么，有没有⼀种对地球任何区域任意两点都普适的求⽅位⾓的⽅法？

答案是肯定的。

现在我们开始计算

第⼀步：在知道AB点经纬度后，要⽤到第⼀个公式，三⾯⾓余弦公式，

A~OC~B是⾯AOC与⾯BOC的⼆⾯⾓，为了⽅便，写成这个样⼦

这⾥我们将已知数据代⼊，公式便写成：

没错，⼆⾯⾓A~OC~B的度数就是两点经度之差

第⼆步：知道了⾓c的余弦值后我们要求得它的正弦值，所⽤的公式就是三⾓函数公式⾥最基本的“扣⽅加赛⽅等于1”的⼀个变形

第三步：求得正弦后，接下来我们要⽤⼀个不太常⽤的公式，球⾯正弦公式

将已知数据代⼊并稍微变形⼀下，公式写为：

这⾥需要注意⼀点，我们⼀开始的假设便是求B点相对于A点的⽅位⾓，因此这⾥是Bj-Aj，不要写反，否则得不到正确结果。

算到这⾥，还没有完，得到的结果并不总符合我们对⽅位⾓的定义，因此要根据B相对于A的位置在四个象限两个轴上进⾏讨论，依据不同情况对计算结果进⾏

不同处理。假设A点固定于原点，则:

B点在第⼀象限，Azimuth=A；

B在第⼆象限，Azimuth=360+A;

B在第三四象限，Azimuth=180-A。

这⾥是我⾃⼰结合实际和程序的计算，得出的坐标系：X轴为经度，Y轴为纬度。

这⾥只说了象限的讨论结果，因为轴上的讨论更复杂些，要结合程序运⾏环境⼀起考虑，考虑的主要因素是系统的计算精度。譬如，在三⾯⾓余弦公式中，

当AB点纬度值相同时，对公式的值起决定作⽤的就是cos(Bj-Aj)这⼀项，当Bj-Aj的值⽐较⼩时，例如0.0001（这在⾚道地区对应的长度为11⽶左右），⽤⼀

般的计算器计算时值为1，这样，后⾯的计算便不可能完成。但是，如果⽤计算机计算则为0.999999999998476913…………。所以，基于以上原因，需要

对轴的“范围进项扩充”，要⽤单⽚机、⼿机运算的尤其要注意。

经过⼀系列计算，最后，就得到了最终结果。

似乎有⼈注意到了，以上的计算都是把地球看成标准的球体，⽽事实是地球是个椭圆，其实，地球的偏⼼率极低，各位可以将此法得到的计

三、第⼆点经纬度的求算

最近在⽹上看到不少⼈在问第⼆点经纬度的求算，所以，这⾥也附加说⼀下求算⽅法。

都应该能想到⼀个最最最笨的⽅法，就是将前⾯两部分⽤的公式联⽴解⽅程。我想只有那些度娘知道⾥的专家会采⽤这种⽅法(因为这种⽅

法费的唾沫最少)。

⾔归正传，解⽅程的⽅法可以，但是运算量极⼤，费时，对于⼀些系统不现实。

另⼀种⽅法，其实就是对⽅位⾓求算的再次运⽤。

已知Aj，Aw，L，R,Azimuth（这⾥的Azimuth依然定为B相对于A的⽅位⾓）。

⾸先求算c，

（注意此处L、R的单位要统⼀）

之后求解a，将已知量代⼊，公式为：

Bw=90-a

Bj=Aj+C

对于某些系统，再单独设定B位于X正负半轴上的值就可以了，有些系统可以返回arctan（X/0）=90。

这种求⽅位⾓的算法代⼊了⼏个坐标值与⾕歌地球⽐对，短距离内偏差在0.1度以下。

最后，关于百度上很多⼈答的需要将WGS84坐标转换成什么北京54、西安80再计算的⾼深⾔论，我认为只要不是搞⼤地测量、⼟⽊⼯程、

导弹发射根本就没有必要，理论分析与计算结果都说明WGS84坐标系统完全可直接⽤于民⽤领域。

按照上⾯的公式我在QT中做了个⼩程序实现，经证明可⽤。