橢球點空間直角坐標轉換的研究

王仲鋒,劉 凱,初鳳婷

(1.長春工程學院 勘查與測繪工程學院，吉林 長春 130021;2.中水東北勘測設計研究有限責任公司,吉林 長春 130021)

目前，測量上用于地面點坐標轉換的模型主要有平面四參數模型、布爾沙(Bursa)模型、三維七參數模型、二維七參數模型和多元回歸模型等[1-2]。平面四參數模型是一種平面相似變換模型，當重合點所在坐標系的橢球長半軸和扁率差異不大、轉換坐標的區域不大且不存在跨帶問題時使用。

Bursa模型也是一種相似變換模型，需要知道重合點在源坐標系和目標坐標系中的大地高。但是目前使用的1954年北京坐標系、1980西安坐標系和2000國家大地坐標系中，只有2000國家大地坐標系中的點可以精確獲得大地高。王解先研究認為，對于地面上100 km×100 km的范圍，即使公共點中地方坐標的高程存在誤差，在求得轉換參數后，轉換出來的點平面坐標仍基本不變，當多個公共點地方坐標高程有誤差時，該結論仍然成立[3]。張銳研究認為：在應用“大地高不確定性主要表現在對高程的影響，而對平面坐標影響較小”這一結論時還是有一定的局限性[4]。張勤和王利研究認為，高程異常誤差對坐標的影響隨經度減小而增大；區域平均高程面的高度對平面坐標轉換精度隨高度增大而降低[5]。針對三維空間和高程未知的二維空間之間的七參數的求解問題，謝鳴宇、姚宜斌提出高程趨近法[6],甄佳寧、吳瓊、康靖玉、楊國東提出基于曲面逼近的空間坐標轉換方法[7]。

高程趨近法和曲面逼近法的思路基本相似，前者取公共點源坐標的大地高初值與目標坐標的初值一致，后者取源坐標與目標坐標的大地高初值均為0。高程趨近法和曲面逼近法均認為通過迭代可求解高精度參數，但在多次迭代后發現，每次迭代后所求參數并無明顯變化。劉根友、朱耀仲、朱才連為解決三維空間和大地高未知(但高程已知)的二維空間之間的七參數的求解問題，提出用高程加高程異常代替大地高，帶入Bursa模型，將高程異常作為未知參數和七參數一并求解的方法[8]。該方法忽略將高程異常作為未知參數和七參數一并求解帶來的另一個問題，即病態方程問題，使得解算結果極不穩定、極不可靠，有時可能無解。

二維七參數模型是在大地微分公式(三維七參數模型)基礎上，舍棄其中的大地高和大地高微分式而得來的一種用于不同坐標系下經緯度之間轉換的模型[1-2]。該理論極其不嚴密，存在用較大的增量代替微分的問題(相關研究成果將專文發表)。在實踐上，該公式也較為復雜，誤差方程呈病態現象也很嚴重。

多元回歸模型分為以平面坐標為變量的多元回歸模型和以經緯度為變量的多元回歸模型，前者適用于小區域不跨帶的平面坐標轉換，后者適用于省級及全國區域的坐標轉換問題。從實踐來看，多元回歸坐標轉換模型是高精度的轉換模型，但并未較好地被接受。

本研究將通過推導，證明不同坐標系下的同名橢球點之三維空間直角坐標之間，利用“布爾沙”模型進行轉換，并敘述轉換的方法步驟及應用實例。

1 轉換原理

由地面點的經度L、緯度B和大地高H可算得該點在三維空間直角坐標系中的坐標X，Y，Z，即:

(1)

當H=0時，即為地面點沿著橢球法線方向在橢球面上的投影點(簡稱橢球點)的空間直角坐標，簡稱橢球點空間直角坐標，記為標X0，Y0，Z0，即有：

(2)

其中：

(3)

e2=(a2-b2)/a2.

(4)

根據空間直角坐標轉換的“布爾沙”模型，點在空間直角坐標系“1”中的空間直角坐標通過3個平移參數ΔX0，ΔY0，ΔZ0，3個旋轉參數ωX，ωY，ωZ和1個尺度比參數m可轉換為空間直角坐標系“2”中的空間直角坐標，即有：

(5)

或：

(6)

利用n個公共點的三維空間直角坐標公式(6)，列出3n個誤差方程，即可利用最小二乘法解出7個坐標轉換參數。但是，當點在空間直角坐標系“1”和“2”中的大地高均未知時，通常認為“布爾沙”模型不再適用，因此有學者推導出了所謂“二維七參數模型”，用于同名點不同坐標系之間經緯度轉換。其中，“二維七參數模型”在理論上是不嚴密的，一是存在用較大的增量代替微分的問題，二是存在舍棄大地高的問題。

以下討論當點在空間直角坐標系“1”和“2”中的大地高均未知時如何利用“布爾沙”模型求解坐標轉換參數，實現同名點在不同坐標系之間經緯度的轉換問題。在以下的推導中，假定同名點在“1”和“2”坐標系中的經度和緯度差極小(實際上均在10秒以下)。

