几何大地测量采用一个同地球外形最为接近的旋转椭球作为参考模型，研究地球椭球数学性质、椭球投影数学变换、椭球面测量计算等。几何大地测量可为航海提供定位、定向、导航信息[1]，也可为目标运动要素解算提供依据。大地主题解算[2]是几何大地测量问题中的一个重要方面，比较经典的方法有直接对大地线微分方程进行积分的方法和利用地图投影理论解算大地问题的方法[3-4]。类似贝塞尔法[5]等一些主题解算方法都是利用地图投影理论，先将椭球面问题对应到球面问题加以解决[6-7]。因此讨论球面距离计算是大地主题解算的一个重要基础。
球面距离是球面上两点之间的最短距离，是经过这两点的大圆在这两点间的一段劣弧的长度。球面距离较常用的计算公式有大圆公式(Great-circle distance formulas)[8-10]、Haversine公式[10-12]等。类似大圆公式的另一种有代表性的计算方法是使用法向量代替经纬度来描述位置，通过采用点积、叉积或两者组合的三维矢量代数进行球面距离求解[13]。另可由弦长得到大圆距离，具体思路为：在地球球体上，通过三维空间连接两点之间的线段即为大圆弦长，两点之间的中心角可以根据弦长确定，大圆距离与弦长成正比。Haversine公式是球面三角学中更通用公式的一个特例，通常其与大圆公式的计算精度相当，但在某些特殊情形，Haversine公式具有更高的计算精度。以上这些大圆距离计算方法本质上都是基于球面三角关系[14-15]得出的，本文从另一个全新的角度首先针对大圆弧固有特性论述了球面大圆弧的A·r(A：Azimuth，r:Radius)约束，然后基于此约束提出球面距离的积分解算模型，并给出模型的解析解，分析奇异点(沿大圆弧的同一走向导致纬度变化量的符号发生改变的点)对模型求解的影响，进一步完善了解算模型。球面距离模型的解析解可应用于大地主题解算，同时可为航海计算及辅助决策等问题提供基础支撑。
1 球面A·r约束
经过球心的平面截得的球面圆为球面大圆。球面大圆的圆弧、劣弧分别称为球面大圆弧、大圆劣弧。通过球面上任意两点(非直径对点)有且仅有一个大圆。任意大圆上的两点A0和B将大圆分为两部分，都叫做以A、B点为端点的大圆弧。
如图1所示，设K为球面大圆弧上一点，其地理坐标为(φ,λ)，球面方位角为A。取KK1为大圆弧弧素dS，则K1点坐标为(φ+dφ,λ+dλ)，球面方位角为(A+dA)，纬圈弧素为KK′；过K、K1分别作子午线切线KT、K1T(可视两切线同交于短轴的延长线于T点)；
微小三角形KK1K′可视为平面直角三角形，于是有：
(1)
式(1)中，R为地球半径。
图1 A·r约束的球面示意图
即:
(2)
由平面直角三角形TKK1可得：A+dA=A+∠K1TK,即dA=∠K1TK。
因为K′十分接近K1，所以可视K′也在切平面K1TK上，于是有：
(3)
在直角三角形OKT中，KT=Rcotφ，代入式(3)得：
dA=sinφdλ
(4)
结合式(2)可得：
dA=tanAtanφdφ
(5)
将平面K1K′O从球体中分离，如图2所示。
图2 纬度变化引起的Rdφ与dr关系示意图
定义纬度增大的方向为“+”，随着纬度的增大，等纬圈半径减小，故dr与dφ异号。
令
为子午圈上一段距离极小的弧素，故可视K′也在K1点的切线K1T上，令K1纬圈半径为r1，得：
sinφ·Rdφ=r1-r=-dr
(6)
故
(7)
即 rcosAdA+sinAdr=0。
积分得A·r约束条件为：
rsinA=C
(8)
式(8)中，C为与大圆弧对应的常量。
因此大圆弧上任意一点均满足该条件，即球面A·r约束为：
r1sinA1=r2sinA2=…=rnsinAn=C
2 基于球面A·r约束的大地距离解算模型及其解析解
如图3所示，pN为北极，p1、p2分别为球面上两点，其球面纬度分别为u1、u2，球面经差为Δω，两点之间的球面距离|p1p2|为δ，α1、α2为对应的球面方位角。
图3 球体大圆
假设在大圆弧p1p2取任意一段弧长为无限小的圆弧
之间的球面纬差为
在过px点的子午线上并与
同纬度，A为px点的球面方位角，三角形
可近似直角三角形，故圆弧
的长度可表示为：
(9)
对式(9)积分，求得大圆弧p1p2的长度为：
(10)
大圆弧上每一点的方位角A不同，但都满足式(8)球面A·r约束条件。
由p1点的球面纬度与球面方位角可求得大圆弧常数C，即：
C=R·cosu1·sinα1
(11)
将式(8)、(11)代入式(10)可得：
(12)
对上式的不定积分进行换元积分可得：
(13)
最终得到大圆弧长计算公式为：
(14)
3 模型奇异点的存在问题
这种积分求取两点间球面距离的方法不需要两点的经度或经度差信息，只需两点的球面纬度及其中一点的大圆弧方位角。但此种方法适用于两点之间大圆弧上不存在大圆顶点的情况，下面对模型存在奇异点的原因进行分析。
如图4所示，pa、pb为球面大圆弧上两点，pv为此大圆弧其中一个顶点，则可得两点间球面距离公式为：
(15)
式(15)中第二部分积分在积分区间[uv,ub]上，du始终为负值，而被积函数式中为正，因此式(15)又可表示为：
(16)
由此得出结论：若大圆弧上的连续点的纬度值非单调，则式(12)中的纬度变化量du时正时负，大圆距离的积分结果会偏离实际值。
图4 奇异点的示意图
因此式(12)需改写为：
(17)
由大圆特性可知，南北半球各有一个大圆顶点，且两者的纬度相同，经度相差180°。假设两点之间大圆弧p1p2上存在一个大圆顶点pv(uv,ωv)，进一步有：
δ=
(18)
在大圆顶点pv处，球面大圆与子午线正交，球面方位角αv=90°，根据式(11)可得pv的纬度有如下关系：
cosuv=cosu1·sinα1
(19)
由球面三角形公式有：
cosα1=cosΔωv1·cos(π-αv)+
(20)
cosα2=cosΔωv2·cos(π-αv)+
(21)
可得大圆顶点pv与点p1、p2的经差分别为：
(22)
(23)
由Δωv1、Δωv2、Δω可知大圆顶点是否在p1p2上，有：
(24)
4 解算模型的有效性检验
综上所述，可最终确定球面大圆距离的求解公式：若大圆顶点pv在两点之间大圆弧上，球面距离求解公式为式(14)；若大圆顶点pv在两点之间大圆弧之外，球面距离求解公式为式(18)。
为验证积分模型的有效性，将本模型的计算结果与球面三角公式解算的大圆距离在Matlab中进行比对。固定其中一点p1(120°E,86°N)，p2经度变化范围为[120°E,122°E]，纬度变化范围为[83.5°N,85.5°N],得大圆距离解算值及解算误差分别如图5、图6所示。此处取地球半径R=6 377 830 m。
图5 经纬度变化区间内的大圆距离
图6 经纬度变化区间内的大圆距离比对差值
通过研究在p2点的经纬度变化范围内积分模型解算的球面距离如图5所示，与球面三角公式的比对差值如图6所示。第二种情况研究了高纬度下球面距离的准确性问题，球面距离最大约为2.7×105 m，随着球面距离的增加，差值略有增大，但差值最大值控制在1×10-7 m左右。可见此模型求解球面距离具有非常高的可靠性。
上述情况中不涉及奇异点问题，为研究奇异点存在情况下模型的解算准确性，选取了如表1所示的几种典型奇异点存在情况，计算球面要素并进行大圆顶点的位置判断，给出模型的解算误差。表1中第一列是在两点纬度相差非常小的情形下得到的距离求解结果，第二列是在高纬度下同纬度两点的解算结果，第三列是两点长距离情形下得到的解算结果，3种情形下距离解算误差都控制在1×10-7 m以内。经验证在两点同经度或同纬度及长距离的情况下积分解算模型依然有效。
表1 几种典型奇异点存在情况下球面元素
p1、p2点坐标(120°E,33°N)、(122°E,32.99°N)(120°E,86°N)、(160°E,86°N)(120°E,33°N)、(160°E,32.99°N)pv点的纬度/(°)33.000 554 170 158 20386.240 515 632 788 30234.642 658 309 940 636pv点与p1点的经差/(°)0.372 859 601 084 75620.000 000 000 000 01119.969 922 780 667 616pv点与p2点的经差/(°)1.627 140 398 915 28920.000 000 000 000 01120.030 077 219 332 377式(18)得到的大圆距离/m1.867 228 601 787 295e+0053.043 549 343 853 731e+0053.711 206 903 348 704e+006应用文献[8]大圆距离公式得到的结果/m1.867 228 601 787 035e+0053.043 549 343 852 938e+0053.711 206 903 348 709e+006大圆距离解算差值/m2.601 882 442 831 993e-0087.933 704 182 505 608e-008-4.190 951 585 769 653e-009
5 结论
提出的基于A·r约束的大地距离解算模型能够有效解算大地距离，且计算精度与球面三角公式相当。球面距离模型应用不受使用条件限制，在高纬度地区或同纬度、同经度等情况下能保证较高的精度。球面距离模型的解析解可应用于大地主题解算，同时可为航海计算及辅助决策等问题提供重要支撑。
[1] 潘瑞鸿,徐胜红.基于几何特性的多无人机协同导航算法[J].兵器装备工程学报,2017,38(10):55-59.
[2] VINCENTY T.Direct and inverse solution of geodesic on the ellipsoid with application of nested equation[J].Survey Review,1975,22:88-93.
[3] 施一民,范业明,朱紫阳.新型大地坐标系中的大地主题解算[J].同济大学学报(自然科学版),2006,34(2):213-216.
[4] 施一民,朱紫阳,方胤祺.测地主题正反解解算[J].测绘工程,2003,12(1):9.
[5] BOERSMA J,GLASSER M L.A differentiation formula for spherical Bessel functions[J].Journal of Physics A:General Physics,2005,38(8):1687-1690.
[6] 史国友,周晓明,贾传荧.贝塞尔大地主题正解的改进算法[J].大连海事大学学报,2008,34(1):16-19.
[7] 周江华，苗育红,成文生,等.贝塞尔大地问题反解的高效率算法[J].测绘学报，2002，31(2)：108-111.
[8] KELLS L M,KERN W F,BLAND J R.Plane and Spherical Trigonometry[M].McGraw Hill Book Company,Inc,pp.323-326.Retrieved July 13,2018.
[9] MONTICONE L C,SNOW R E,BOX F,et al.Minimizing Great-circle Distance Ratios of Undesired and Desired Signal Paths on a Spherical earth[J].IEEE Transactions on Vehicular Technology,2009,58(9):4868-4877.
[10] SHEBA S,RAMADOSS B,BALASUNDARAM S R.Geo distance-based event detection in social media[J].International Journal of Computational Intelligence Studies,2015,4(1):87-101.
[11] ANISYA,SWARA G Y.Implementation of harversine formula and best first search method in searching of tsunami evacuation route[J].IOP Conference Series:Earth and Environmental Science,2017,97(1):1-7.
[12] MARKOU M T,KASSOMENOS P.Cluster analysis of five years of back trajectories arriving in Athens,Greece[J].Atmospheric Research,2010,98(2):438-457.
[13] GADE K.A non-singular horizontal position representation[J].The Journal of Navigation.Cambridge University Press.2010,63(3):395-417.
[14] MAKOW D M .A spherical triangle computer for marine and air navigation[J].IEEE Transactions on Aerospace and Navigational Electronics,1963：324-328.
[15] YOICHI M.Seven types of random spherical triangle in Sn and their probabilities[J].Lecture Notes in Computer Science,2008,4535(1):119-126.