在閱讀 RTKLIB的源碼時,發現了ECEF和大地坐標系的相互轉換的函數,大地坐標系(φ,λ,h)轉成ECEF(X,Y,Z)與所看書籍(GPS原理與接收機,謝剛,電子工業出版社)的公式是一樣的,而ECEF轉成大地坐標系的公式則與上述書籍和 RTKLIB的使用手冊(P135)中的公式都不一樣。簡而言之,遇到的問題如下圖所示。
一開始想着能不能直接在某個文獻中直接找到與程序算法一致的處理流程,但找了很久都沒能如願。最后,只能想着徹底從平面幾何上理解這個數學推導過程,從而再嘗試理解為什么程序會那么寫。
在查找介紹大地坐標系的文獻中,總能看到酉卯圓曲率半徑這個名詞。百度百科上給的定義是:過橢球面上一點的法線,可作無限個法截面,其中一個與該點子午面相垂直的法截面同橢球面相截形成的閉合的圈稱為卯酉圈,如下圖中的PEE′所示。
圖 1
投影到二維平面中,就是下面這幅圖。其中,Pn就稱為酉卯圓曲率半徑,PT則為橢圓的切線,角B則是酉卯圓曲率半徑與長軸的夾角,稱為緯度(也就是說通常意義上某個點的緯度並不是該點到地心的連線與長軸的夾角!)。
圖 2
由橢圓切線的斜率公式,可知
(1)
聯立上面兩式,可得
(2)
其中,。再將上式帶入到橢圓的標准方程
中,即可解得
(3)
另外,在△PnA中通過簡單的幾何關系可知,。所以有
(4)
又因為,所以有
(5)
注意,上面這個關於PQ、Qn的公式在后面推導ECEF和大地坐標系的相互轉換公式時,是很重要的。
為了說明坐標轉換的一般性,這里的待轉換點P是在地球外部的,示意圖如下。圖 3在理解(φ,λ,h)→XYZ中的X、Y坐標的變換公式時較為重要。
圖 3
而要想真正簡單、清晰地理解大地坐標系(φ,λ,h)與ECEF(X,Y,Z)如何轉換,還是要從二維平面圖中着手。下圖可以看作是橢球體的主視圖,XY平面壓縮成了水平X軸。
圖 4
- 大地坐標系(φ,λ,h)轉成ECEF(X,Y,Z)。
OD是點P在XY平面上的投影向量,其長度為
(6)
再結合圖3中XY平面中經度角λ的示意圖,可得
(7)
而關於z坐標,則有
(8)
- ECEF(X,Y,Z)轉成大地坐標系(φ,λ,h)。
首先是經度角λ,在圖3的XY平面中可以很清晰地看到
(9)
對於緯度角φ和高度h,可以在在圖4中分析其幾何關系。
(10)
其中,
(11)
從上述公式中可以看到,在計算φ時還是會用到緯度角φ的,所以不能直接用上述公式來計算緯度角。這里給出RTKLIB中ecf2pos函數中的算法步驟:
① 假設PD=PE,計算出夾角φ’。
② 根據式(4)和式(11)再計算N’和PE’。
③ 對比②中的PE’與①中的PE之間的差別是否小於截斷因子。不符合條件,將②中的PE’代入到①中,繼續循環計算;符合條件,則說明此時φ’=φ,可以跳出循環。
④ 將最終的PE代入到式(10)中,從而計算出緯度角φ和高度h。
至此,關於RTKLIB的源碼中的pos2ecf和ecf2pos函數的計算過程的理論公式推導過程,就證明完了。
這個事情證明,不會推公式的算法程序員是不合格的。找到的現成理論依據在復雜編程問題時可能會不夠高效,或者當算法流程稍作改變時就會看不懂別人寫的程序了。