平面距离只需要计算连接两点线段的长度，通过二维直角坐标使用勾股定理可以很容易地计算出来。由于只涉及乘法、加法与开平方运算，在计算机中也可以保留很高的精度。而球面距离计算则要复杂得多。球面在直角坐标系中实际上是三维结构，除非距离相对于半径非常小，否则不能当作平面处理。

两点球面距离计算是天文或地理学中最常用的计算之一，是目标查找、锥形检索、交叉证认等计算中最基础的运算。比如查找距离某个位置最近的山峰，需要先对附近的山峰都计算一遍距离，取出距离最小的一个。而锥形检索是虚拟天文台的数据检索协议之一，用于在球面上查找与某点的距离小于指定半径的目标列表，几乎所有提供了星表查询功能的系统都有类似功能^[1]。基于位置的星表交叉证认则可视为对一个星表A中所有的源做一遍锥形检索取得与其最近的另一星表B的源，这样就完成了A，B两个星表的交叉^[2]。这些应用的基础都是球面距离计算。

两点球面距离计算，计算的是两点在球面上的最短距离，即球心与两点组成的大圆上两点之间较小的那段圆弧的长度，天文上习惯使用该短圆弧对应的大圆圆心角，称为角间距(Angular Separation)或角距离。球面距离需要使用球面几何进行求解，要大量使用三角函数。而计算机的数值计算精度有限，在对三角函数、反三角函数进行计算时容易出现舍入误差。多个函数的误差积累之后，将导致严重的结果偏离，甚至无法得到正确的结果。前人已经推导了很多数学公式求解球面距离，但在何种情况下使用哪个公式，公式的精度如何，却鲜有人讨论。因此，本文考察天文技术中常用的球面距离计算公式，给出它们的算法，并对比它们的精度, 为球面距离算法的选择提供参考。

1 球面几何方法

球面距离计算可直接通过球面几何推导。较常用的计算公式有大圆(Great-circle)公式^①、Haversine公式^②等。另外还有Vincenty公式^③被应用到Astropy等代码库中，但其复杂度过高，本文在第3.1节中简要提及。

① https://en.wikipedia.org/wiki/Great-circle_distance

② https://en.wikipedia.org/wiki/Haversine_formula

③ https://en.wikipedia.org/wiki/Vincenty%27s_formulae

设有两点赤道坐标为${{p}_{1}}({{\alpha }_{1}}, \text{ }{{\delta }_{1}})~, {{p}_{2}}({{\alpha }_{2}}, \text{ }{{\delta }_{2}})$，求它们的球面角距离d。其中大圆公式为(1)式；当两点距离很小时，也有人使用简化的(2)式^[2-6]；Haversine公式如(3)式。

$ d=\text{arccos}[\text{sin}{{\delta }_{1}}\text{sin}{{\delta }_{2}}+\text{cos}{{\delta }_{1}}\text{cos}{{\delta }_{2}}\text{cos}({{\alpha }_{1}}-{{\alpha }_{2}})]\text{ }, \text{ } $

(1)

$ d=\sqrt{\left[ ({{\alpha }_{1}}-{{\alpha }_{2}})\text{cos}\left( \frac{{{\delta }_{1}}+{{\delta }_{2}}}{2} \right){{~}} \right]^{2}+{{({{\delta }_{1}}-{{\delta }_{2}})}^{2}}}\text{ }, $

(2)

$ d=2\text{arcsin }\!\!~\!\!\text{ }\sqrt{\text{si}{{\text{n}}^{2}}~\left( \frac{{{\delta }_{1}}-{{\delta }_{2}}}{2} \right)~\text{ }+\text{cos}{{\delta }_{1}}\text{cos}{{\delta }_{2}}\text{sin}{{\text{ }\!\!~\!\!\text{ }}^{2}}~\left( \frac{{{\alpha }_{1}}-{{\alpha }_{2}}}{~2} \right)}\text{ }. $

