樊東衛(wèi)，何勃亮，李長華，韓 軍，許允飛，崔辰州

(中國科學(xué)院國家天文臺，北京 100101)

平面距離只需要計算連接兩點線段的長度，通過二維直角坐標(biāo)使用勾股定理可以很容易地計算出來。由于只涉及乘法、加法與開平方運算，在計算機(jī)中也可以保留很高的精度。而球面距離計算則要復(fù)雜得多。球面在直角坐標(biāo)系中實際上是三維結(jié)構(gòu)，除非距離相對于半徑非常小，否則不能當(dāng)作平面處理。

兩點球面距離計算是天文或地理學(xué)中最常用的計算之一，是目標(biāo)查找、錐形檢索、交叉證認(rèn)等計算中最基礎(chǔ)的運算。比如查找距離某個位置最近的山峰，需要先對附近的山峰都計算一遍距離，取出距離最小的一個。而錐形檢索是虛擬天文臺的數(shù)據(jù)檢索協(xié)議之一，用于在球面上查找與某點的距離小于指定半徑的目標(biāo)列表，幾乎所有提供了星表查詢功能的系統(tǒng)都有類似功能[1]?；谖恢玫男潜斫徊孀C認(rèn)則可視為對一個星表A中所有的源做一遍錐形檢索取得與其最近的另一星表B的源，這樣就完成了A，B兩個星表的交叉[2]。這些應(yīng)用的基礎(chǔ)都是球面距離計算。

兩點球面距離計算，計算的是兩點在球面上的最短距離，即球心與兩點組成的大圓上兩點之間較小的那段圓弧的長度，天文上習(xí)慣使用該短圓弧對應(yīng)的大圓圓心角，稱為角間距(Angular Separation)或角距離。球面距離需要使用球面幾何進(jìn)行求解，要大量使用三角函數(shù)。而計算機(jī)的數(shù)值計算精度有限，在對三角函數(shù)、反三角函數(shù)進(jìn)行計算時容易出現(xiàn)舍入誤差。多個函數(shù)的誤差積累之后，將導(dǎo)致嚴(yán)重的結(jié)果偏離，甚至無法得到正確的結(jié)果。前人已經(jīng)推導(dǎo)了很多數(shù)學(xué)公式求解球面距離，但在何種情況下使用哪個公式，公式的精度如何，卻鮮有人討論。因此，本文考察天文技術(shù)中常用的球面距離計算公式，給出它們的算法，并對比它們的精度，為球面距離算法的選擇提供參考。

1 球面幾何方法

球面距離計算可直接通過球面幾何推導(dǎo)。較常用的計算公式有大圓(Great-circle)公式①https：//en.wikipedia.org/wiki/Great-circle＿distance、Haversine公式②https：//en.wikipedia.org/wiki/Haversine＿formula等。另外還有Vincenty公式③https：//en.wikipedia.org/wiki/Vincenty%27s＿formulae被應(yīng)用到Astropy等代碼庫中，但其復(fù)雜度過高，本文在第3.1節(jié)中簡要提及。

設(shè)有兩點赤道坐標(biāo)為p1(α1，δ1)，p2(α2，δ2)，求它們的球面角距離d。其中大圓公式為(1)式；當(dāng)兩點距離很小時，也有人使用簡化的(2)式[2-6]；Haversine公式如(3)式。

計算機(jī)中通常使用IEEE754二進(jìn)制浮點數(shù)算術(shù)標(biāo)準(zhǔn)處理浮點數(shù)，其中單精度32位(float或single)真值可精確到小數(shù)點后6～9位，雙精度64位(double)真值可精確到小數(shù)點后15～17位④https：//en.wikipedia.org/wiki/Single-precision＿floating-point＿format。本文使用C＋＋語言編寫程序?qū)?個公式的精度進(jìn)行了對比，結(jié)果如表1。本文使用的單位均為度(°)，表中有兩行數(shù)據(jù)時，上方的值為單精度計算的結(jié)果，下方的值為雙精度計算的結(jié)果。

