[转载]根据两点的经纬度求方位角和距离，等

最近自己做的一个小东西要用到经纬度方面的计算，查遍中文网页见到的要么基本上是一帮惜字如金装大爷的“砖家”，要么就是像贴膏药一样，啪，一大堆代码往上一贴，一点说明都没有，让人看不懂，有的看了半天看懂了，结果他用的公式要么有使用局限（但没有半点声明）要么根本就是个错的。所以现在将自己几天学习来的在这里总结一下，方便后来人少走弯路。

这里主要解决四个问题：

1、已知两点经纬度，求一点相对于另一点方位角；

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

3、已知一点经纬度及与另一点距离和方位角，求另一点经纬度；

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

注：简化算法的运算量和对系统的运算精度要求都大大降低，但只在短距离内（高纬地区建议10km以下）可以保证精度，除简化算法之外的算法可适用于地球上任意两点。这里只是出于便于理解的目的来解释“原理”，具体到不同的编程环境还要自己做化简和注意单位。

 

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

此处设定求B相对于A的方位角，即A为当前位置，B为目标位置

Aj：A点经度

Aw：A点纬度

Bj：B点经度

Bw：B点纬度

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

经纬度使用度，DDD.DDDDDD°，非度分或度分秒。

度数未加说明均采用角度制

R：地球平均半径

Azimuth：方位角，以真北为0度起点，由东向南向西顺时针旋转360度

A,B,C表示球面上的三个点及球面上“弧线”在该点处所夹的角

a,b,c表示A,B,C三点的对“弧”两端点与地心连线所夹的角（其实这里解释成ABC三点对弧的弧度更方便）

O为球心

L为AB两点间球面距离

（注：因我考虑欠缺，没有注意字母C大小写较难分辨，所以此处提醒读者在后面的公式中注意C的大小写。）

 

一、方位角的求算

已知A、B两点经纬度，如何求出B相对于A的方位角？

——————————————————卖关子环节，可选择性跳过———————————————————

说到这里，人们或许会首先想因为地球是个球体，如果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。

这里只说了象限的讨论结果，因为轴上的讨论更复杂些，要结合程序运行环境一起考虑，考虑的主要因素是系统的计算精度。譬如，在三面角余弦公式中，当AB点纬度值相同时，对公式的值起决定作用的就是cosBj-Aj)这一项，当Bj-Aj的值比较小时，例如0.0001（这在赤道地区对应的长度为11米左右），用一般的计算器计算时值为1，这样，后面的计算便不可能完成。但是，如果用计算机计算则为0.999999999998476913…………。所以，基于以上原因，需要对轴的“范围进项扩充”，要用单片机、手机运算的尤其要注意。

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

似乎有人注意到了，以上的计算都是把地球看成标准的球体，而事实是地球是个椭圆，其实，地球的偏心率极低，各位可以将此法得到的计算结果与谷歌地球（WGS84坐标系统，我说的不是谷歌地图）上的结果进行对比，偏差是非常小的（我测的几个值，最大偏差0.5度）。

二、距离的求算

其实，“眼尖”的或许已经注意到了，第一步的余弦值结果就可以直接用来求算AB两点间的球面距离，用反余弦函数求得c的度数，再将度数转换为弧度，乘以地球半径就得到了两点间的球面距离。

公式为

这里要注意，L的单位与R的单位一致，单位不同的不要忘记换算

短距离（例如100米，30米）使用这个公式，计算出的结果与谷歌地球给出的距离偏差在0.5%以下，长距离计算时，偏差则可以降至0.01%以下。求算的距离越大，偏差越小，就是这个公式的特点，原因不说自明。

PS:对于一些GPS接收机，其数据格式为NMEA-0183，经纬度数据为DDDMM.MMMM，需要将它转换为度，公式为：

经纬度（度）=DDD+MM.MMMM/60

三、第二点经纬度的求算

最近在网上看到不少人在问第二点经纬度的求算，所以，这里也附加说一下求算方法。

都应该能想到一个最最最笨的方法，就是将前面两部分用的公式联立解方程。我想只有那些度娘知道里的专家会采用这种方法因为这种方法费的唾沫最少)。

言归正传，解方程的方法可以，但是运算量极大，费时，对于一些系统不现实。

另一种方法，其实就是对方位角求算的再次运用。

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

首先求算c，

（注意此处L、R的单位要统一）
之后求解a，将已知量代入，公式为：

求得a之后我们求解C，

写到这里应该都恍然大悟了吧，

Bw=90-a

Bj=Aj+C

PS：对于上面两个三角公式的理解，可以想成把A点作为北极点，相当于把方位角求算中的公式原封不动的再用一遍（其实就是再用了一遍）。

四、简化算法

上面讲的算法对于运算精度低的系统简直就像是噩梦。所以这里不得不讲讲简化算法。简化算法的结果在短距离内可以保证精确度，但是在长距离时，因为地球球面的特性，是万万不能用的。

简化算法的基本思路就是将以经纬度表示的球坐标转换成三维直角坐标，再利用平面几何知识去解决。

求算距离：

设：Xa、Ya、Za为三维直角坐标下A点的坐标，B点坐标同样式，

    Ha为A点海拔高度，Hb为b点海拔高度

则

Xa=R+Ha)×cosAw)×cosAj)

Ya=R+Ha)×cosAW)×sinAj)

Za=R+Ha)×sinAw)

Xb=R+Hb)×cosBw)×cosBj)

Yb=R+Hb)×cosBw)×sinBj)

Zb=R+Hb)×sinBw)

（注：此处坐标转换为诱导公式化简后的形式，关于球坐标转直角坐标的原公式可点此查看：球坐标转直角坐标）

 

ΔX=Xa-Xb    
  ΔY=Ya-Yb  
 
  ΔZ=Za-Zb
知道三个坐标轴方向上的差值后再用勾股定理就可以求出两点间距离了，即：

这里说明一下，海拔高度H可有可无，如果有的话，注意H与R单位要统一。普通GPS接收机给出的海拔一般很不准确，所以不推荐使用。另外NMEA规范报文中有两个量涉及到海拔，注意计算。
 
 
求算方位角
 
如果想直接由经纬度求算方位角，则可以避开上面的坐标转换，直接这样求：

上面这个式子的基本思路就是将经度和纬度差转化成地面距离再运用平面几何知识求解，化简之后即为上式（过程从略，可自行推导）
这样，当

B点在第一象限及Y轴正半轴，Azimuth=A；

B在第二象限，Azimuth=360+A；

B在第三四象限及Y轴负半轴，Azimuth=180+A。

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

这种求方位角的算法代入了几个坐标值与谷歌地球比对，短距离内偏差在0.1度以下。

 

最后，关于百度上很多人答的需要将WGS84坐标转换成什么北京54、西安80再计算的高深言论，我认为只要不是搞大地测量、土木工程、导弹发射根本就没有必要，理论分析与计算结果都说明WGS84坐标系统完全可直接用于民用领域。