(3)

计算机中通常使用IEEE754二进制浮点数算术标准处理浮点数，其中单精度32位(float或single)真值可精确到小数点后6~9位，双精度64位(double)真值可精确到小数点后15~17位^④。本文使用C++语言编写程序对3个公式的精度进行了对比，结果如表 1。本文使用的单位均为度(^°)，表中有两行数据时，上方的值为单精度计算的结果，下方的值为双精度计算的结果。

④ https://en.wikipedia.org/wiki/Single-precision_floating-point_format

由表 1可见，单精度下，两点之间距离较近时，大圆公式返回的是全0，出现了严重的舍入误差, 而距离较大时误差也较大。双精度下，大圆公式能够返回结果，但比Haversine公式的精度要差。大圆公式的简化形式多数时候与Haversine公式的精度相当。这3个公式共有的一个问题，它们均有可能略微超过准确值，当严格限定一个边界值的时候，可能导致漏源。因而在实际应用时，需要考虑将边界值略微外扩，以将因为计算精度而漏掉的源包含在内。

值得注意的是，大圆公式的简化形式在极点附近时误差非常大。虽然它在上述公式中是计算复杂度最低的，并且在两极之外、距离较小时精度与Haversine公式相近，但是在实际使用该公式时仍应当非常小心。

2 直角坐标系方法

在基于赤道坐标系进行求解之外，也可以使用直角坐标系，通过三维坐标向量对两点球面距离进行计算。可视天球半径为单位1，以天球球心为三维直角坐标系原点，赤道坐标系北天极点方向为z轴方向，天赤道面上赤经为0点的方向为x正向，遵循右手坐标系与x-z相交的轴为y轴方向。从而可将赤经、赤纬转换为三维直角坐标(x, y, z)，其转换关系如(4~6)式⑤。

⑤ https://arxiv.org/abs/cs/0701171

$ x=\text{cos}\delta \text{cos}\alpha , $

(4)

$ y=\text{cos}\delta \text{sin}\alpha , $

(5)

$ z=\text{sin}\delta \text{ }. $

(6)

这样两点距离计算可以直接使用(8)式进行计算。该公式可简单地从Haversine公式中推得，其中Haversine公式根号中的部分(7式)实际上是两点的三维空间直线距离一半的平方。

$ \begin{align} &A=\text{si}{{\text{n}}^{2}}~\left( \frac{{{\delta }_{1}}-{{\delta }_{2}}}{~2} \right)\text{ }+\text{cos}{{\delta }_{1}}\text{cos}{{\delta }_{2}}\text{si}{{\text{n}}^{2}}\left( \frac{{{\alpha }_{1}}-{{\alpha }_{2}}}{~2} \right) \\ &=\frac{1}{{{2}^{2}}}~[{{({{x}_{2}}-{{x}_{1}})}^{2}}+{{({{y}_{2}}-{{y}_{1}})}^{2}}+{{({{z}_{2}}-{{z}_{1}})}^{2}}] \\ &=\text{si}{{\text{n}}^{2}}~\left( \frac{d}{2} \right), \\ \end{align} $

(7)

$ d=2\text{arcsin}~\sqrt{A}\text{ }=2\text{arcsin }\!\!~\!\!\text{ }\left[ \frac{1}{2}\sqrt{{{\left( {{x}_{2}}-{{x}_{1}} \right)}^{2}}+{{\left( {{y}_{2}}-{{y}_{1}} \right)}^{2}}+{{\left( {{z}_{2}}-{{z}_{1}} \right)}^{2}}} \right]. $

(8)

直角坐标系的计算精度方面，本文依旧使用C++语言在单精度及双精度两种环境下对公式进行计算，并与Haversine公式的结果对比，结果如表 2。从表 2可以看出，直角坐标系方法的精度与Haversine公式几乎完全一致，而直角坐标系方法有时甚至优于Haversine公式(见42, 43两行数据)。

