球面距離計算方法及精度比較*

樊東衛，何勃亮，李長華，韓 軍，許允飛，崔辰州

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

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

兩點球面距離計算是天文或地理學中最常用的計算之一，是目標查找、錐形檢索、交叉證認等計算中最基礎的運算。比如查找距離某個位置最近的山峰，需要先對附近的山峰都計算一遍距離，取出距離最小的一個。而錐形檢索是虛擬天文臺的數據檢索協議之一，用于在球面上查找與某點的距離小于指定半徑的目標列表，幾乎所有提供了星表查詢功能的系統都有類似功能[1]。基于位置的星表交叉證認則可視為對一個星表A中所有的源做一遍錐形檢索取得與其最近的另一星表B的源，這樣就完成了A，B兩個星表的交叉[2]。這些應用的基礎都是球面距離計算。

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

1 球面幾何方法

球面距離計算可直接通過球面幾何推導。較常用的計算公式有大圓(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被應用到Astropy等代碼庫中，但其復雜度過高，本文在第3.1節中簡要提及。

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

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

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

由表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。

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

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

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

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

3 計算庫

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

3.1 Astropy

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

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

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

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

3.2 分層三角網格

分層三角網格(Hierarchical Triangular Mesh，HTM)[9]是一種對球面的多層次、遞歸三角分割方法。分層三角網格目前公開提供下載的軟件包⑦http：//www.skyserver.org/htm/index.html中包含了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的轉換。

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

表4 C#環境下Haversine與HTM fDistanceEq計算結果對比(單位：度)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實現，其目標是讓天文工作者更簡便快捷地編寫精確可靠的天體測量程序。雖然它是由Fortran 77實現的，但是已經有人為其編寫了Python接口pyslalib⑨https：//github.com/scottransom/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進一步將其轉換為法線向量⑩https：//en.wikipedia.org/wiki/N-vector，即與，再使用(11)式[10]進行計算。使用向量的主要優勢是向量代數可取代部分三角函數的計算，其數值計算穩定性更好，能保持較好的精度，并且在邊界點(如南北天極、0～360度邊界)也無異常。

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

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

4 數據庫

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

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

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

4.1 PostgreSQL數據庫插件

數據庫中除了可以直接使用公式之外，也可以使用各種數據庫擴展插件進行相關計算，減少使用者編程的麻煩。 如 PostgreSQL 數據庫有 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是專為地理信息系統設計的，雖然也可以在天文上使用，但是需要將米、千米等單位轉換回弧度，太過繁瑣低效。而pgphere，Q3C[11]，H3C則是專門為天文檢索設計的，更為簡潔。下文演示了3種插件的使用方法，表7對這3個插件的結果進行了對比。需要注意的是Q3C和H3C直接使用度進行計算，而pgsphere需要先轉成弧度進行計算，再將結果轉回角度。

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

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

5 總 結

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

已有軟件包如Astropy，HTM，SLALIB，SOFA以及部分數據庫擴展也將球面距離計算模塊包含在內。它們所使用的計算方法不盡相同，結果精度也非常高。當需要進行球面距離計算時，可以直接使用這些庫，而無須自行實現過于復雜的公式，以免重復工作且易因考慮不周而導致出錯。