表1 使用不同的點測試常用球面距離公式的精度(單位：度)Table 1 Accuracy testing at different points for widely used formulas(unit： degree)

由表1可見，單精度下，兩點之間距離較近時，大圓公式返回的是全0，出現(xiàn)了嚴(yán)重的舍入誤差，而距離較大時誤差也較大。雙精度下，大圓公式能夠返回結(jié)果，但比Haversine公式的精度要差。大圓公式的簡化形式多數(shù)時候與Haversine公式的精度相當(dāng)。這3個公式共有的一個問題，它們均有可能略微超過準(zhǔn)確值，當(dāng)嚴(yán)格限定一個邊界值的時候，可能導(dǎo)致漏源。因而在實際應(yīng)用時，需要考慮將邊界值略微外擴(kuò)，以將因為計算精度而漏掉的源包含在內(nèi)。

值得注意的是，大圓公式的簡化形式在極點附近時誤差非常大。雖然它在上述公式中是計算復(fù)雜度最低的，并且在兩極之外、距離較小時精度與Haversine公式相近，但是在實際使用該公式時仍應(yīng)當(dāng)非常小心。

2 直角坐標(biāo)系方法

在基于赤道坐標(biāo)系進(jìn)行求解之外，也可以使用直角坐標(biāo)系，通過三維坐標(biāo)向量對兩點球面距離進(jìn)行計算?？梢曁烨虬霃綖閱挝?，以天球球心為三維直角坐標(biāo)系原點，赤道坐標(biāo)系北天極點方向為z軸方向，天赤道面上赤經(jīng)為0點的方向為x正向，遵循右手坐標(biāo)系與x-z相交的軸為y軸方向。從而可將赤經(jīng)、赤緯轉(zhuǎn)換為三維直角坐標(biāo)(x，y，z)，其轉(zhuǎn)換關(guān)系如(4～6)式⑤https：//arxiv.org/abs/cs/0701171。

這樣兩點距離計算可以直接使用(8)式進(jìn)行計算。該公式可簡單地從Haversine公式中推得，其中Haversine公式根號中的部分(7式)實際上是兩點的三維空間直線距離一半的平方。

直角坐標(biāo)系的計算精度方面，本文依舊使用C＋＋語言在單精度及雙精度兩種環(huán)境下對公式進(jìn)行計算，并與Haversine公式的結(jié)果對比，結(jié)果如表2。從表2可以看出，直角坐標(biāo)系方法的精度與Haversine公式幾乎完全一致，而直角坐標(biāo)系方法有時甚至優(yōu)于Haversine公式(見42，43兩行數(shù)據(jù))。

表2 Haversine與直角坐標(biāo)系計算結(jié)果對比(單位：度)Table 2 Result comparison between Haversine and Cartesian method(unit： degree)

直角坐標(biāo)系方法在進(jìn)行距離比較時還可以進(jìn)一步優(yōu)化，如(9)式：設(shè)兩點直角坐標(biāo)為p1(x1，y1，z1)，p2(x2，y2，z2)，其優(yōu)點是可以預(yù)先計算x，y，z與threshold(設(shè)最大角距離為θ)，從而避免在實際進(jìn)行距離計算時使用非常耗時的三角函數(shù)，而只需要三次減法、三次乘法、兩次加法、一次比較，大大減少了計算量，非常適用于大規(guī)模數(shù)據(jù)計算。相應(yīng)地，它的缺點是需要額外占用存儲空間來保存x，y，z三個值，在對內(nèi)存占用要求較苛刻的環(huán)境中可能會成為問題。因而，實際的公式選擇，仍應(yīng)視數(shù)據(jù)規(guī)模、計算環(huán)境而定。另外，基于三維向量的計算方法，還有法線向量的形式，但其計算過程比本節(jié)要復(fù)雜，本文將在第3.3節(jié)中進(jìn)行簡要描述。