直角坐标系方法在进行距离比较时还可以进一步优化，如(9)式：设两点直角坐标为${{p}_{1}}\left( {{x}_{1}}, \text{ }{{y}_{1}}, \text{ }{{z}_{1}} \right)~, ~{{p}_{2}}\left( {{x}_{2}}, \text{ }{{y}_{2}}, \text{ }{{z}_{2}} \right) $，其优点是可以预先计算x, y, z与threshold (设最大角距离为θ)，从而避免在实际进行距离计算时使用非常耗时的三角函数，而只需要三次减法、三次乘法、两次加法、一次比较，大大减少了计算量，非常适用于大规模数据计算。相应地，它的缺点是需要额外占用存储空间来保存x, y, z三个值，在对内存占用要求较苛刻的环境中可能会成为问题。因而，实际的公式选择，仍应视数据规模、计算环境而定。另外，基于三维向量的计算方法，还有法线向量的形式，但其计算过程比本节要复杂，本文将在第3.3节中进行简要描述。

$ {{2}^{2}}A={{({{x}_{2}}-{{x}_{1}})}^{2}}+{{({{y}_{2}}-{{y}_{1}})}^{2}}+{{({{z}_{2}}-{{z}_{1}})}^{2}}\le {{2}^{2}}\text{si}{{\text{n}}^{2}}~\left( \frac{\theta }{2} \right)=\text{threshold }\!\!~\!\!\text{ }. $

(9)

3 计算库

许多天文软件包也包含了球面距离计算功能。本部分选择了几个安装、使用较便捷，应用也比较广泛的库，演示了其使用方法，并进行测试。由于难以分辨软件包中实际使用的是单精度还是双精度，此部分不再对此进行区分，而直接与Haversine双精度结果进行比较。

3.1 Astropy

当前天文界已经非常流行使用Python进行数据处理，其中Astropy^[7]对Python的推广传播功不可没。Astropy进行球面距离计算需要使用astropy.units及astropy.coordinates.SkyCoord计算。调用方法如下，其中d为两点之间的角距离(单位是度)：

c₁ = SkyCoord(ra₁*units.rad, dec₁*units.rad)

c₂ = SkyCoord(ra₂*units.rad, dec₂*units.rad)

d = c₁.separation(c₂).deg

为公平起见，在64位Windows10中的64位python3.6.3重新实现了Haversine算法，并与Astropy的距离计算函数进行了对比，对比结果见表 3，可以看到两种方法基本一致，不分伯仲。

阅读Astropy的源代码可发现其距离计算函数separation实际调用的是astropy/coordinates/angle_utilities.py^⑥中的angular_separation(lon1, lat1, lon2, lat2)函数，其代码实现使用Vincenty公式，即(10)式。该式在所有位置的距离计算都更稳定，但计算复杂度也更高^[8]。从上述计算对比结果来看，Astropy并没有处处比Haversine的精度更高，如(180, 0)一行；当然，Astropy在(42, 43)上的精度也比Haversine更好。但是两者的差距都非常小，几乎可以忽略不计。从比较结果上看，Haversine已经相当精确。当然，若是对计算稳定性有更高要求，可以考虑Vincenty公式。

⑥ https://github.com/astropy/astropy/blob/master/astropy/coordinates/angle_utilities.py

$ d=\text{arctan}~\frac{\sqrt{{{[\text{cos}{{\delta }_{2}}\text{sin}({{\alpha }_{1}}-{{\alpha }_{2}})]}^{2}}+{{[\text{cos}{{\delta }_{1}}\text{sin}{{\delta }_{2}}-\text{sin}{{\delta }_{1}}\text{cos}{{\delta }_{2}}\text{cos}({{\alpha }_{1}}-{{\alpha }_{2}})]}^{2}}~}}{\text{sin}{{\delta }_{1}}\text{sin}{{\delta }_{2}}+\text{cos}{{\delta }_{1}}\text{cos}{{\delta }_{2}}\text{cos}({{\alpha }_{1}}-{{\alpha }_{2}})\text{ }}. $

