﻿ 球面距离计算方法及精度比较

Research on Spherical Distance Computation and Accuracy Comparison
Fan Dongwei, He Boliang, Li Changhua, Han Jun, Xu Yunfei, Cui Chenzhou
National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
Key words: Numerical computation    Spherical distance    Angular separation    Calculation accuracy

1 球面几何方法

 $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)

 p1 p2 简化的大圆公式 大圆公式 Haversine公式 0, 0 0.01, 0 0.009999998845160010.01000000000000000 0.000000000000000000.00999999998265327 0.009999998845160010.01000000000000000 0, 0 0, 0.01 0.009999998845160010.01000000000000000 0.000000000000000000.00999999998265327 0.009999998845160010.01000000000000000 180, 0 180.01, 0 0.009985735639929770.00999999999999938 0.000000000000000000.00999999998265327 0.009985735639929770.00999999999999938 42, 43 42.01, 43 0.007310598157346250.00731353701619125 0.000000000000000000.00731353703719986 0.007310598157346250.00731353701187370 42, 43 42, 43.01 0.009999396279454230.00999999999999302 0.000000000000000000.01000000001909975 0.009999396279454230.00999999999999302 0, 90 180, 89.99 0.018634965643286710.01862095887437574 0.000000000000000000.00999999998265327 0.010001217015087600.01000000000000639 0, 0 0.1, 0 0.099999994039535520.10000000000000001 0.100870549678802490.10000000000019954 0.09999999403953552 0.10000000000000001

2 直角坐标系方法

 $x=\text{cos}\delta \text{cos}\alpha ,$ (4)
 $y=\text{cos}\delta \text{sin}\alpha ,$ (5)
 $z=\text{sin}\delta \text{ }.$ (6)

 \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)

 p1 p2 Haversine公式 直角坐标系方法 0, 0 0.01, 0 0.009999998845160010.01000000000000000 0.009999998845160010.01000000000000000 0, 0 0, 0.01 0.009999998845160010.01000000000000000 0.009999998845160010.01000000000000000 180, 0 180.01, 0 0.009985735639929770.00999999999999938 0.009985735639929770.00999999999999938 42, 43 42.01, 43 0.007310598157346250.00731353701187370 0.007312425877898930.00731353701187507 42, 43 42, 43.01 0.009999396279454230.00999999999999302 0.009999873116612430.00999999999999879 0, 90 180, 89.99 0.010001217015087600.01000000000000639 0.010001217015087600.01000000000000639 0, 0 0.1, 0 0.099999994039535520.10000000000000001 0.09999999403953552 0.10000000000000001

 ${{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 计算库

3.1 Astropy

d = c1.separation(c2).deg

 p1 p2 Haversine Astropy 0, 0 0.01, 0 0.01 0.01 0, 0 0, 0.01 0.01 0.01 180, 0 180.01, 0 0.009999999999999377 0.009999999999990905 42, 43 42.01, 43 0.007313537011873697 0.007313537011872698 42, 43 42, 43.01 0.009999999999993016 0.009999999999994572 0, 90 180, 89.99 0.010000000000006394 0.010000000000006394 0, 0 0.1, 0 0.1 0.1

 $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 分层三角网格

 p1 p2 Haversine HTM fDistanceEq 0, 0 0.01, 0 0.01 0.01 0, 0 0, 0.01 0.01 0.01 180, 0 180.01, 0 0.00999999999999938 0.00999999999999938 42, 43 42.01, 43 0.0073135370118737 0.00731353701187507 42, 43 42, 43.01 0.00999999999999302 0.00999999999999879 0, 90 180, 89.99 0.0100000000000064 0.0100000000000064 0, 0 0.1, 0 0.1 0.1

3.3 SLALIB与SOFA

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

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度边界)也无异常。

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

 p1 p2 Haversine pyslalib SOFA 0, 0 0.01, 0 0.01 0.01 0.01000000000000000 0, 0 0, 0.01 0.01 0.01 0.01000000000000000 180, 0 180.01, 0 0.009999999999999377 0.009999999999999377 0.00999999999999938 42, 43 42.01, 43 0.007313537011873697 0.007313537011875421 0.00731353701187542 42, 43 42, 43.01 0.009999999999993016 0.0099999999999971 0.00999999999999710 0, 90 180, 89.99 0.010000000000006394 0.010000000000006394 0.01000000000000639 0, 0 0.1, 0 0.1 0.1 0.10000000000000001

4 数据库

 p1 p2 MS SQL Server PostgreSQL MySQL 0, 0 0.01, 0 0.00999999999999926 0.01 0.01 0, 0 0, 0.01 0.00999999999999926 0.01 0.01 180, 0 180.01, 0 0.00999999999999811 0.00999999999999938 0.009999999999999377 42, 43 42.01, 43 0.00731353701187193 0.0073135370118737 0.007313537011873697 42, 43 42, 43.01 0.00999999999999697 0.00999999999999302 0.009999999999993016 0, 90 180, 89.99 0.0100000000000109 0.0100000000000064 0.010000000000006394 0, 0 0.1, 0 0.0999999999999995 0.1 0.1

4.1 PostgreSQL数据库插件

 p1 p2 Q3C H3C pgsphere 0.0, 0.0 0.01, 0.0 0.01 0.01 0.00999999998265327 0.0, 0.0 0.0, 0.01 0.01 0.01 0.00999999998265327 180.0, 0.0 180.01, 0.0 0.00999999999999091 0.00999999999999938 0.00999999998265327 42.0, 43.0 42.01, 43.0 0.0073135370118727 0.0073135370118737 0.00731353703719986 42.0, 43.0 42.0, 43.01 0.00999999999999801 0.00999999999999302 0.0100000000190997 0.0, 90.0 180.0, 89.99 0.0100000000000064 0.0100000000000064 0.00999999998265327 0.0, 0.0 0.1, 0.0 0.1 0.1 0.1000000000002

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

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

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

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)