3 計算庫

許多天文軟件包也包含了球面距離計算功能。本部分選擇了幾個安裝、使用較便捷，應(yīng)用也比較廣泛的庫，演示了其使用方法，并進(jìn)行測試。由于難以分辨軟件包中實際使用的是單精度還是雙精度，此部分不再對此進(jìn)行區(qū)分，而直接與Haversine雙精度結(jié)果進(jìn)行比較。

3.1 Astropy

當(dāng)前天文界已經(jīng)非常流行使用Python進(jìn)行數(shù)據(jù)處理，其中Astropy[7]對Python的推廣傳播功不可沒。Astropy進(jìn)行球面距離計算需要使用astropy.units及astropy.coordinates.SkyCoord計算。調(diào)用方法如下，其中d為兩點之間的角距離(單位是度)：

為公平起見，在64位Windows10中的64位python3.6.3重新實現(xiàn)了Haversine算法，并與Astropy的距離計算函數(shù)進(jìn)行了對比，對比結(jié)果見表3，可以看到兩種方法基本一致，不分伯仲。

表3 Windows 64位Python3環(huán)境下Haversine與Astropy計算結(jié)果對比(單位：度)Table 3 Result comparison between Haversine and Astropy on Windows 64bit Python3(unit： degree)

閱讀Astropy的源代碼可發(fā)現(xiàn)其距離計算函數(shù)separation實際調(diào)用的是astropy/coordinates/angle＿utilities.py⑥https：//github.com/astropy/astropy/blob/master/astropy/coordinates/angle＿utilities.py中的angular＿separation(lon1，lat1，lon2，lat2)函數(shù)，其代碼實現(xiàn)使用Vincenty公式，即(10)式。該式在所有位置的距離計算都更穩(wěn)定，但計算復(fù)雜度也更高[8]。從上述計算對比結(jié)果來看，Astropy并沒有處處比Haversine的精度更高，如(180，0)一行；當(dāng)然，Astropy在(42，43)上的精度也比Haversine更好。但是兩者的差距都非常小，幾乎可以忽略不計。從比較結(jié)果上看，Haversine已經(jīng)相當(dāng)精確。當(dāng)然，若是對計算穩(wěn)定性有更高要求，可以考慮Vincenty公式。

3.2 分層三角網(wǎng)格

分層三角網(wǎng)格(Hierarchical Triangular Mesh，HTM)[9]是一種對球面的多層次、遞歸三角分割方法。分層三角網(wǎng)格目前公開提供下載的軟件包⑦h(yuǎn)ttp：//www.skyserver.org/htm/index.html中包含了C＃代碼庫及一個SQL Server擴(kuò)展包。其中Spherical.Htm.Sql.cs文件中包含了函數(shù)fDistanceEq(ra1，dec1，ra2，dec2)，可用于計算距離。需要注意的是，fDistanceEq的返回值單位是arcmin，不是degree。針對分層三角網(wǎng)格，本文使用C＃實現(xiàn)了Haversine函數(shù)并與SphericalHTM v3.1.2的fDistanceEq函數(shù)進(jìn)行對比。如表4，其結(jié)果基本相同。查看源代碼發(fā)現(xiàn)fDistanceEq實際上也使用了Haversine公式進(jìn)行計算，只是多了degree到arcmin的轉(zhuǎn)換。

除了C＃程序中可以直接調(diào)用fDistanceEq外，HTM軟件包提供的SQL Server數(shù)據(jù)庫擴(kuò)展包中也帶有此函數(shù)，在數(shù)據(jù)庫中使用非常方便。

表4 C#環(huán)境下Haversine與HTM fDistanceEq計算結(jié)果對比(單位：度)Table 4 Result comparison between Haversine and HTM fDistanceEq on C#(unit： degree)

3.3 SLALIB與SOFA