(10)

3.2 分层三角网格

分层三角网格(Hierarchical Triangular Mesh, HTM) ^[9]是一种对球面的多层次、递归三角分割方法。分层三角网格目前公开提供下载的软件包^⑦中包含了C#代码库及一个SQL Server扩展包。其中Spherical.Htm.Sql.cs文件中包含了函数fDistanceEq (ra1, dec1, ra2, dec2)，可用于计算距离。需要注意的是，fDistanceEq的返回值单位是arcmin，不是degree。针对分层三角网格，本文使用C#实现了Haversine函数并与SphericalHTM v3.1.2的fDistanceEq函数进行对比。如表 4，其结果基本相同。查看源代码发现fDistanceEq实际上也使用了Haversine公式进行计算，只是多了degree到arcmin的转换。

⑦ http://www.skyserver.org/htm/index.html

除了C#程序中可以直接调用fDistanceEq外，HTM软件包提供的SQL Server数据库扩展包中也带有此函数，在数据库中使用非常方便。

3.3 SLALIB与SOFA

SLALIB即Starlink library of positional astronomy routines^⑧，使用Fortran 77实现，其目标是让天文工作者更简便快捷地编写精确可靠的天体测量程序。虽然它是由Fortran 77实现的，但是已经有人为其编写了Python接口pyslalib^⑨，因而可以在已有的Python环境中获得Fortran 77的计算精度。

⑧ http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?slalib.sys

⑨ https://github.com/scottransom/pyslalib

Pyslalib中的slalib.sla_dsep直接对应了SLALIB中的sla_DSEP(A1, B1, A2, B2)球面距离计算函数，在64位的Ubuntu 16.04.4操作系统中使用64位Python3.5.2环境将其与Haversine函数的计算结果进行了对比，如表 5。首先可以看到Linux版的Python Harversine程序计算结果与Windows一致，表明Python3在不同的操作系统中的表现一致。其次pyslalib的计算结果，除了(42, 43)两行有轻微差别，与Haversine完全一致，精度也非常高。查看SLALIB的源代码，可以看到其使用了第2节中的直角坐标系三维向量作为计算单元。但与第2节不同，SLALIB进一步将其转换为法线向量^⑩，即${{\overrightarrow{n~}}_{1}} $与$ {{\overrightarrow{n~}}_{2}}$，再使用(11)式^[10]进行计算。使用向量的主要优势是向量代数可取代部分三角函数的计算，其数值计算稳定性更好，能保持较好的精度，并且在边界点(如南北天极、0~360度边界)也无异常。

⑩ https://en.wikipedia.org/wiki/N-vector

$ d=\text{arctan}\left( \frac{{{\overrightarrow{n~}}_{1}}\times {{\overrightarrow{n~}}_{2}}}{{{\overrightarrow{n~}}_{1}}\cdot {{\overrightarrow{n~}}_{2}}} \right). $

(11)

天体测量中常用的SOFA库^⑪的球面计算函数iauSeps ^⑫也同样使用了(11)式。SOFA即Standards of Fundamental Astronomy，是国际天文学联合会建立并维护的一个权威的基本天文学标准库，包含有ANSI C及Fotran 77两个版本。表 5的最右侧列出了双精度位C++程序调用SOFA计算的结果。

⑪ http://www.iausofa.org

⑫ http://www.iausofa.org/2018_0130_C/SeparationPA.html

4 数据库

关系型数据库均支持SQL语句，因而可以在不同的系统中直接使用SQL语句实现Haversine公式进行距离计算。由于每种数据库使用的数学函数、计算精度略有差别，本文在Microsoft SQL Server 2008，PostgreSQL 9.4.5，MySQL 5.7.18上进行了测试，测试结果如表 6，从表 6可以看到，3个数据库的差距不大，MySQL默认保留的精度更多一些。