5 总结

 [1] 张海龙, 冶鑫晨, 李慧娟, 等. 天文数据检索与发布综述[J]. 天文研究与技术, 2017, 14(2): 212–228 [2] 张海龙, 聂俊, 赵青, 等. 新疆天文台在线交叉证认服务[J]. 天文研究与技术, 2017, 14(3): 347–355 [3] 张彦霞.多波段天体物理中的自动分类方法研究[D].北京: 中国科学院国家天文台, 2003. http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y557147 [4] 高丹.海量天文数据融合系统的开发与数据挖掘算法的研究[D].北京: 中国科学院国家天文台, 2008. [5] 赵青.面向海量数据的高效天文交叉证认的研究[D].天津: 天津大学, 2010. http://cdmd.cnki.com.cn/Article/CDMD-10056-1011266617.htm [6] 杜鹏.海量天文星表数据的交叉与融合[D].济南: 山东学, 2013. http://cdmd.cnki.com.cn/Article/CDMD-10422-1013222095.htm [7] ROBITAILLE T P, TOLLERUD E J, GREENFIELD P, et al. Astropy:a community Python package for astronomy[J]. Astronomy & Astrophysics, 2013, 558:A33(9pp). [8] VINCENTY T. Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations[J]. Empire Survey Review, 1975, 23(176): 88–93. DOI: 10.1179/sre.1975.23.176.88 [9] KUNSZT P Z, SZALAY A S, THAKAR A R. The Hierarchical Triangular Mesh[C]//Proceedings of the MPA/ESO/MPE Workshop. 2001: 631-637. [10] KENNETH G. A non-singular horizontal position representation[J]. The Journal of Navigation, 2010, 63(3): 395–417. DOI: 10.1017/S0373463309990415 [11] KOPOSOV S, BARTUNOV O. Q3C, Quad Tree Cube-the new sky-indexing concept for huge astronomical catalogues and its realization for main astronomical queries (cone search and Xmatch) in open source database PostgreSQL[C]//Astronomical Data Analysis Software and Systems XV ASP Conference Series. 2006: 735-738.

0

#### 文章信息

Fan Dongwei, He Boliang, Li Changhua, Han Jun, Xu Yunfei, Cui Chenzhou

Research on Spherical Distance Computation and Accuracy Comparison

Astronomical Research and Technology, 2019, 16(1): 69-76.