SLALIB即Starlink library of positional astronomy routines⑧http：//stsdas.stsci.edu/cgi-bin/gethelp.cgi?slalib.sys，使用Fortran 77實現(xiàn)，其目標(biāo)是讓天文工作者更簡便快捷地編寫精確可靠的天體測量程序。雖然它是由Fortran 77實現(xiàn)的，但是已經(jīng)有人為其編寫了Python接口pyslalib⑨https：//github.com/scottransom/pyslalib，因而可以在已有的Python環(huán)境中獲得Fortran 77的計算精度。

Pyslalib中的slalib.sla＿dsep直接對應(yīng)了SLALIB中的sla＿DSEP(A1，B1，A2，B2)球面距離計算函數(shù)，在64位的Ubuntu 16.04.4操作系統(tǒng)中使用64位 Python3.5.2環(huán)境將其與Haversine函數(shù)的計算結(jié)果進(jìn)行了對比，如表5。首先可以看到Linux版的Python Harversine程序計算結(jié)果與Windows一致，表明Python3在不同的操作系統(tǒng)中的表現(xiàn)一致。其次pyslalib的計算結(jié)果，除了(42，43)兩行有輕微差別，與Haversine完全一致，精度也非常高。查看SLALIB的源代碼，可以看到其使用了第2節(jié)中的直角坐標(biāo)系三維向量作為計算單元。但與第2節(jié)不同，SLALIB進(jìn)一步將其轉(zhuǎn)換為法線向量⑩https：//en.wikipedia.org/wiki/N-vector，即與，再使用(11)式[10]進(jìn)行計算。使用向量的主要優(yōu)勢是向量代數(shù)可取代部分三角函數(shù)的計算，其數(shù)值計算穩(wěn)定性更好，能保持較好的精度，并且在邊界點(如南北天極、0～360度邊界)也無異常。

天體測量中常用的SOFA庫?http：//www.iausofa.org的球面計算函數(shù)iauSeps?http：//www.iausofa.org/2018＿0130＿C/SeparationPA.html也同樣使用了(11)式。SOFA即Standards of Fundamental Astronomy，是國際天文學(xué)聯(lián)合會建立并維護(hù)的一個權(quán)威的基本天文學(xué)標(biāo)準(zhǔn)庫，包含有ANSI C及Fotran 77兩個版本。表5的最右側(cè)列出了雙精度位C＋＋程序調(diào)用SOFA計算的結(jié)果。

表5 Linux 64位Python3環(huán)境下Haversine與pyslalib及雙精度C++版本SOFA計算結(jié)果對比(單位：度)Table 5 Result comparison among Haversine and pyslalib on Linux 64bit Python3 and double precision C++version SOFA(unit：degree)

4 數(shù)據(jù)庫

關(guān)系型數(shù)據(jù)庫均支持SQL語句，因而可以在不同的系統(tǒng)中直接使用SQL語句實現(xiàn)Haversine公式進(jìn)行距離計算。由于每種數(shù)據(jù)庫使用的數(shù)學(xué)函數(shù)、計算精度略有差別，本文在Microsoft SQL Server 2008，PostgreSQL 9.4.5，MySQL 5.7.18上進(jìn)行了測試，測試結(jié)果如表6，從表6可以看到，3個數(shù)據(jù)庫的差距不大，MySQL默認(rèn)保留的精度更多一些。

表6 不同數(shù)據(jù)庫中使用Haversine公式計算距離的精度對比(單位：度)Table 6 Accuracy testing for pure SQL Haversine on different database(unit： degree)

除了Haversine公式，數(shù)據(jù)庫中使用直角坐標(biāo)系方法實際上更方便。數(shù)據(jù)庫對I/O進(jìn)行了大量優(yōu)化，可以更好地對數(shù)據(jù)進(jìn)行調(diào)度。因而，直角坐標(biāo)系方法雖然需要增加x，y，z 3個字段，但并不會給數(shù)據(jù)庫增加太大的負(fù)擔(dān)，還可以大大降低計算的復(fù)雜度，更快地取得海量數(shù)據(jù)的計算結(jié)果。