除了Haversine公式，数据库中使用直角坐标系方法实际上更方便。数据库对I/O进行了大量优化，可以更好地对数据进行调度。因而，直角坐标系方法虽然需要增加x, y, z 3个字段，但并不会给数据库增加太大的负担，还可以大大降低计算的复杂度，更快地取得海量数据的计算结果。

4.1 PostgreSQL数据库插件

数据库中除了可以直接使用公式之外，也可以使用各种数据库扩展插件进行相关计算，减少使用者编程的麻烦。如PostgreSQL数据库有PostGIS ^⑬, pgsphere ^⑭, Q3C^⑮, H3C ^⑯等多个插件支持球面距离计算。其中PostGIS是专为地理信息系统设计的，虽然也可以在天文上使用，但是需要将米、千米等单位转换回弧度，太过繁琐低效。而pgphere, Q3C^[11], H3C则是专门为天文检索设计的，更为简洁。下文演示了3种插件的使用方法，表 7对这3个插件的结果进行了对比。需要注意的是Q3C和H3C直接使用度进行计算，而pgsphere需要先转成弧度进行计算，再将结果转回角度。

⑬ https://postgis.net/

⑭ https://github.com/akorotkov/pgsphere

⑮ https://github.com/segasai/q3c

⑯ http://cds.u-strasbg.fr/resources/doku.php?id=h3c

SELECT a ra1, b dec1, c ra2, d dec2,

q3c_dist(a, b, c, d) q3c,

h3c_dist(a, b, c, d) h3c,

DEGREES(spoint(RADIANS(a), RADIANS(b))<->spoint(RADIANS(c), RADIANS(d))) pg

FROM (VALUES(0.0, 0.0, 0.01, 0.0), (0.0, 0.0, 0.0, 0.01), (180.0, 0.0, 180.01, 0.0), (42.0, 43.0, 42.01, 43.0), (42.0, 43.0, 42.0, 43.01), (0.0, 90.0, 180.0, 89.99), (0.0, 0.0, 0.1, 0.0))AS v(a, b, c, d)

从三者的结果来看，精度大抵相当，没有太大差距，在极点的计算也都没有太大的偏差。经阅读源代码，可以看到h3c_dist函数使用了Haversine公式，q3c_dist函数使用了Haversine的一种变体，pgsphere的情况则较为复杂，它在距离大时使用大圆公式，距离较小时使用直角坐标公式。从计算结果来看，这3个插件都不失为成熟可用的扩展，具体要使用哪个，还要看其他的需求。如Q3C, H3C均提供了join函数(分别是q3c_join, h3c_join)，对两星表交叉证认进行了优化。而pgsphere提供了丰富的球面几何计算，如球面形状的面积计算，球面形状的交、并计算等。

5 总结

本文探讨了大圆公式、简化的大圆公式、Haversine公式等球面距离计算方法，并对比了它们在计算机中进行单精度及双精度数值计算的结果。在此基础上，引申出了基于三维直角坐标系的距离计算方法。结果表明大圆公式在单精度下舍入误差非常严重，而大圆公式的简化公式在两极或两点距离较大时误差很大。Haversine公式与直角坐标系方法的精度相近，都非常高，而直角坐标系方法的计算量较小，适合大数据的计算。此外，还需要注意所有这些公式在计算时，有时候计算结果会略微大于准确值。因而在实际应用时，可能需要将边界值取得略大一点，以将因为舍入误差漏掉的源包含在内。

已有软件包如Astropy, HTM, SLALIB, SOFA以及部分数据库扩展也将球面距离计算模块包含在内。它们所使用的计算方法不尽相同，结果精度也非常高。当需要进行球面距离计算时，可以直接使用这些库，而无须自行实现过于复杂的公式，以免重复工作且易因考虑不周而导致出错。