單玉浩，楊曉東，邢世宏

(海軍潛艇學院，山東 青島 266000)

幾何大地測量采用一個同地球外形最為接近的旋轉橢球作為參考模型，研究地球橢球數學性質、橢球投影數學變換、橢球面測量計算等。幾何大地測量可為航海提供定位、定向、導航信息[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)約束，然后基于此約束提出球面距離的積分解算模型，并給出模型的解析解，分析奇異點(沿大圓弧的同一走向導致緯度變化量的符號發(fā)生改變的點)對模型求解的影響，進一步完善了解算模型。球面距離模型的解析解可應用于大地主題解算，同時可為航海計算及輔助決策等問題提供基礎支撐。

1 球面A·r約束

經過球心的平面截得的球面圓為球面大圓。球面大圓的圓弧、劣弧分別稱為球面大圓弧、大圓劣弧。通過球面上任意兩點(非直徑對點)有且僅有一個大圓。任意大圓上的兩點A0和B將大圓分為兩部分，都叫做以A、B點為端點的大圓弧。

(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φ異號。

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 球體大圓

(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)中第二部分積分在積分區(qū)間[uv,ub]上，du始終為負值，而被積函數式中為正，因此式(15)又可表示為：

(16)

由此得出結論：若大圓弧上的連續(xù)點的緯度值非單調，則式(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 經緯度變化區(qū)間內的大圓距離

圖6 經緯度變化區(qū)間內的大圓距離比對差值

通過研究在p2點的經緯度變化范圍內積分模型解算的球面距離如圖5所示，與球面三角公式的比對差值如圖6所示。第二種情況研究了高緯度下球面距離的準確性問題，球面距離最大約為2.7×105m，隨著球面距離的增加，差值略有增大，但差值最大值控制在1×10-7m左右?？梢姶四Ｐ颓蠼馇蛎婢嚯x具有非常高的可靠性。

上述情況中不涉及奇異點問題，為研究奇異點存在情況下模型的解算準確性，選取了如表1所示的幾種典型奇異點存在情況，計算球面要素并進行大圓頂點的位置判斷，給出模型的解算誤差。表1中第一列是在兩點緯度相差非常小的情形下得到的距離求解結果，第二列是在高緯度下同緯度兩點的解算結果，第三列是兩點長距離情形下得到的解算結果，3種情形下距離解算誤差都控制在1×10-7m以內。經驗證在兩點同經度或同緯度及長距離的情況下積分解算模型依然有效。

表1 幾種典型奇異點存在情況下球面元素

5 結論

提出的基于A·r約束的大地距離解算模型能夠有效解算大地距離，且計算精度與球面三角公式相當。球面距離模型應用不受使用條件限制，在高緯度地區(qū)或同緯度、同經度等情況下能保證較高的精度。球面距離模型的解析解可應用于大地主題解算，同時可為航海計算及輔助決策等問題提供重要支撐。