4.1 PostgreSQL數(shù)據(jù)庫插件

數(shù)據(jù)庫中除了可以直接使用公式之外，也可以使用各種數(shù)據(jù)庫擴(kuò)展插件進(jìn)行相關(guān)計算，減少使用者編程的麻煩。 如 PostgreSQL 數(shù)據(jù)庫有 PostGIS?https：//postgis.net/， pgsphere?https：//github.com/akorotkov/pgsphere， Q3C?https：//github.com/segasai/q3c， H3C?http：//cds.u-strasbg.fr/resources/doku.php?id＝h3c等多個插件支持球面距離計算。其中PostGIS是專為地理信息系統(tǒng)設(shè)計的，雖然也可以在天文上使用，但是需要將米、千米等單位轉(zhuǎn)換回弧度，太過繁瑣低效。而pgphere，Q3C[11]，H3C則是專門為天文檢索設(shè)計的，更為簡潔。下文演示了3種插件的使用方法，表7對這3個插件的結(jié)果進(jìn)行了對比。需要注意的是Q3C和H3C直接使用度進(jìn)行計算，而pgsphere需要先轉(zhuǎn)成弧度進(jìn)行計算，再將結(jié)果轉(zhuǎn)回角度。

表7 PostgreSQL數(shù)據(jù)庫Q3C，H3C，pgsphere插件的球面距離計算精度比較(單位：度)Table 7 Accuracy testing for different PostgreSQL extensions(unit： degree)

從三者的結(jié)果來看，精度大抵相當(dāng)，沒有太大差距，在極點的計算也都沒有太大的偏差。經(jīng)閱讀源代碼，可以看到h3c＿dist函數(shù)使用了Haversine公式，q3c＿dist函數(shù)使用了Haversine的一種變體，pgsphere的情況則較為復(fù)雜，它在距離大時使用大圓公式，距離較小時使用直角坐標(biāo)公式。從計算結(jié)果來看，這3個插件都不失為成熟可用的擴(kuò)展，具體要使用哪個，還要看其他的需求。如Q3C，H3C均提供了join函數(shù)(分別是q3c＿join，h3c＿join)，對兩星表交叉證認(rèn)進(jìn)行了優(yōu)化。而pgsphere提供了豐富的球面幾何計算，如球面形狀的面積計算，球面形狀的交、并計算等。

5 總 結(jié)

本文探討了大圓公式、簡化的大圓公式、Haversine公式等球面距離計算方法，并對比了它們在計算機(jī)中進(jìn)行單精度及雙精度數(shù)值計算的結(jié)果。在此基礎(chǔ)上，引申出了基于三維直角坐標(biāo)系的距離計算方法。結(jié)果表明大圓公式在單精度下舍入誤差非常嚴(yán)重，而大圓公式的簡化公式在兩極或兩點距離較大時誤差很大。Haversine公式與直角坐標(biāo)系方法的精度相近，都非常高，而直角坐標(biāo)系方法的計算量較小，適合大數(shù)據(jù)的計算。此外，還需要注意所有這些公式在計算時，有時候計算結(jié)果會略微大于準(zhǔn)確值。因而在實際應(yīng)用時，可能需要將邊界值取得略大一點，以將因為舍入誤差漏掉的源包含在內(nèi)。

已有軟件包如Astropy，HTM，SLALIB，SOFA以及部分?jǐn)?shù)據(jù)庫擴(kuò)展也將球面距離計算模塊包含在內(nèi)。它們所使用的計算方法不盡相同，結(jié)果精度也非常高。當(dāng)需要進(jìn)行球面距離計算時，可以直接使用這些庫，而無須自行實現(xiàn)過于復(fù)雜的公式，以免重復(fù)工作且易因考慮不周而導(dǎo)致出錯。