坐標轉換 張嘉強 1 GPS 測量基準 -WGS84 衛星測量之基準由於著重全球性大地參考系統之觀點, 而與一般之區域性大地基準在定義上有所不同 一般而言, 其提供之資訊包括地球之幾何形狀及基本參考框架 地球之重力場模式以及不同區域性大地系統對應至地心固定坐標系統之平均相應位置 對 GPS 之全球性應用而言, 其作為基準所建立之地心系統稱為 WGS84(World Geodetic Syste, 1984) WGS84 坐標系統之原點位在地球質量中心, 其 Z 軸與 BIH(Bureau International de l Heure) 所定義之慣用地形北極 (Conventional Terrestrial Pole, CTP) 方向相平行,X 軸係與 BIH 所定義之零子午圈平面相平行之參考子午面在 CTP 赤道平面之交線方向,Y 軸則係與 X 軸垂直而在赤道面上所形成之右旋地球固定直角坐標系統 ( 如圖 1) BIH 所定義之 CTP(1984.0) BIH 所定義之零子午圈 Z 軸 地球質量中心 X 軸 Y 軸 圖 1 WGS-84 坐標系統 1
由於橢球體通常會配合一組坐標系統予以引用, 以便將地球上之點位位置以慣用之地理或地心坐標之形式展現, 因此 WGS84 所採用之橢球體相關幾何與物理常數現摘述於表 1 之中 其中之參數值除了第二階帶諧係數之外, 大都與 GRS80 橢球體之相關係數相同 GRS80 是 IUGG 所採用之最新地心橢球體, 其與 WGS84 之間的微量差異僅在其幾何扁率上, 這是因為 WGS84 所採用之第二階帶諧係數 (C2,0) 在取用 GRS80 之地球動力因子 (J2) 作數學聯繫時只截取到 8 位數字之結果 表 1 WGS84 參數及常數 參數及常數 符號 數值 長半徑 a 6378137±2 地球引力常數 GM (3986005±0.6) 10 8 3 s -2 正規化之第二階帶諧 C2,0 (-484.16685±0.00130) 10-6 引力係數 地球自轉角速度 ω (7292115±0.1500) 10-11 rad s -1 真空中之光速 C 299792458±1.2 s -1 扁率 f 1 / 298.257223563 2 全球地心參考框架 -ITRF 近年來隨著空間定位技術的日益精進, 多項高精度定位之應用已需要一組定義良好之全球性坐標系統 假定不同空間定位技術所建立之各個資料組能予以綜合分析, 並將各組之坐標成果資料合併形成單一之資料組, 則此一包含測站坐標與地球定向參數之地形參考系統即可定義完成並視為一個大地應用所需之參考框架 目前最重要之一組參考框架係由 IERS(International Earth Rotation Service) 所提供之 ITRF 2
(IERS Terrestrial Reference Frae) IERS 這個組織是由 IAU(International Astronoical Union) 與 IUGG 在 1987 年所聯合建立的, 並在 1988 年初開始運作 在定義 ITRF 時是採用一組基本測站所組成之全球網, 該類測站之坐標值則必須採用最精確之空間定位技術予以良好之測定 在考慮點位會受地殼變動之影響而產生坐標改變之現象,ITRF 除透過一組地面測站坐標值之公佈外, 測站坐標變化之速度場量亦會加以提供 以新的台灣大地基準 (TWD97) 定義時所使用之一組 ITRF 系列為例, 其稱之為 ITRF94, 其中各測站分類成 A B C 三個等級, 且分別提供其在 1993.00 時刻之坐標及其相應之速度場 定義 ITRF94 之全球測站則係分別採用 VLBI SLR DORIS 及 GPS 所進行觀測之坐標成果資料 ITRF94 所依據之基準則係採用下述之定義 : 其原點是 SLR 與 GPS 成果資料之加權平均 ; 其尺度是 VLBI SLR 與 GPS 成果資料之加權平均 ; 其方位與 ITRF92 在 1988.0 之時刻相一致 ; 其隨時間之變量與 NNR-NUVEL1A 之地球物理模式相一致 而當前最新一組之 ITRF 系列稱之為 ITRF2000, 由該基準之定義可知, 其尺度是 VLBI 與 SLR 成果之加權平均, 原點則是 SLR 成果之加權平均, 而方位及其速率則是由 50 個高品質之測點 ( 至少 3 年之連續觀測 位於板塊穩固區且遠離變形帶 速度場量之誤差優於 3/y 至少三組成果間之速度場量殘差值小於 3/y) 所決定, 其中之方位與 ITRF97 在 1997.0 時刻之值相一致, 方位速率則與 NNR-NUVEL1A 相一致 ITRF2000 所提供高精度三維坐標值及其速度場量之全球測點分布可繪如圖 2 所示 3
ITRF 與 GPS 之互動關係自 1992 年 IGS(International GPS Service) 開始運作後變得日益密切 透過互相之合作,ITRF 可以利用 IGS 公佈之全球 GPS 測站之連續追蹤資料, 而 IGS 分析中心亦可採用 ITRF 所公佈之坐標進行 GPS 之軌道計算 現階段在利用 GPS 進行高精度之定位應用時, 由於 GPS 必須要能夠提供公分級之全球性測站坐標資料, 因此 ITRF 坐標參考框架之採用便顯示出其優異性與必然性 圖 2 ITRF2000 之主要測站 3 GPS 測定之高程高程的定義通常是由一參考面向上量測之線量距離, 在傳統之大地測量中, 代表高程之正高 (orthoetric height) 係採用水準測量, 從近似大地水準面之平均海水面 (ean sea level,msl) 起算而量得 以大地水準面為基準之正高, 因其可由大地位數 (geopotential nuber) 加以定義而具有物理上的意義, 即 H C= W0 W= gdh 0 (1) 式中 C 表大地位數, W 表地面點之重力位, W0 表大地水準面之重力位, g 表觀測之重力值, dh 表各段之高差值 4
而正高 (H) 之定義則為 H = C (2) g 其中 g 是大地水準面至地面點之間沿垂線 (plub line) 上之重力平均 值 近年來隨著 GPS 之發展, 以橢球面為基準之橢球高 (ellipsoidal height) 也常被加以應用, 這是因為 GPS 定位所決定之地面點三維直角 坐標, 可以轉換為參考橢球體 ( 如 WGS84) 上大地坐標之緣故 橢球高 是一種幾何高度, 它表示從參考之橢球面沿法線 (noral line) 量至地 面點之高度距離量 橢球高 (h) 與正高 (H) 之差稱為大地水準面高 (geoidal height) 或是 大地起伏 (geoidal undulation), 通常以 N 表示 ( 如圖 3), 兩者間之近 似關係為 H = h N (3) P Topography Geoid Ellipsoid h H N 圖 3 大地起伏之定義 以全球觀點來看, 相應於一個最密切參考橢球體上之大地起伏值約 為 ±100 之等級 ( 如圖 4) 對台灣地區而言, 大地起伏之平均值約為 22 ( 如圖 5), 換言之, 利用 GPS 所測得之幾何高可減去約 22, 5
即可概略轉換為以平均海水面起算之正高值, 惟對精度需求較高之高程 間轉換而言, 則需利用重力資料配合大地水準面模式加以計算獲得 ( 我 國內政部已提供較精密之大地起伏值內差計算程式 ) 圖 4 全球之大地起伏值 圖 5 台灣地區之大地起伏值 6
以 GPS 之觀點而言, 一旦大地起伏值之決定可達足夠之精度, 則 GPS 測高 (GPS heighting) 之作法, 將可輕易的透過 GPS 測定之橢球高轉換 為正高的過程而予以達成 以長距離之高程測定而言,GPS 測高之優點 除了 GPS 定位之野外作業較傳統水準測量便利外,GPS 所決定兩點間橢 球高差之精度亦較水準測量所得之正高差為佳 惟一要解決的問題即為 高程系統轉換所需之大地起伏值, 其目前所能提供之精度與水準測高之 誤差相比仍嫌過高, 這是猶待努力之課題 由於 GPS 定位與計算大地起伏值之精度受限之緣故, 相對性的 GPS 測高模式即可加以考量, 期望經由差分模式之運作, 來決定短距離內二 個測站間的精密高差值, 其方法為 或者表為 H = h N (4) H H = ( h h ) ( N N ) (5) 2 1 2 1 2 1 事實上, h 可從 GPS 之差分模式中獲得,N 則需靠重力法大地水準 面模式加以計算獲得 此時由於零階項大地起伏值之主要誤差已不存 在, 而決定 N 之其它誤差亦可有效減小, 因此對於短距離所決定之 H 會 有比較高的精度 4 坐標型態之轉換由前述介紹的幾種坐標相關系統可知, 當地球表面 ( 或空間中 ) 的任一點位欲表達其空間坐標位置時, 在 WGS84 系統中所顯現之坐標會是架構在橢球體上的曲面坐標 ( 如圖 6), 其參數是以經度 緯度 高程 (λ ϕ h) 為主的大地坐標 ( 如圖 7), 而在 ITRF 系統中所顯現之坐標則是架構在一組參考框架上的卡氏坐標, 其參數是以三維空間 (X, Y,Z) 為主的直角坐標 ( 如圖 8), 另針對平面坐標系統而言, 則又有 7
以橫麥卡托 (TM) 投影為主之坐標表達方式 ( 如圖 9), 其參數常以二 維平面 (x,y) 或 (N,E) 來表示 圖 6 球面坐標系統 圖 7 大地 ( 地理 ) 坐標系統 8
圖 8 卡式地心直角坐標系統 圖 9 地圖投影示意圖 ( 中間一組即為 TM 使用之圓柱投影 ) 上述三種坐標型態 ( 大地坐標 地心直角坐標 平面坐標 ) 之間可以利用簡單之數學公式予以轉換, 其中大地坐標化算為地心直角坐標之方式如圖 10, 地心直角坐標化算為大地坐標之直接計算方式如圖 11, 大地坐標化算為 TM 平面坐標之方式如圖 12,TM 平面坐標化算為大地坐標之方式如圖 13 9
圖 10 大地坐標化算為地心直角坐標之計算方式 圖 11 地心直角坐標化算為大地坐標之計算方式 10
圖 12 大地坐標化算為 TM 平面坐標之計算方式 11
圖 13 TM 平面坐標化算為大地坐標之計算方式 12
5 坐標基準之轉換 在太空大地測量 (space geodesy) 技術尚未普及而不易進行全球性 聯測之時期, 世界各國大都採行最密合 (the best fitting) 該地區之 橢球體, 進而建立其地區性的坐標系統 ( 如臺灣地區採用 GRS67 橢球體 之虎子山坐標系統 大陸地區採用 1924 年國際橢球體之南京坐標系 統 ), 而目前則以建立全球性之地心坐標系統為走向, 以期構建能與世 界各國接軌之坐標系統 ( 如臺灣地區最新採用 GRS80 橢球體之 TWD97 坐 標系統 ) 有關大陸與台灣地區在近代所使用之大地基準資料可列如表 2 所示 表 2 台灣及大陸地區近代大地基準 坐標系 橢球體 長半徑 () 扁率 (1/f) 台灣虎子山坐標系 GRS 67 6378160.0 298.247167427 台灣 TWD97 坐標系 GRS 80 6378137.0 298.257222101 GPS 坐標系 WGS 84 6378137.0 298.257223563 大陸南京坐標系 International 1924 6378388.0 297.0 大陸 1954 北京坐標系 Krassovsky 1940 6378245.0 298.3 大陸 1980 西安坐標系 GRS 75 6378140.0 298.257 當一個地區更換一組新的大地基準時, 地圖所呈現之坐標資料勢必伴隨著坐標系統之更動而改變 以我國政府所公布之一組新的國家坐標系統 (TWD97) 為例, 各坐標點成果便與舊有之虎子山坐標系統 (TWD67) 產生一定程度之差異 由於此時大地基準之改變必將造成地圖更新其坐標系統之相關問題, 有效的解決方法之一就是利用基準轉換的方式, 使新舊之坐標系統能夠相互關聯 因此, 當某一地區在某一時間點上必須 13
進行新舊坐標之間的轉換時, 只要其系統轉換模式 ( 轉換參數 ) 能夠確 定, 即可在二組新舊之坐標系統內進行適當之演算程序, 以將地圖之坐 標成果展現在所需要之坐標系統內 ( 如圖 14) 圖 14 基準轉換之執行程序 在決定基準轉換模式及參數之過程中, 通常須同時具備若干點位在兩種坐標系統內之坐標值, 以期經由數量足夠與適當分布的此類地控點 ( 共同點 ), 先行確定坐標系統間的轉換關係, 以供地圖內之其他點位進行坐標轉換 因此, 地控點坐標資料獲取之正確性及其精度如何, 以及地控點分布位置之幾何條件如何, 便成為此項基準轉換之關鍵性因素 由於不同大地基準間進行坐標轉換時, 轉換公式及參數會與所選定區域的大小 地控點分布位置以及數量等因素有關, 進而會影響到坐標轉換後的成果精度 因此, 基準轉換模式之選擇亦扮演相當重要之角色, 一般而言可採行之轉換模式包括 : 14
( 一 ) 七參數轉換在三度空間內同一點的兩組卡氏坐標進行轉換時, 在沒有網形變形的理想情形下, 只需利用七參數即可進行坐標的正交轉換, 其中七參數包括一個尺度參數 三個原點平移參數以及三個軸旋轉參數, 因模式由七個參數組成, 故至少需三個共同之地控點坐標才能加以求解 七參數轉換法較常使用的二種標準模式分別是 : 1. Bursa-Wolf 模式 : (1) 1 1 1 + = z y x S Z Y X Z Y X X Y X Z Y Z 其中包含平移參數 ( X, Y, Z ) 尺度參數 (S) 與旋轉參數 ( X, Y, Z ), 並可假設旋轉參數的值很小 而新 舊兩組坐標係分別以 (X,Y,Z) 與 (x,y,z) 表示 2. Molodensky-Badekas 模式 : 為避免平移與旋轉參數產生高度相關的現象, 可先將坐標原點平移至轉換區域重心後再進行轉換, 即 (2) 1 1 1 + + = X Y X Z Y Z z z y y x x S z y x Z Y X Z Y X 其中 (,, ) 即為轉換區域之重心坐標 x y z 15
( 二 ) 微變公式轉換由於使用七參數轉換至少需要三個地控點同時擁有新 舊系統之坐標值, 因此若在地控點缺少之情形下, 亦可考量使用微變公式法 (differential ethod), 如此只需一個基點之新舊坐標資料以及新舊橢球體之元素資料, 即可求出各點位舊坐標值轉換到新坐標系統上之坐標值, 即 δx δy δz 0 0 0 + δφ J δλ + δh δa B = δf 0 (3) 其中, δ x 0 y0 δ δz0 為基點在新舊坐標系統原點平移量之差值, δ a δ f 為兩基準所使用橢球體之長半徑 扁率之差值, J B 為由轉換點大地坐標及曲率半徑值所組成之轉換式微分矩陣, δφ δλ δ h 為坐標點由舊坐標系轉換至新坐標系之大地坐標改變量 此一模式之轉換雖為美國國防部所採行 ( 亦稱之為標準 Molodensky 轉換公式, 如圖 15), 但由於此種模式基本上只利用一個基點之新舊坐標資料來獲取轉換參數, 且有原理上新舊坐標之間基準差異大能太大之限制, 故其轉換後之坐標精度是否能與七參數模式轉換之成果相對等, 仍是可以探討之課題 16
圖 15 標準 Molodensky 轉換計算方式 在領土廣大的國家 ( 如美國 澳洲 ), 通常政府會對不同地區內提供其各自適用之坐標轉換方法 ( 另如多項式轉換法 曲面擬合法等 ) 及參數, 以進行較精準需求之基準轉換使用 針對台灣地區而言, 內政部亦已正式公告可經由申請程序獲得約 30 c 轉換精度之 (TWD67 與 TWD97 間 ) 七參數基準轉換程式 17