由式(1)，(2)知：

(7)

記：

(8)

或

(9)

則：

(ΔX12)2+(ΔY12)2+(ΔZ12)2=(H2cosB2cosL2-H1cosB1cosL1)2+

(H2cosB2sinL2-H1cosB1sinL1)2+(H2sinB2-H1sinB1)2=

而：

(cosB2cosL2cosB1cosL1+cosB2sinL2cosB1sinL1+sinB2sinB1)=

cosB2cosB1(cosL2cosL1+sinL2sinL1)+sinB2sinB1=

cosB2cosB1cos(L2-L1)+sinB2sinB1≈cosB2cosB1+sinB2sinB1=

cos(B2-B1)≈1.

故有：

(ΔX12)2+(ΔY12)2+(ΔZ12)2=(H2-H1)2.

(10)

取H1=H2=0時，式(6)可改寫為：

(11)

顯然，利用n個公共點在“1”、“2”空間直角坐標系中于橢球面上的投影點的三維空間直角坐標和式(11)，列出3n個方程，如果不考慮觀測誤差的話，利用最小原理求解的7個坐標轉換參數為滿足[ΔXΔX+ΔYΔY+ΔZΔZ]=min=0的解。實際上，由于H1=H2=0時已經使得ΔX12=ΔY12=ΔZ12=0，故利用最小二乘原理求解的7個坐標轉換參數實際上就是滿足誤差平方和等于最小的解。

由此可見，同名點在不同空間直角坐標系下沿法線方向投影到橢球面上的點的空間直角坐標之間的轉換，同樣可使用“布爾沙”模型。

當H1=H2≠0時，如果不考慮觀測誤差的話，利用最小二乘原理求解的7個坐標轉換參數為滿足[ΔXΔX+ΔYΔY+ΔZΔZ]≈0的解。

由式(2)知，橢球點的空間直角坐標，求出點的經度和緯度，即：

(12)

因此，同名點在不同空間直角坐標系下，沿法線方向投影到橢球面上的點的空間直角坐標之間的轉換，實際上也是同名點在不同坐標系下的經緯度之間的轉換。

2 轉換參數計算步驟

1)選取n個公共點，要求點位均勻分布地覆蓋測區，且點數不少于6個；

2)計算公共點在“1”，“2”空間直角坐標系中的橢球點坐標，即：

(i=1,2,…,n).

3)列誤差方程：

(12)

(13)

解式(13)，便可得要解算的7個參數。

(14)

3 算 例

如圖1所示，香港地區74個點(覆蓋面積月1 000 km2)，其在HK80坐標系和WGS84坐標系下的橢球點三維空間直角坐標見表1(根據香港地政署測繪處網站下載資料計算與整理)，試求HK80坐標系下的橢球點三維空間直角坐標向和WGS84坐標系下的橢球點三維空間直角坐標轉換的7個參數。

圖1 香港地區74個控制點分布示意圖

表1 香港地區74個控制點的橢球點之三維空間直角坐標(節選)

解：為了證明本文所述方法的可行性，取表中63-74號點作為建模點(如圖1中“口”形符號所示)，1-62號點作為檢驗點(如圖1中“◇”形符號所示)。

1)解算結果。利用中心化坐標求解，得3個平移參數-157.095 6 m、-267.509 8 m和-187.720 0 m；3個旋轉參數分別為2.863 667 s、7.079 517 s、-6.658 663 s；尺度比為-0.604 489 ppm。

2)內符合精度。反應建模內符合精度的點位殘差見圖2。其中空間點位差最小者為70號點，差值為0.10 cm，最大者為64號點，差值為0.64 cm，平均為0.38 cm。

圖2 內符合檢驗空間點位差散點

3)外符合精度。利用61個檢驗點進行外符合檢驗，相應的空間點位差分布情況見圖3。其中空間點位差最小者為49號點，差值為0.07 cm，最大者為53號點，差值為1.26 cm，平均為0.46 cm。

圖3 外符合檢驗空間點位差散點

4 結 論

本文推導不同坐標系下同名橢球點的三維坐標之間相互轉換的原理，敘述轉換的步驟和方法，并用實例驗證了其可行性。通過研究分析可得如下結論：

1)用“布爾沙”模型進行橢球點的三維坐標相互轉換理論嚴密、方法可行，適用于任何大小的區域，可代替“二維七參數模型”通過橢球點的三維空間直角坐標轉換來間接實現不同坐標系下地面點的經緯度轉換；

2)取同名點在“1”、“2”坐標系下的大地高H1=H2時(如H1代表北京54或西安80坐標系下的大地高(通常為未知值)，H2代表CGCS2000坐標系下的大地高(為已知值)，可取H54 or 80=H2 000)，其三維空間直角坐標也可用“布爾沙”模型進行相互轉換，但如此轉換帶有一點“近似”的性質。