航行中相对位置设计与计算精度 – Aeroinformatics / 航空信息研究

相对方位与距离是测量和领航中极为常用的位置描述格式。航行中对航迹的定义或障碍物的描述，经常涉及根据一个已知点及其相对方位和距离求另一个点坐标的计算。这时，已知点的位置往往是从公布的航行资料或测绘的成果表格中取得的经纬度(地理坐标)形式，而方位和距离则是来自VOR/DME导航设施(公布)或全站仪(实测)定义的局部平面坐标系下获取的。显然不能在经纬度上加上平面距离增量，而飞行程序设计实务中会用到一个手手相传的MFC三无exe工具来实现相关转换，背后原理与可靠性不得而知。

与之相关联的，是设计图件在AutoCAD (.dwg)与Google Earth (.kml)之间的转换。实务中，见到过一个名为acad2kml的AutoLISP插件，但转换往往只能是单向的，难以将kml转换和导入进AutoCAD。后者自带的地理位置功能在常用的版本中看起来只能以一个点为参照，将图纸映射到某个自定义的投影坐标系中的对应位置。

对于相对经纬度一定方位和距离的点的计算，考虑三种朴素的思路：

先把已知经纬度投影到合适的UTM坐标系中，再在UTM投影坐标系平面内计算三角函数关系，后再逆投影回经纬度地理坐标系。这需要用到ArcGIS/QGIS等大型软件，或ogr2ogr命令行工具，或proj4等编程函数库，有一定门槛且步骤繁琐。
先将距离按平均半径6371km转化为球面弧度，再运用球面三角学知识(如球面正弦定理)，在球面上根据相对方位，计算出距离弧度对应的经纬度增量。这需要对球面三角学有个基础认识，找到合适材料¹认真学一小时应该搞得清楚。
如果要考虑参考椭球模型，那就只能请出Vincenty (1975)的正向测地线公式²了。精度够高，只有接近对跖点的情形不太适用，当然也还有相关的精确方法。因为有各种函数库实现³，就暂且忽略其中的数学细节。

考虑AIP China(题外话，最近情报中心开放了AIP访问，可以成为一个公开的权威数据源)中的一个算例：南翔NDB(PK)位于虹桥ARP磁方位3度9630米处，同时也公布了精确到0.1分的经纬度，可以用来验证。尽管数据是公开的且加偏的，此处还是模糊处理不直接出现经纬度坐标点了。

解法一我的实际操作是用了epsg.io的坐标转换，外加自己手按计算器。

解法一变量	值(东)	值(北)
参考点投影坐标	341365.7316m	3452591.9855m
相对位置增量	-464.8273m	9618.7751m
相对位置投影坐标	340900.9042m	3462210.7606m
相对位置经纬增量	-0.006404°	0.086689°
相对位置误差	-0.08426′	0.00133′

解法二我的实际操作是用自己粗浅的球面三角现学现卖，推导了一下公式（三角函数和反函数都是度，非弧度）：

\( Δ\mathrm{lng}=\arcsin\cosθ\sin({πρ}/{180R}) \)

\( Δ\mathrm{lat}=\arcsin\sinθ\sin({πρ}/{180R}) \)

解法二变量	值(东)	值(北)
相对位置经纬度增量	-0.004180°	0.086504°
相对位置误差	0.04918′	-0.00978′

解法三采用pyproj.Geod.fwd实现，可以看到在东西方向上精度直接高了一个数量级，南北方向上也比球面优一些。

解法三变量	值(东)	值(北)
相对位置经纬度增量	-0.004882°	0.086754°
相对位置误差	0.00710′	0.00526′

需要注意的是，以上个例不见得具有代表性，比如纬度低、远离投影中央经线。同事，所谓误差是与AIP公布值的对比，而公布值本来就是有误差(舍入和人为)的，南翔这个误差(公布值与实务工具软件)在0.01′左右，是凑巧南翔台确实在接近0.1分的位置，让这个算例更能体现出各个思路的计算精度。也就是说，在AIP公布的0.1分的精度下，这三个思路基本都可以满足需要。从实际操作的方便程度而言，如果手工计算可以选球面正弦定理的方式，计算简明、精度适中；如果编程开发可以选Vincenty公式的有关实现，开箱即用、精度较高。