國 立 交 通 大 學
土 木 工 程 學 系
碩 士 論 文
以 GPS 決定船姿態角改善船載海面高
GPS determination of ship attitude for improving
shipborne sea surface height
研 究 生：翟邦和
指導教授：黃金維
I
以 GPS 決定船姿態角改善船載海面高
研究生：翟 邦 和 指導教授：黃 金 維 國立交通大學土木工程學系 摘 要 本研究於船測重力計畫之船上建立一多天線 GPS 系統獲取取樣間隔 1 Hz 之 觀測資料，以 TRACK 進行動態差分定位解算，而後進行船隻運動姿態解算，解 算方法為姿態角直接解法(Direct Computation Formulas, DCF)；而後以此姿態角 成果對天線垂距作改正，得時變天線垂距改正量，再應用時變天線垂距於海面高 (Sea Surface Height, SSH)改正。本研究撰寫了一 FORTRAN 程式，處理資料時間 比對、解算天線相位中心間距離、解算船隻運動姿態、使用姿態成果取得時變天 線垂距並以其改正海面高。TRACK 定位成果在靜態長距離基線中精度可達約 1 公分；動態船隻基線則可達約 2 公分。姿態角成果與同時於船上施測之 INS/GPS 整合系統所得成果進行比較，其成果顯示多天線 GPS 系統與 INS/GPS 整合系統 之航向角差異甚小，故知此多天線 GPS 系統搭配姿態角直接解法於求取船隻運 動姿態上是可行之方法。時變天線垂距改正海面高成果則以交叉點分析、與 DTU10 之差值分析進行比較：與平均海水面 DTU10 之差值，於直線航線上標準 差約改善 0.7 公分，改善比率約為 3.45%，於轉向航線上標準差約改善 1.4 公分 ，改善比率約為 4.07%；交叉點分析方面，方均根值改善約 0.5 公分，比率為 3.31% 。 關鍵字：船載海面高；姿態；多天線 GPS 系統；INS/GPS 整合系統II
GPS determination of ship attitude for improving
shipborne sea surface height
Student：Pang Ho Jai Advisor：Dr. Cheinway Hwang Department of Civil Engineering
National Chiao Tung University
Abstract
In this research, a GPS multi-antenna system has been built in the campaigns of collecting shipborne gravity with 1-Hz data rate. The positioning package TRACK is used for kinematic relative positioning (KRP). The assessment of baseline length suggests that the accuracy of KRP from TRACK reaches 1cm in a static condition and 2 cm on a ship. A computer program has been built to compute ship attitudes from kinematic baselines by the direct computation formulas (DCF). The estimated attitudes are used to correct for the antenna height to the sea surface due to attitude effect. The comparison between GPS and GPS/INS-derived attitudes suggest that GPS attitude system and DCF function properly. With the attitude effect accounted, the RMS crossover differences of the corrected SSHs are improved by 0.5 cm, or 3.31% relatively, compared to the case of raw SSHs. With the ocean tide corrected, the difference between corrected SSHs and DTU10 modeled values are reduced by 0.7 cm in standard deviation along straight survey lines, and by 1.4 cm along non-straight lines, corresponding to relative improvements of 3.45% and 4.07%, respectively.
Keywords: shipborne sea surface height; attitude; kinematic relative positioning;
III
誌謝
6 年，隨著碩士論文的完成，在交大當學生的日子也跟著結束了，在這裡 畢業了兩次，大學和碩班。碩班這兩年，要感謝很多人，首先就是指導教授黃金 維老師，在他適時的指引方向和悉心指導下，讓我能確定目標和前進的方向，他 也常與我討論關於我未來的人生規劃，並給予我許多的建議，讓我在專業領域以 及生活態度上都收穫良多。非常感謝台北大學的陳國華老師，他辛苦的往來新竹 台北指導並給予許多建議，為我節省了許多摸索的時間。還有成功大學的江凱偉 老師與小黑學長，教導我並提供我重要的研究數據，讓這篇論文能夠順利完成。 最要感謝的是爸媽和家人們的體諒、支持和期盼，你們的關心是我努力的動 力，你們的拉拔和栽培造就出現在的我，未來我會努力的不讓你們失望的。 兩年時間，啟訓學長是我想感謝最多的人，他陪我走過低潮，也陪著我一起 研究，如果沒有他，我就沒有現在的成果，我從他那邊得到了許多，真的是由衷 的感謝，希望未來有機會能夠好好的報答學長。感謝逸如學姊，我回到交大的時 候正好是她要離開的時候，但後續她仍給了我許多的幫助和解答我許多的疑惑。 善解人意的綉雯，陪我聊了很多也幫了我很多忙；辛苦耐勞的恩銘學長，在我碩 一的時候幫了我許多的忙；聰明絕頂的詠升，除了課業上的幫忙，也讓大家常常 笑開懷，沒想到當初我講的小鄭博居然一語成讖了；開朗大方又美麗的雅琦，我 對妳的感謝妳懂的；謝謝鐙凱和元旎的許多幫忙；真的很感謝研究室的大家，許 多開心的回憶和大家庭一般的感覺，真的很棒！ 謝謝女友小葉，她是支持我的最大原動力，每天的陪伴和打氣加油是我最後 關頭能一直堅持的大功臣，把我從沮喪和混亂中一次次的扶起，還有她對我的信 心比我自己還多讓我能夠一直督促著自己，辛苦妳了，謝謝妳。聖翔，對你的感 謝就不在這裡多說了。碩班的同學們，世涵、薇帆、柏溶、伶蒨、郁珊、昭儀、 詠升，很懷念那些一起趕作業、專題報告和拼命寫大抄的日子，謝謝你們，祝大 家鵬程萬里。IV
目錄
摘 要 I Abstract II 誌謝 III 目錄 IV 圖目錄 VI 表目錄 IX 第一章 緒論 1 1.1 研究動機與目的 1 1.2 文獻回顧 2 1.3 論文架構 4 第二章 理論基礎 5 2.1 座標系統 52.1.1 地心地固座標框架 (Earth-Centered Earth-Fixed, ECEF) 5
2.1.2 區域水平座標框架(Local level frame coordinate system, LLF) 6
2.1.3 載體座標框架(Vehicle frame coordinate system, VF) 6
2.1.4 天線座標框架(Antenna body frame coordinate system, ABF) 7
2.2 座標框架轉換 8
2.2.1 旋轉矩陣介紹與區域水平座標框架、天線座標框架轉換 8
2.2.2 地心地固座標框架之直角坐標系、大地坐標系轉換 11
2.2.3 地心地固座標框架與區域水平座標框架轉換 12
2.3 姿態角解算: 直接解法(Direct Computation Formulas, DCF) 12
2.4 姿態角改正天線垂距 14
第三章 船動態定位原理與方法 17
3.1 全球衛星定位系統 17
V 3.3 TRACK GPS 動態差分定位解算 21 3.4 TRACK 動態定位精度檢核分析 25 3.4.1 靜態測試 25 3.4.2 外部精度評估 27 3.4.3 TRACK 內部精度評估 33 3.5 改正模式 34 3.5.1 海潮改正模式 34 3.5.2 固體潮改正模式 34 3.5.3 濾波 35 第四章 GPS 與 INS/GPS 資料計算 36 4.1 研究區域簡介 36 4.1.1 施測範圍、儀器、測量船隻簡介 36 4.2 施測與計算方法 41 4.2.1 GPS 41 4.2.2 整合式 INS/GPS 系統 45 4.3 GPS 船載海面高 49 第五章 研究成果及分析 51 5.1 姿態角解算成果 51 5.2 海面高改正成果 63 5.2.1 與平均海水面模式 DTU10 MSS 差異分析 63 5.2.2 交叉點分析 76 第六章 結論與建議 78 6.1 結論 78 6.2 建議 79 第七章 參考文獻 81
VI
圖目錄
圖 2.1 當地水平座標框架(LOCAL LEVEL FRAME COORDINATE SYSTEM)示意圖 ... 6
圖 2.2 載體座標框架(THE VEHICLE FRAME COORDINATE SYSTEM)示意圖 ... 7
圖 2.3 天線座標框架(ANTENNA BODY FRAME COORDINATE SYSTEM)示意圖 ... 8
圖 2.4 船載 GPS 天線垂距示意圖 ... 15 圖 2.5GPS 姿態角改正天線垂距示意圖 ... 16 圖 3.1 GPS 相對定位示意圖 ... 19 圖 3.2 TRACK 動態差分處理流程圖 ... 22 圖 3.32012 年 DOY135 基線 VR01-SHMN 之殘差 ... 26 圖 3.42012 年 DOY135 基線 HL01-SHMN 之殘差 ... 27 圖 3.52012 年 DOY135 基線 HL01-VR01 之殘差 ... 27 圖 3.6TRACK 定位成果計算所得 DC 基線長與全測儀觀測值之差異。 ... 30 圖 3.7TRACK 定位成果計算所得 DA 基線長與全測儀觀測值之差異。 ... 31 圖 3.8TRACK 定位成果計算所得 AC 基線長與全測儀觀測值之差異。 ... 32 圖 3.9TRACK 解算 2012DOY136 之橢球高、天頂距對流層延遲量及軟體內部精 度... 33 圖 4.12012 年近岸船載重力計畫航線規劃圖(摘自內政部 101 年度臺灣本島近岸 船載重力計畫書) ... 37 圖 4.22012 年近岸船載重力計畫實際施測航線圖 ... 38 圖 4.3 尖再發號船體外觀 ... 39 圖 4.4 上為實驗儀器於船隻上架設位置圖、下為全測儀施測數據圖，紅色圓圈及 A、B、C、D 為 GPS 接收天線位置，2 綠色圓圈為 IMU 及其 GNSS 天線位 置。... 40 圖 4.5 岸基站樁標埋設示意圖 ... 42 圖 4.6GPS 資料處理流程圖 ... 43 圖 4.7GPS 動態定位成果應用、改正及分析流程 ... 44 圖 4.8CAINS-21 軟體使用之改良閉合式鬆耦合架構(摘自內政部 100 年度發展與 應用多平台遙測製圖技術工作案總報告書，江凱偉等人，2011) ... 48 圖 4.9 船載海面高測量與衛星測高原理圖(陳逸如,2010) ... 49 圖 5.1DOY135GPS 資料時間匹配後所繪航線圖 ... 52 圖 5.2DOY135 解算所得之船隻航向角：藍線為 GPS 系統解算所得；紅線為 INS/GPS 系統解算所得 ... 52 圖 5.3DOY135GPS 系統解算所得之船隻俯仰角 ... 53 圖 5.4DOY135GPS 系統解算所得之船隻側傾角 ... 53 圖 5.5DOY136GPS 資料時間匹配後所繪航線圖 ... 54 圖 5.6DOY136 解算所得之船隻航向角：藍線為 GPS 系統解算所得；紅線為 INS/GPS 系統解算所得 ... 54
VII 圖 5.7DOY136GPS 系統解算所得之船隻俯仰角 ... 55 圖 5.8DOY136GPS 系統解算所得之船隻側傾角 ... 55 圖 5.9DOY137GPS 資料時間匹配後所繪航線圖 ... 56 圖 5.10DOY137 解算所得之船隻航向角：藍線為 GPS 系統解算所得；紅線為 INS/GPS 系統解算所得 ... 56 圖 5.11DOY137GPS 系統解算所得之船隻側傾角 ... 57 圖 5.12DOY137GPS 系統解算所得之船隻側傾角 ... 57 圖 5.13DOY138GPS 資料時間匹配後所繪航線圖 ... 58 圖 5.14DOY138 解算所得之船隻航向角：藍線為 GPS 系統解算所得；紅線為 INS/GPS 系統解算所得 ... 58 圖 5.15DOY138GPS 系統解算所得之船隻俯仰角 ... 59 圖 5.16DOY138GPS 系統解算所得之船隻側傾角 ... 59 圖 5.17DOY139GPS 資料時間匹配後所繪航線圖 ... 60 圖 5.18DOY139 解算所得之船隻航向角：藍線為 GPS 系統解算所得；紅線為 INS/GPS 系統解算所得 ... 60 圖 5.19DOY139GPS 系統解算所得之船隻俯仰角 ... 61 圖 5.20DOY139GPS 系統解算所得之船隻側傾角 ... 61 圖 5.21DOY135GPS 系統與 IMU 系統解算所得之航向角差值 ... 62 圖 5.22DOY136GPS 系統與 IMU 系統解算所得之航向角差值 ... 62 圖 5.23DOY137GPS 系統與 IMU 系統解算所得之航向角差值 ... 62 圖 5.24DOY138GPS 系統與 IMU 系統解算所得之航向角差值 ... 63 圖 5.25DOY139GPS 系統與 IMU 系統解算所得之航向角差值 ... 63 圖 5.262012 年 DOY135~139 船測編輯後所得直線航線圖 ... 66 圖 5.27 所選航線 ROUTE 01 經姿態角改正天線垂距後所得之海水面高成果 ... 67 圖 5.28 所選航線 ROUTE 04 經姿態角改正天線垂距後所得之海水面高成果 ... 67 圖 5.29 所選航線 ROUTE 06 經姿態角改正天線垂距後所得之海水面高成果 ... 67 圖 5.30 所選航線 ROUTE 08 經姿態角改正天線垂距後所得之海水面高成果 ... 68 圖 5.31 所選航線 ROUTE 09 經姿態角改正天線垂距後所得之海水面高成果 ... 68 圖 5.32 所選航線 ROUTE 10 經姿態角改正天線垂距後所得之海水面高成果 ... 68 圖 5.33 所選航線 ROUTE 11 經姿態角改正天線垂距後所得之海水面高成果 ... 69 圖 5.34 所選航線 ROUTE 12 經姿態角改正天線垂距後所得之海水面高成果 ... 69 圖 5.35 所選航線 ROUTE 13 經姿態角改正天線垂距後所得之海水面高成果 ... 69 圖 5.36 所選航線 ROUTE 14 經姿態角改正天線垂距後所得之海水面高成果 ... 70 圖 5.37 所選航線 ROUTE 15 經姿態角改正天線垂距後所得之海水面高成果 ... 70 圖 5.38 所選航線 ROUTE 16 經姿態角改正天線垂距後所得之海水面高成果 ... 70 圖 5.39 所選航線 ROUTE 17 經姿態角改正天線垂距後所得之海水面高成果 ... 71 圖 5.40 所選航線 ROUTE 18 經姿態角改正天線垂距後所得之海水面高成果 ... 71 圖 5.41 所選航線 ROUTE 20 經姿態角改正天線垂距後所得之海水面高成果 ... 71
VIII 圖 5.42 所選航線 ROUTE 22 經姿態角改正天線垂距後所得之海水面高成果 ... 72 圖 5.43 所選航線 ROUTE 24 經姿態角改正天線垂距後所得之海水面高成果 ... 72 圖 5.44 所選航線 ROUTE 26 經姿態角改正天線垂距後所得之海水面高成果 ... 72 圖 5.45 所選航線 ROUTE 27 經姿態角改正天線垂距後所得之海水面高成果 ... 73 圖 5.46 所選航線 ROUTE 28 經姿態角改正天線垂距後所得之海水面高成果 ... 73 圖 5.47 所選航線 ROUTE 29 經姿態角改正天線垂距後所得之海水面高成果 ... 73 圖 5.482012 年 DOY136 船測編輯後所得船轉向航線圖 ... 75 圖 5.49 交叉點差異量分析示意圖 ... 76
IX
表目錄
表 1.1 近期船載測高研究之海水面高精度分析 ... 3 表 3.1HL01、SHMN 、VR01、 VR02 之經度(LON)、緯度(LAT)及高(H)之殘差值 統計... 25 表 3.2 測站間組基線之基線長殘差值統計 ... 26 表 3.3 全測儀基線長、TRACK 平均基線長、兩者差值 ... 29 表 3.4TRACK 基線長與全測儀所得基線長差異值統計 ... 29 表 4.1 船載重力計畫選用船隻─尖再發號 ... 39 表 4.2 測線起迄座標(摘自內政部 101 年度臺灣本島近岸船載重力計畫書) ... 41 表 4.3GPS、INS 及 INS/GPS 整合系統之特色比較 ... 45 表 4.4 一般 INS/GNSS 整合架構的簡易比較(CHIANG,2004) ... 47表 5.12012 年 DOY135~139 航向角(YAW)GPS 系統成果與 IMU 成果比較 ... 62
表 5.22012 年 DOY135~139 船測姿態角改正前後分析 ... 66
表 5.32012 年 DOY136 船測轉彎處姿態角改正成果分析 ... 75
1
第一章 緒論
1.1 研究動機與目的
海面高(sea surface height, SSH)為海洋及地球重力場應用中重要的一環。地球
重力在許多研究領域中佔有著重要的位置，它能夠反映地球形狀和地殼結構的重 要資訊；在大地測量領域，地球重力場除了能研究地球樣貌外，還可以確定高程 座標；在地球物理領域，地球重力場的觀測資料可以用來探勘地層中可能的組成 以發掘礦藏資源，還能夠對地震成因進行分析，此外也有學者嘗試利用觀測重力 場的變化情形進行地震預測；在國防應用領域上也有相當重要的作用(施宣昶, 2004)。除了以相對或絕對重力儀等科技搜集陸地上的重力資訊外，海洋重力測 量的資料亦可提供做為國家高精度基本測量、海洋資源開發、海面與海底工程設 計和施工、海底地殼活動和海洋潮汐之研究等運用，故海洋重力測量具有相當之 重要性。 船載測高與衛星測高的目的一致，主要為測定海洋大地水準面與海洋地形等 資訊。衛星測高的原理為應用衛星測高儀取得瞬間海水面至衛星之間的高度，在 選定一固定的參考橢球體之後，經過海潮改正，則可以獲得海水面與參考橢球體 之間的海面高度，然而衛星測高在近岸地區，雷達波的雜訊相對放大許多，導致 近岸區域的測高品質不佳，而船載測高則無此項問題，因此在近岸海面高資訊的 獲取還是需要依靠船載測高技術。 本研究應用內政部「臺灣本島近岸船載重力測量」計畫中的船載定位GPS 觀測資料，船載GPS獲得的橢球高度，經過天線垂距、海潮、固體潮等改正後， 可以獲得與衛星測高相當性質的海面高。動態差分定位(Kinematic Relative Positioning, KRP)與精密單點定位(Precise Point Positioning, PPP)目前在海測上所 得的定位精度表現相當接近(陳逸如, 2010)，故本研究選擇以動態差分方法進行 GPS定位解算。
2 時，便以GPS定位所得之橢球高減去此定值天線垂距；然而船隻在航行過程中必 定會因船速及海況等因素而產生姿態的變化，此一變化使得船身傾斜，進而改變 了GPS訊號接收天線相位中心到海水面的垂直距離，此垂直距離即為天線垂距改 正所需改正之數值，倘若能得到船隻運動姿態的變化並將此姿態成果用於改正天 線垂距，應能改善海水面高之成果。本研究為解算船隻於施測時的運動姿態，除 原船載重力計畫所裝設的GPS接收儀外，另外增設3組GPS接收儀，並使用姿態 解直接解法進行姿態角的求解，然後將此姿態角成果應用於改正天線垂距，使其 不再是一定值，而是一時變數值，期望能夠改善海面高之成果。
1.2 文獻回顧
船載測高，依Chadwell and Bock (2001)所提出之實驗結果，短基線GPS浮標 測高經平滑化後精度為2.4公分；Rocken et al. (2005)提出的實驗成果，使用長基 線船載儀器精度為10公分；Foster et al. (2009)所提出之實驗結果，使用TRACK 解算1 Hz GPS資料之基線解海水面高結果，與tidal height比較差異值RMS為9.3 公分，使用濾波五分鐘達8.3公分；與model SSH做比較，差異值RMS為56公分， 使用五分鐘濾波罩窗差異值RMS達16.1公分；與model SSH做比較並引入驗潮站 潮位資料改正後差異值RMS為69公分，使用五分鐘濾波罩窗差異值RMS達13.3 公分。Bouin et al. (2009)提出的研究成果，使用GAMIT解算1Hz船載GPS資料之 交叉點海面高差異量為10~14公分，以Sea Level Anomalies, SLA改正後，相對於 DNSC08所得RMS為27.3公分，相對於EGM08為21.1公分；陳逸如(2010)提出的 實驗成果，使用Bernese PPP解算1 Hz船載GPS資料經平差後之交叉點海面高差異 量，精度為2~16公分。相關文獻整理如表1.1所示。
而在GPS 多天線系統解算姿態角方面， Lu(1995)提出以直 接解法 (Direct Computation Formulas, DCF)或最小二乘配至法(Implicit Least Squares, LSQ)進行 姿態角解算，兩種方法各有其利弊，DCF之優點在於其解算快速，但只能針對已 選取之3組GPS資料進行姿態角求解； LSQ之優點在於能對任意3組GPS資料進行
3 求解，多於之GPS觀測資料則可用作多餘觀測量，但其解算較為繁複且耗時。 表 1.1 近期船載測高研究之海水面高精度分析 案例 精度(公分) 經濾波或平 滑化處理精 度 (公分) 參考資料 短基線 GPS 浮筒(buoy)研 究
- 2.4 Chadwell and Bock,
2001 長基線船載研究設備 10 - Rocken et al., 2005 1 Hz 船載 GPS 以動態基線 解算海水面高，相對於所量 測到的海潮訊號為參考面 9.3 8.3 Foster et al., 2009 1 Hz 船載 GPS 以動態基線 解算海水面高，相對於一 SSH 模型為參考面 56 16.1 Foster et al., 2009 1 Hz 船載 GPS 以動態基線 解算海水面高，相對於一 SSH 模型為參考面以及使 用驗潮資料改正 69 13.3 Foster et al., 2009 1Hz 船載 GPS 以動態基線 解算交叉點海面高差異量 10~14 - Bouin et al., 2009 1Hz 船載 GPS 以動態基線 解算並以 SLA 改正，相對 於 DNSC08 27.3 - Bouin et al., 2009 1Hz 船載 GPS 以動態基線 解算並以 SLA 改正，相對 於 EGM08 21.1 - Bouin et al., 2009 1Hz 船載 GPS 以精密單點 定位解算交叉點海面高差 異量 - 2~16 陳逸如, 2010
4
1.3 論文架構
第一章：緒論 包括研究動機與目的、文獻回顧與本文章節架構之說明 第二章：座標框架轉換與姿態角解算方法 本章共分4個部分：說明本研究中所使用之各座標框架、各座標框架間相互 轉換之方法、講述姿態角直接解法之理論、姿態角改正天線垂距之方法。 第三章：船動態定位原理與方法 本章主要重點在於說明GPS動態定位之原理和對定位成果進行精度分析：對 全球衛星定位系統作一簡單介紹；講述動態差分定位之原理；介紹TRACK，本 研究所使用之GPS動態差分解算軟體；TRACK定位成果精度分析；改正模式之 介紹。 第四章：GPS與INS/GPS資料計算 本章說明2012年船載計畫之研究區域、使用之儀器；並簡介GPS系統、 INS/GPS整合系統和各自作業流程；最後講述GPS船載海面高之基本原理。 第五章：研究成果及分析 此章共分兩大主題：姿態角解算成果與海面高改正成果。姿態角改正成果方 面，先將GPS系統所得航向角與INS/GPS所得航向角進行比較，而後對於GPS系 統之姿態角成果進行討論；海面高改正成果方面，以與側高衛星所得之平均海水 面參考場DTU10比較、交叉點海面高差異值分析，共兩項進行討論。 第六章：結論與建議5
第二章 理論基礎
2.1 座標系統
在定位、導航中，座標系統非常重要，因此將在此對本研究中所使用之座標 系統進行介紹。
2.1.1 地心地固座標框架 (Earth-Centered Earth-Fixed, ECEF)
地心地固座標框架與協議地球座標系統(Conventional Terrestrial System, CTS) 之三軸相符合，亦稱為平均地球座標系統(Average Terrestrial System, AT)，為全 球性之地心座標系統。依表達方式可分為直角座標系及大地座標系。直角座標系
之原點為地球質量中心，系統 Z 軸平行於 BIH(Bureau International de l'Heure)定
義之 CTP(Conventional International Pole)，X 軸為地心指向 BIH 定義之零子午圈 (格林威治子午圈)與赤道面交點之方向，Ｙ軸則是與上述Ｘ、Ｚ軸正交面構成一 右旋正交系統(Seeber, 2003)，座標一般以(X, Y, Z)表示。目前美國的 GPS(Global Positioning System)所建立的 WGS84(World Geod,etic System 1984)座標系統即是 協議地心地固座標系統的具體實現，是美國國防部的 GPS 廣播星曆和精密星曆 的參考座標系統(Misra and Enge, 2011)。GPS 系統和慣性導航中所提到的地固座 標框架通常是使用 WGS84 座標系統。(胡智祐, 2009)
大地座標系(Geodetic coordinate system)之定義為地球橢球的中心與地球質
心重合，橢球的短軸與地球自轉軸相合，一般以( ,φ λg g, h )g 分別表示大地緯度、
經度及橢球高。大地緯度ϕg定義為過地面點之橢球法線與橢球赤道面的夾角，
大地經度λg則是過地面點的橢球子午面與格林威治平大地子午面之間的夾角，
6
2.1.2 區域水平座標框架(Local level frame coordinate system, LLF)
當地水平座標框架的原點為測站的中心，在本文中為多天線系統中的主天 線，Z 軸的方向定義為測站的橢球面法線方向，Y 軸指向地平北方(橢球子午圈 方向)，X 軸則與上述 Z 軸、Y 軸組成一右旋正交座標框架，此三軸分別指向北、 東、地，所以稱為 NED(North-East-Down)座標系統。在此篇論文中所使用的則 是另一三軸定義為東、北、天(East-North-Up, ENU)的座標系統(圖 2.1)，此座標 系統在慣性導航系統中是被廣泛使用的當地水平座標框架，與 NED 座標系統不 同的是其 Z 軸指向為橢球面法線的反方向，也就是當高程增加時 Z 軸是正向的。
圖 2.1 當地水平座標框架(Local level frame coordinate system)示意圖
2.1.3 載體座標框架(Vehicle frame coordinate system, VF)
載體座標框架為一正交的座標框架，本文中所使用之載體為船，其 X 軸、
Y 軸分別指向載體的右側與前方，Z 軸則與 X 軸、Y 軸構成一右旋座標系統，此 三軸即為載體的三正交軸，分別對應的旋轉角分別稱為 pitch、roll 及 yaw(heading)， 定義逆時鐘方向為正，如圖 2.2 所示。
7 Z Y (forward) X pitch roll yaw (heading) platform
圖 2.2 載體座標框架(The vehicle frame coordinate system)示意圖
2.1.4 天線座標框架(Antenna body frame coordinate system, ABF)
基於空間中三個點可以定義一平面，天線座標框架就是選定三個 GPS 天線 來組成一個面，當天線平面組成後，就可以定義出天線座標系統。假設載具上架 設了三個 GPS 接收儀(天線)，分別為 Antenna 1(Ant1)、Antenna 2(Ant2)和 Antenna 3(Ant3)，將原點設定在 Ant1，意即其座標為(0, 0, 0)；Y 軸為沿 Ant1 至 Ant2 組
成的基線(baseline)方向，因此 Ant2 之座標為(0, L12, 0)，當中 L12為 Ant1 至 Ant2
之長度；X 軸則是在由 Ant1、Ant2 和 Ant3 組成的平面上，正交於 Y 軸並指向
右方，Ant3 之座標是由向量 Ant1Ant2、Ant1Ant3 及餘弦定理求出α 角後定義
出而得(L13sinα, L13cosα, 0)；最後 Z 軸與上述 X 軸、Y 軸組成一右旋座標系統。
這個座標系統也稱作自體座標框架(body frame coordinate system)。圖 2.3 為天線 座標框架示意圖。
8
圖 2.3 天線座標框架(Antenna body frame coordinate system)示意圖
2.2 座標框架轉換
在慣性導航中，座標框架的轉換非常重要。本研究中所使用之 GPS 與 IMU 皆須經由座標框架的轉換，而後進行各種解算與比較。如本研究以 TRACK 所解 算出之地心地固座標需轉換至當地水平座標框架下並利用天線座標系統才能進 行姿態角的解算，而後修正海水面高時則須要使用到載體座標系統之概念。IMU 所量測到的值為慣性導航中最原始的觀測量 b f 、ω ，其所使用的座標框架為感ibb 測器座標系統，但在導航中，習慣上將座標轉換到當地水平座標框架或是地心地 固座標框架上。2.2.1 旋轉矩陣介紹與區域水平座標框架、天線座標框架轉換
本研究中所使用之座標系統皆為卡氏座標系統(Cartesian coordinates system)， 要將兩個卡氏座標系統結合需經過兩種轉換，一是平移(translation)，另一是旋轉 (rotation)，在姿態角的解算中，平移對座標系統的三軸方向不會有任何影響，故 非此研究所關注之重點；旋轉則是座標軸系統轉換和求解姿態角的關鍵。以下將 針對座標軸轉換的旋轉矩陣(rotation matrix)進行介紹。
9
假設兩卡氏座標系統原點相重合，但各自之三軸皆指向不同方向，兩者可透 過 3 個旋轉矩陣分別針對三軸進行座標系統旋轉以達成目標 (Wertz et al., 1978; Schwarz and Krynski, 1992)
( )
1 1 1 1 1 1 1 0 0 0 cos sin 0 sin cos R θ θ θ θ θ = − ，對 x 軸旋轉 (2.1a)( )
2 2 2 2 2 2 cos 0 sin 0 1 0 sin 0 cos R θ θ θ θ θ − = ，對 y 軸旋轉 (2.1b)( )
3 3 3 3 3 3 cos sin 0 sin cos 0 0 0 1 R θ θ θ θ θ = − ，對 z 軸旋轉 (2.1c) (2.1a)、(2.1b)及(2.1c)分別為針對右旋系統之 x、y 及 z 軸的旋轉矩陣依照旋轉順 序的不同可以有各種旋轉矩陣。舉例來說，先對 x 軸做旋轉，接著對y′軸旋轉， 最後對z′′軸旋轉，這組旋轉矩陣以 1-2-3 (R123) 為名，即 R1(θ1)R2(θ2)R3(θ3)。而 在本研究中所使用之旋轉矩陣之排序為 3-1-2(R312)，這是在航海上及慣性導航中 被廣泛使用的(Loncarevic, 1993; Wong, 1988)： 312( , , ) 2( ) 1( ) 3( ) R ψ φ λ =R λ R φ R ψ =cos( ) cos( ) sin( ) sin( ) sin( ) cos( ) sin( ) sin( ) sin( ) sin( ) sin( ) cos( ) cos( ) sin( ) cos( ) cos( ) sin( ) sin( ) cos( ) cos( ) sin( ) sin( ) sin( ) sin( ) cos( ) sin( ) cos( ) cos( ) cos( )
λ ψ λ φ ψ λ ψ λ φ ψ λ φ φ ψ φ ψ φ λ ψ λ φ ψ λ ψ λ φ ψ λ φ − + − − + − (2.2)
其中的旋轉角φ 、λ、ψ 分別為尤拉角(Euler angles)中的俯仰角φ (pitch)、側傾角
10 而各旋轉角與旋轉矩陣中元素 R(i, j)之間的關係為： 1 (2,1) tan ( ) (2, 2) R R ψ = − − (2.3a) 1 sin ( (2, 3))R φ = − (2.3b) 1 (1, 3) tan ( ) (3, 3) R R λ= − − (2.3c) 旋轉矩陣(2.2)可將參考座標中之向量轉換至目標座標系當中，譬如以區域水平座 標系(LLF)之 l l l T x y z 乘上旋轉矩陣後即可重合至自體座標框架(BF)上 T b b b x y z ，即： 2( ) 1( ) 3( ) b l b l b l x x y R R R y z z λ φ ψ = (2.4) 旋轉矩陣有各種組合方式，但它們有共通的特性：正交性(Orthogonality) ( , , ) T( , , ) Rψ φ λ R ψ φ λ = ; I R−1( , , )ψ φ λ =RT( , , )ψ φ λ (2.5) 由其正交性可以清楚的知道旋轉矩陣中的 3 個元素φ 、λ及ψ 皆獨立，且其逆矩 陣與轉置矩陣相等。若R( , , )ψ φ λ 為將參考座標系轉至體座標系之旋轉矩陣， ( , , ) T R ψ φ λ 即為將體座標系轉至參考座標系之旋轉矩陣。 此外，旋轉矩陣之乘積亦為正交矩陣，如(2.6)所示： I I G H G H R =R R (2.6)
11 當中 I G R 為將 G 座標框架轉至 I 座標框架之旋轉矩陣， G H R 為將 H 座標框架轉至 G 座標框架之旋轉矩陣，而 I H R 即是將 H 座標框架轉至 I 座標框架。
2.2.2 地心地固座標框架之直角坐標系、大地坐標系轉換
地面任一點之地固座標皆可以直角坐標系與大地坐標系展現，一般直角坐標 系以(X, Y, Z)表示；大地座標系以( ,φ λg g, h )g 表示，兩者轉換關系如下： g g g g g 2 g g ( ) cos cos ( h ) cos sin [ (1 ) ]sin X N h Y N Z N e φ λ φ λ φ φ = + = + = − + (2.7) 當中 N 為參考橢球之卯酉圈曲率半徑，e 為參考橢球的第一偏心率： 2 2 2 2 a b e a − = (2.8) 1 2 2 2 g (1 sin ) a N e φ = − (2.9) a、b 分別為參考橢球之長、短半徑。 直角座標轉換為大地座標： 1 2 1 g 2 2 1 g 2 2 g tan (1 ) tan h cos Z N e N h X Y Y X X Y N φ λ ϕ − − − = − + + = + = − (2.10) 其中φ φ 與 g hg需漸進疊代計算。12
2.2.3 地心地固座標框架與區域水平座標框架轉換
從地心地固座標框架轉換到區域水平座標框架如(2.11)：
sin cos sin sin cos
sin cos 0
cos cos cos sin sin
LLF ECEF R ϕ λ ϕ λ ϕ λ λ ϕ λ ϕ λ ϕ − − = − − − − (2.11) 由旋轉矩陣正交性可得區域水平座標框架至地心地固座標框架：
sin cos sin cos cos
( ) sin sin cos cos sin
cos 0 sin ECEF LLF T LLF ECEF R R ϕ λ λ ϕ λ ϕ λ λ ϕ λ ϕ ϕ − − − = = − − − (2.12)
2.3 姿態角解算: 直接解法
(Direct Computation Formulas, DCF)依照 2.1.2.1 節所示，解算姿態角參數於數學理論上的必要條件為每條基線 都要有一組參考座標系和組成一天線座標系。GPS 多天線系統解算載具姿態，所 使用的參考座標系統為區域水平座標系統，由 GPS 動態差分定位完成後轉換所 得，而天線座標系統則是如 2.1.1.5 節定義一主天線為原點，再由另兩個 GPS 接 收天線共同組成天線座標框架。在公式的推導上需要區域座標系及天線座標系， 然而在實際解算上將只需要區域座標系之定位座標，此節將詳細解說。 此節中以Aib = xib yib zib 代表天線座標框架之座標、T Ail = xil yil zilT 代表區域水平座標框架之座標。天線座標框架的定義如 2.1.1.節及圖 2.X 所示， 乃以 Antenna 1 為原點、Antenna 1 至 Antenna 2 基線方向為 y 軸；區域水平座標 框架之定義則如 2.1.1.3 節中所述，東、北、天(East-North-Up, ENU)的右旋座標 系統。而根據以上座標框架定義可得：
13
[
]
2 2 2 2 0 12 0 T T b b b b A =x y z = L[
]
3 3 3 3 13sin( ) 13cos( ) 0 T T b b b b A =x y z = L α L α[
]
2 2 2 2 2 2 2 T T l l l l A =x y z = E N U[
]
3 3 3 3 3 3 3 T T l l l l A =x y z = E N U 此兩座標系之間的轉換則如(2.13)所示： 2( ) 1( ) 3( ) b l i i A =R λ R φ R ψ A (2.13) 將 2b A 及 2l A 帶入(2.13)後可得到下式 2 12 2 12 2 12 cos( ) sin( ) cos( ) cos( ) sin( ) E L N L U L φ ψ φ ψ φ − = (2.14) 由(2.14)中可推得航向角ψ 及俯仰角φ 之公式： 1 2 2 tan (E /N ) ψ = − − (2.15) 1 2 2 2 2 2 tan (U / E N ) φ = − − + (2.16) 當航向角ψ 及俯仰角φ 求得後，接著以 3 3( ) l A R ψ 將A3l對 z 軸旋轉ψ 得A ′3l ，而 後以A R3l′ 1( )φ 對 x′ 軸旋轉φ 得 3 3 3 3 T l A′′=E′′ N′′′ U ′′ ，最後將 3 l A ′′乘以R2( )λ 便可14 得到A3b =
[
L13sin( )α L13cos( )α 0]
T，而由於R1( )φ 、R2( )λ 及R3( )ψ 同樣皆為正 交矩陣，故 1 2 ( ) 2 ( ) T R − λ =R λ ，由此得： 3 13 3 12 13 3 cos( ) sin( ) cos( ) sin( ) sin( ) E L N L L U λ α α λ α ′′ ′′ = − ′′ (2.17) 1 3 3 tan (U /E ) λ= − − ′′ ′′ (2.18) 方程式(2.15)、(2.16)及(2.18)即為求解姿態角之直接解法(DCF)。2.4 姿態角改正天線垂距
船載 GPS 對海水面高之測量需要經過各種改正，包括定位精度的確定、海 潮模式改正、固體潮改正及天線垂距改正等。本研究著重於增進天線垂距改正， 以求改善海面高之測量成果。 現今一般處理天線垂距改正的方法是於船出航前對 GPS 天線相位中心至海 水面之距離進行量測，而後以此量測所得之值為一常數定值作為天線垂距改正依 據，這在理想狀況上是可行的；然而實際上，儘管船在航行施測期間已維持慢且 穩定之船速，仍會因風浪造成船隻有幅度大小不一的擺動。三維空間中的擺動共 分三方向：改變航行方向所造成的 Z 軸航向角ψ 旋轉變化、因風浪及航行造成 船身於 X 軸俯仰角φ 、Y 軸側傾角λ的旋轉，如圖 2.4 所示。當中航向角ψ 的變 動並不會影響天線垂距，然而俯仰角φ 和側傾角λ則會使得 GPS 天線相對 Z-X 平面及 Z-Y 平面產生變動，進而改變真實天線垂距。15 圖 2.4 船載 GPS 天線垂距示意圖 一般求海面高(SSH)的天線改正公式： con SSH = −h H (2.19) 當中 h 為大地座標系之橢球高，即 TRACK 動態差分解算所得； con H 表示船出航 前所量得之天線垂距定值改正常數。而實際上船隻在運動狀態時，其天線垂距不 可能為一定值，因此本研究藉由多天線 GPS 系統求得姿態角參數後，便可對天 線垂距進行改正以取得各時刻之真實天線垂距，其改正示意圖如圖 2.5， con H 表 示船出航前所量得之天線垂距定值常數， att i H 則是經過姿態角修正後之天線垂距， 紅色天線為該時刻 GPS 天線真實狀況，與 Z 軸有一夾角θ，而 θ 為天線經過姿 態角φ 、λ的旋轉後與 Z 軸之夾角，黑色天線則表示船出航前之狀況。以下將論 述如何透過船隻運動姿態角對天線垂距進行改正進而求得正確海面高(Jekeli,
16 2001)： corr att i i SSH = −h H (2.20)
[
0 0 1][
312]
cos( ) cos( ) att con i H = R �E =H λ φ (2.21) corr i SSH ：第 i 秒之改正後(corrected)海水面高，是利用橢球高 h 減掉第 i 秒經姿 態 角 改正 之天 線垂 距 att i H 後所 得 ； R312 為 公 式(2.2)之 3*3 旋轉矩陣；E： 0 0 Hcon T 。將(2.21)帶入(2.20)後得： cos( ) cos( ) corr con i SSH = −h H λ φ (2.22) 本研究在求得姿態角後即是利用(2.22)進行每一時刻之真實天線垂距改正， 改正所得海面高成果將於第五章進行分析討論。θ
X
Y
GPS Antenna Sea Surface (XY plane) Perpendicular Distance to Sea Surface After Attitude Correction Before Attitude Correction Center of Mass att iH
conH
圖 2.5 GPS 姿態角改正天線垂距示意圖17
第三章 船動態定位原理與方法
3.1 全球衛星定位系統
全球衛星定位系統 GPS(Global Positioning System)為美國國防部為軍事目地 所發展之全球定位系統，而後開放給學術和商業等方面進行研究與應用，目前由 30 顆衛星分佈於 6 軌道面組成系統，其中 24 顆為操作衛星，6 顆為主動備用(active spare)衛星，意即除主衛星失效和加強星群幾何分佈的任務外，平時也可以由此 6 顆衛星的訊號進行定位。GPS 的定位原理是由地面接收儀接收 GPS 衛星所發 射的無線電訊號，進行單程測距，解算未知點(接收儀)之三維空間座標。由於 GPS 定位除三維座標外，還有衛星和接收儀間的時錶誤差共 4 個未知數，因此至 少需同時觀測 4 顆衛星才能用空間後方交會法進行求解，而 GPS 衛星的軌道設 計正可使全球各地在任何時刻同時收到 4 顆衛星的訊號。 GPS 衛星到接收儀天線間的距離可由兩種訊號進行解算：電碼和載波相位。 電碼觀測是利用衛星發射電碼訊號到接收儀間的時間差，得到衛星和接收儀天線 之間的距離(Leick, 1994)。時間是利用接收儀自行產生與衛星發射訊號相同的虛 擬訊號(Pseudo-Random Noise, PRN)，與接收衛星的訊號進行相關分配，當兩訊 號達到最大相關十，將得到傳播訊號的時間延遲(time delay)，距離可由光速與所 得知時間差的乘積得到。另一方面，距離也可利用載波訊號從衛星到接收儀間所 傳遞的總波數，包含整數部分與非整數部分，與波長的乘積獲得(Langley, 1993)。 非整數部分的波數可由接收儀接收的衛星訊號及接收儀內部振盪器的相位差得 知。由衛星到接收器波數之整數部份我們稱之為整數週波值，為一未知且必須求 得的數。若接收儀同時有載波及電碼訊號，通常兩者都會使用。 電碼與載波的頻率皆由基本頻率 10.23MHz 而來。P 電碼隻頻率為基本頻率； C/A 電碼(Coarse/Acquisition)則為基本頻率的十分之一(1.023MHz)。載波訊號 L1 頻率為基本頻率的 154 倍；L2 則為基本頻率的 120 倍。因此，C/A 電碼波長約 為 293 公尺，P 電碼波長則為 29.3 公尺。載波訊號 L1 及 L2 之波長分別為 19 公
18 分及 24 公分。由於載波訊號波長比較電碼訊號為短，因此載波訊號在較要求高 精密的應用上是比較被選擇使用的(El-Rabbany, 2002)。 GPS 定位方法有二大類：相對定位(relative positioning)，為本研究所使用之 定位方法，將在 3.2 節中詳細介紹；單點定位(point positioning)，以單一觀測站 所接收之衛星訊號進行定位，通常用於導航定位和即時定位。精密單點定位 (Precise Point Positioning, PPP)是近來的一項趨勢，此技術廣泛應用於 GPS 動態 定位中(Kouba and Herox, 2001; Satirapod and Homniam, 2006; Ge etal., 2008)，使 用了國際 GPS 服務機構(Global Navigation Satellite System, GNSS)所提供之 GPS 精密星曆和精密時錶誤差，以無電離層影響的載波相位和電碼組合觀測值為觀測 資料，對測站的位置、接收儀時錶誤差、對流層延遲以及組合後的整數週波未定 值等參數進行估算。使用者僅需使用一含雙頻雙碼 GPS 接收儀就可以作高精度 的定位，唯 IGS 所提供之精密星曆需在數據接收後兩週才開放下載。PPP 必須採 用精密的模型加以改正和使用輔助參數進行估計，例如衛星天線相位中心偏差修 正、固體潮修正等。 GPS 之定位誤差來源共三大類，分別是衛星端誤差、訊號傳播誤差及接收設 備相關誤差(包含人為)。衛星端誤差包含衛星軌道誤差及衛星時錶誤差，前者可 透過使用高精度精密星曆消除，後者則透過差分運算消除。訊號傳播誤差包含電 離層延遲誤差、對流層延遲誤差及多路徑效應，電離層延遲可利用雙頻觀測或軟 體中的電離層模式進行消除；對流層延遲誤差又包含了乾、濕分量延遲，乾分量 延遲約佔其中的 90%，一般的處理方法包括避免採用低角度之衛星觀測量，因此 一般會設定在 15 度以上，另外會採用對流層數學模式加以修正；多路徑效應來 自於接收天線週邊環境的反射訊號，透過架設天線時對周遭環境的注意、增加觀 測時間或多餘觀測等即可將此誤差降低。接收設備誤差包含有觀測誤差、接收儀 時錶誤差、週波未定值與週波脫落、天線相位中心位置偏差、天線高量測等，觀 測誤差屬偶然誤差，透過適當增加觀測量即可減弱其影響；接收儀時錶誤差可由
19 差分處理消除；週波未定值與週波脫落經由 GPS 定位計算軟體即可偵測與解決； 天線相位中心則須在施測前先行統一率定；天線高的測量是依靠人為觀測，只需 小心觀測即可。
3.2 相對定位：動態差分定位
相對定位的基本原理在於使用兩部衛星訊號接收儀，其中一部需設置在已 知座標點上，也就是衛星控制點，稱為主站或基站(base station)，而另一部則設 置在待測座標點上，是為待測站或稱移動站(rover)。測得兩測站間的基線向量後， 再根據主站的已知座標進行計算得到移動站座標，如圖 3.1 所示。若待測站在接 收衛星訊號的過程中相對一固定點無任何位移量，即為靜態觀測；反之，若對一 固定點處於運動狀態，即為動態觀測，此類差分定位稱為動態差分定位(Kinematic Relative Positioning, KRP)。 圖 3.6 GPS 相對定位示意圖 差分定位系統採用相對定位的原理，即同時對 2 個或 2 個以上的測點進行訊 號接收和聯合解算，利用觀測方程式之間的相減消除多餘參數，如時錶誤差和電 離層延遲誤差，再以此相減後之方程式進行必要參數的計算，如座標值，差方的 方式可以藉此消除大部分測站間的共同誤差而獲得較高精度之定位成果。差分方 式包含一次差、二次差、三次差，當中一次差可消除衛星時錶誤差及軌道誤差； 二次差可更進一步消除接收儀時錶誤差；三次差則可消除週波未定值。以上差分
20 模式皆有助於減低固體潮影響、電離層以及對流層的延遲誤差(陳逸如, 2010)。 而根據誤差改正量是即時或非即時傳給待測站，可分為即時動態差分和後處理動 態差分，本研究所使用之方法即為後處理動態差分。 後處理動態差分將移動測站所收集之 GPS 數據寫入儲存設備，待測量任務 完成後取出並使用軟體將基站與移動站之資料進行處理解算以達到差分目的。為 了得到高精度的定位成果，可使用高精度的精密星曆進行解算，另外也可以對觀 測資料的品質進行篩選，包括檢查載波相位資料品質、電碼資料品質、觀測資料 平滑化、偵測週波脫落與週波脫落補償等(陳逸如, 2010)。本研究所使用之高精
度 精 密 星 曆 可 由 IGS(International GNSS Service, http://igscb.jpl.nasa.gov/) 或
CODE(Center for Orbit Determination in Europe, ftp://ftp.unibe.ch/aiub/CODE/)自行
下載，IGS 及 CODE 免費提供各級精密星曆供各界使用，唯高精度精密星曆須於 觀測當日起 12~18 日後方得下載。對於船載重力測量、地層下陷研究等非即時性 但需要高精度的測量工作，後處理動態差分是較理想的選擇；其缺點為無及時通 訊數據傳輸，無法即時獲得定位成果以檢視，意即若有操作疏失或數據資料品質 不佳，沒辦法及時調整，必須重新施測。依據參考站數量分為兩類，單參考站及 多參考站：單參考站為使用單一固定站作為相對定位的參考點(基站)，與移動站 進行二次差分解算，亦為本研究所使用之方法；多參考站為使用多個固定站作為 相對定位的參考點，將觀測資料組成網形解算，此法有能夠穩定解算的特點，特 別是在於靜態解算時能夠相當有效的提升待測點座標的精度，然而用於動態解算 時，因待測站為持續移動變換座標位置之移動站，若要將所有移動路徑都組成網 形進行解算，便會造成資料量龐大、運算耗時等問題(陳逸如, 2010)。
21
3.3 TRACK GPS 動態差分定位解算
本研究所使用之 GPS 定位解算軟體為 GAMIT/GLOBK 下之動態定位模組
TRACK 。 GAMIT/GLOBK 為 美 國 麻 省 理 工 學 院 (Massachusetts Institute of Technology, MIT)、美國加州大學聖地雅哥分校海洋科學研究所(Scripps Institution of Oceanography, SIO)與 Harvard-Smithsonian 天文物理中心(Harvard-Smithsonian Center for Astrophysics, CfA)共同研發之高精度 GPS 科學解算軟體(Herring et al., 2008)，而 TRACK 則是麻省理工學院地球、大氣與行星科學系(Department of Earth Atmospheric and Planetary Sciences, EAPS)所發展之解算動態 GPS 數據之後處理 分析計算模組，作業環境需架構於 UNIX/Linux。所有 GAMIT、GLOBK 與 TRACK
相關使用資料與手冊皆可至其官網下載(http://www-gpsg.mit.edu/~simon/gtgk/)。
以下關於 TRACK 的介紹主要參考官網提供之文件檔並結合自行解算經驗。 TRACK 處理程序共分兩大部分，預處理與主要處理。TRACK 的執行是依 靠 一 指 令 批 次 檔 (command file) ， 預 處 理 部 分 首 先 會 需 要 讀 取 精 密 星 曆 檔 (SP3/NAV file)以取得衛星軌道和時錶資訊，接著讀取觀測資料檔(Rinex Obs. file) 處理其檔頭轉換為 TRACK 所需檔案格式並輸出 X file，接著先計算得一初步解 以作後續應用；主要處理第一步先讀取觀測資料檔以及預處理中所得之 X file， 而後先求解電碼解、週波未定值解算及週波脫落偵測、模式選定處理(AIR、 SHORT、LONG)、載波相位求解、定位成果分析(Bad data)、成果輸出，最後除 定位成果外同時也會有 TRACK 對解算成果之評估數據。處理流程如圖 3.2。
22
Read in Command
Initialize Solution Process Obs. Head File Process Sat. Information
Read Obs. Data
Code KF Process
WL Ambiguity Check Cycle Slip Check
Model
Phase KF Process
Post Data Analysis Cycle Slip, Bad Data...
Output
Ambiguity Search SP3 / NAV file
Rinex Obs. file
X file
Pre-processing
Main processing
Loop per epoch
Result & Summary of Analysis
23 TRACK 與其他 GPS 定位解算軟體的最大不同處在於它會預讀(pre-read)所有 待處理數據資料，並透過這個步驟先進行快速求解，而後將此初步解用於卡曼濾 波(Kalman Filter, KF)進行完整 GPS 動態資料解算。KF 主要原理為結合現在時刻 的觀測量以及前一時刻對現在時刻的預測值，而後整合推估現在時刻未知數的最 佳估計值，並遞迴更新估計值的變方矩陣，直至最末筆觀測資料。 動態差分的解算為使用一靜態基站座標解算待測之移動站座標，影響解算成 功率的關鍵在於基站與待測站間之基線長度。一般而言，基線小於 10 km 者解算 相當容易，因基站與移動站之大氣條件可視為相同，故其電離層誤差可設定為 0， 而後分別求解 L1 與 L2 之整數週波未定值；而在基線大於 10 km 時，基站與移 動站之大氣條件不在相同，而為了有效消除電離層延遲誤差，本研究以雙頻無電 離層線性組合(Ionosphere Free Linear Combination, L3)觀測量進行相對差分定位， 並使用軟體中 Melbourne-Wubena wide lane(MW-WL)模組進行 L1-L2 之整數週波 未定值求解。本研究選擇使用 TRACK 的“relative-rank” (RR)演算法，並應用 Chi-squared 統計檢定方式訂定最佳的 L1 及 L2 整數週波未定值之接受門檻，再 與次佳解進行比較，藉由統計測試過程決定最終的最佳解。Chi-squared 的統計 測試將由以下二種指標決定： (1) LC (Linear Combination)線性組合的解；(2) MW-WL average value 寬巷平均值之解。對流層延遲誤差部分，TRACK 採用的 大氣參數(條件)為標準大氣參數，其設定為溫度 25°C，壓力 1013.25 Mbar, 濕度 50 %。而標準大氣模式未改正完善的部分，將由對流層修正模式之附加參數吸收 (黃金維，2011)。 利用地面氣象資料(如壓力、溫度、濕度…等)代入初始的標準對流層折射誤 差修正模式，如 Modified Hopfield、Saastmoinan…等，先進行對流層折射誤差的 初步修正，而後選擇適當的歸算函數(Mapping Function)，配合附加適當數量的對 流層附加參數吸收殘留的對流層誤差，將這些附加的參數與待求解的參數(點位 坐標、速度及整數週波未定值)一起解算，可有效修正對流層折射誤差，提昇 GPS
24
高程測量的定位精度(Dodson, 1996)。本研究使用標準的 Modified Hopfield 改正 模式(Goad and Goodman, 1974)，並引用附加對流層改正參數(每個測站每 1 小時 附加 1 個改正參數)修正對流層折射之影響量，Modified Hopfield 改正模式如下： 9 , 6 1 10 k i Trop k i i k N r k α ϑ − = ∆ = ⋅ ⋅ ⋅
∑
， i=0~1 (3.6) 其中，∆ 為對流層誤差修正量，ϑ Trop i N 、α 、k i, rik分別為 Modified Hopfield Model 之使用參數。 因此，對觀測站 k 與衛星 i 而言，對流層折射延遲效應之修正量 i k ϑ ∆ 可表示 成如下形式： , , 0 , 0 , ik ( ) ( ) ( ) i i i apr k apr k z k k z k z z f z f z t ϑ = ϑ = ϑ = ∆ = ⋅ ∆ + ⋅ ∆ (3.7) , , 0 apr k z ϑ = ∆ ：為初始模式於測站k天頂方向之修正量，如果使用標準的大氣資 料(如標準氣壓、平均溫度、濕度等)，而未有實際之氣象資料，則此修正量與時 間無關，僅與測站的高程有關。 i k z ：為衛星 i 於測站 k 之天頂距。 aprf ：為初始模式的數學模式(Saastamoinen Model 或 Modified Hopfield
Model)。 , 0( ) k z t ϑ = ∆ ：為測站k在天頂方向 Mapping Function 的附加參數，與時間有關。 ( ik) f z ：為 Mapping Function 的數學函數式(例如： 1 cos( )z )。 (黃金維，2011)
25
3.4 TRACK 動態定位精度檢核分析
3.4.1 靜態測試
為確立 TRACK 解算動態差分定位之能力，以確保其後續應用於本研究之可 行性，故本研究選定 3 個一等衛星控制點進行靜態定位測試，分別為苗栗縣後龍 (HL01)、新北市石門(SHMN)及彰化縣和美(VR01)，以新竹市南寮(SHJU)為主站 進行解算，並分別就各自定位成果及組成基線長變化進行分析。所選用 GPS 資 料時間為自 2012 年 DOY 135 之 00 時 00 分 00 秒開始約 6 小時 15 分鐘的數據， 資料取樣頻率為 1 Hz。 各自定位成果乃指以所選各衛星控制點自身經 TRACK 解算所得之經度、緯 度及橢球高評估每一控制點之定位精度。將以上各站、各項目與各自平均值相減 得其殘差值(residual)後，殘差值中最大值(Max)、最小值(Min)、平均值(Mean)及 RMS 如表 3.1 所示。 表 3.2 HL01、SHMN 、VR01、 VR02 之經度(Lon)、緯度(Lat)及高(h)之殘差值 統計Station Coordinate Max Min Mean RMS
HL01
Lon (degree) 3.155E-07 -3.045E-07 -6.772E-12 9.105E-08
Lat (degree) 3.691E-07 -2.809E-07 -1.718E-12 9.805E-08
h (m) 0.400 -0.215 3.607E-14 0.099
SHMN
Lon (degree) 3.859E-07 -3.640E-07 -4.085E-12 7.681E-08
Lat (degree) 2.799E-07 -2.601E-07 -1.723E-12 9.900E-08
h (m) 0.230 -0.166 -1.813E-13 0.070
VR01
Lon (degree) 4.115E-07 -3.685E-07 6.710E-11 1.178E-07
Lat (degree) 4.932E-07 -3.768E-07 1.027E-10 1.086E-07
26 將各站之間組成基線並計算基線長，得其平均基線長後再以每秒之成果與平 均基線長相減得殘差分析如表 3.2 及圖 3.3 至 3.5 之成果。而由表 3.2 可知，在靜 態狀況下，基線長度小於等於 169 公里左右時，其 RMS 表現都有 1 公分左右之 精度。此部份成果將於下節中與船測動態成果作一綜合討論。 表 3.3 測站間組基線之基線長殘差值統計 Baseline Mean
Baseline Length (m) Max (m) Min (m) Mean (m) RMS (m)
VR01-SHMN 169115.456 0.039 -0.036 -4.546E-10 0.011 HL01-SHMN 107026.305 0.029 -0.038 -2.861E-09 0.008 HL01-VR01 65144.738 0.035 -0.033 -7.476E-11 0.011 圖 3.8 2012 年 DOY135 基線 VR01-SHMN 之殘差 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.05 1 10 57 21 13 31 69 42 25 52 81 63 37 73 93 84 49 95 05 10 561 11 617 12 673 13 729 14 785 15 841 16 897 17 953 19 009 20 065 21 121 m et er DOY 135 SEC
VR01-SHMN
VR01-SHMN27 圖 3.9 2012 年 DOY135 基線 HL01-SHMN 之殘差 圖 3.10 2012 年 DOY135 基線 HL01-VR01 之殘差 3.4.2 外部精度評估 外部精度的驗證上，本研究以全測儀量測出的接收儀座標計算得出接收天線 間的距離為比較基準。全測儀型號為 Trimble 601M，其測距精度為 2 mm + 2 ppm， 經過如公式 3.8 的計算後可得出此次施測精度： -0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0 5000 10000 15000 20000 25000 m et er DOY 135 SEC
HL01-SHMN
HL01-SHMN -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0 5000 10000 15000 20000 25000 m et er DOY 135 SECHL01-VR01
HL01-VR0128 ( ) M = ± + ×a b D (3.8) 其中 a 為儀器固定誤差，即 2 mm，主要是由儀器加常數的測定誤差、對中誤差、 測相誤差造成的，與測量的距離無關，意即不管施測的實際距離有多遠，全測儀 都將存在不大於該值的固定誤差； b D× 為比例誤差，即 2 ppm×D(km)，主要是 由儀器頻率誤差、大氣折光差所引起，D 為全測儀實際施測的距離，單位為公里， 即比例誤差會隨著實際施測距離以一定比例增加。本次實驗施測距離小於 10 公 尺，經公式(3.8)的計算所得之施測精度約為 2.02 mm。
TRACK 的動態差分定位成果方面，將 2012 年 DOY135 至 DOY139 所得定 位座標計算後得到 A、C、D 三天線相位中心之間的距離，共三個基線線段 DC、 DA、AC，其平均值與全測儀量測之值約有 0.2~2.1 cm 的差值，由表 3.3 中可發 現基線 DA 解算成果之平均值皆僅有 mm 級之差值，而基線 DC、AC 兩者解算 成果之平均值卻與全測儀所得普遍有 1cm 以上的差異，因此推測可能是接收儀 C 所得資料品質稍差於 A 和 D。與比較基準(全測儀)相減後得到之殘差分析如表 3.4，RMS 的表現介於 1.6 cm 至 3.4 cm，平均約 2.1 cm；標準差則介於 1.3 cm 至 2.6 cm，平均約 1.9 cm，由此分析可知其精度可達公分級。圖 3.6 至圖 3.8 分別 為 DOY135~DOY139 共 5 日由 TRACK 定位成果計算所得 DC、DA 及 AC 基線 長減去全測儀觀測量之殘差值。 將 3.4.1 節中靜態長基線(百公里以上)之成果與此節動態短基線的成果作一 比較分析：靜態長基線之 RMS 精度約為 1 公分，而動態短基線之 RMS 則約 2 公分，比較得短基線之成果較長基線差，這與一般長短基線解算成果之比較相反， 推估原因為長基線雖達 169 公里，但其為靜態觀測，反觀短基線卻是在船載動態 的狀況下，每一時刻皆在變動，因此使得短基線的精度表現不如長基線，也說明 動態觀測條件下，精度表現會受到影響。
29
表 3.4 全測儀基線長、TRACK 平均基線長、兩者差值
Baseline DOY Total Station (m) TRACK mean value (m) Difference (m)
DC 135 1.989 2.007 0.018 136 2.004 0.015 137 2.002 0.013 138 1.999 0.010 139 2.000 0.011 DA 135 3.312 3.310 -0.002 136 3.309 -0.003 137 3.307 -0.004 138 3.307 -0.005 139 3.308 -0.004 AC 135 2.075 2.054 -0.021 136 2.057 -0.018 137 2.056 -0.019 138 2.080 -0.012 139 2.060 -0.015 表 3.5 TRACK 基線長與全測儀所得基線長差異值統計
DOY Baseline Max Min Mean Std. RMS
DOY135 DC 0.110 -0.053 0.018 0.013 0.022 DA 0.062 -0.086 0.002 0.018 0.018 AC 0.092 -0.125 -0.021 0.021 0.030 DOY136 DC 0.149 -0.083 0.015 0.014 0.020 DA 0.171 -0.131 -0.003 0.016 0.022 AC 0.246 -0.178 -0.018 0.023 0.025 DOY137 DC 0.082 -0.076 0.013 0.014 0.015 DA 0.114 -0.103 -0.004 0.017 0.021 AC 0.090 -0.155 0.019 0.023 0.024 DOY138 DC 0.097 -0.078 0.016 0.015 0.020 DA 0.066 -0.109 0.005 0.019 0.020 AC 0.090 -0.145 0.005 0.026 0.034 DOY139 DC 0.058 -0.067 0.011 0.014 0.017 DA 0.069 -0.073 -0.004 0.015 0.016 AC 0.064 -0.104 -0.015 0.018 0.024
30 圖 3.11 TRACK 定位成果計算所得 DC 基線長與全測儀觀測值之差異。 由上至下為 DOY135~DOY139。 -0.1 -0.05 0 0.05 0.1 0.15 40000 50000 60000 70000 80000 m et er DOY135 SEC
DOY135 DC
Diff. -0.1 0 0.1 0.2 0 20000 40000 60000 80000 m et er DOY136 SECDOY136 DC
Diff. -0.1 -0.05 0 0.05 0.1 0 10000 20000 30000 40000 m et er DOY SECDOY 137 DC
Diff. -0.12 -0.08 -0.04 0 0.04 0.08 0.12 25000 40000 55000 70000 85000 m et er DOY SECDOY 138 DC
Diff. -0.05 0 0.05 0.1 25000 27000 29000 31000 33000 35000 m et er DOY139 SECDOY139 DC
Diff.31 圖 3.12 TRACK 定位成果計算所得 DA 基線長與全測儀觀測值之差異。 由上至下為 DOY135~DOY139。 -0.1 -0.05 0 0.05 0.1 40000 50000 60000 70000 80000 m et er DOY135 SEC
DOY135 DA
Diff. -0.2 -0.1 0 0.1 0.2 0 20000 40000 60000 80000 m et er DOY136 SECDOY136 DA
Diff. -0.2 -0.1 0 0.1 0.2 0 10000 20000 30000 40000 m et er DOY SECDOY 137 DA
Diff. -0.15 -0.1 -0.05 0 0.05 0.1 25000 40000 55000 70000 85000 m et er DOY SECDOY 138 DA
Diff. -0.1 -0.05 0 0.05 0.1 25000 27000 29000 31000 33000 35000 m et er DOY139 SECDOY139 DA
Diff.32 圖 3.13 TRACK 定位成果計算所得 AC 基線長與全測儀觀測值之差異。 由上至下為 DOY135~DOY139。 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 40000 50000 60000 70000 80000 m et er DOY135 SEC
DOY135 AC
Diff. -0.2 -0.1 0 0.1 0.2 0.3 0 20000 40000 60000 80000 m et er DOY136 SECDOY136 AC
Diff. -0.2 -0.1 0 0.1 0.2 0 10000 20000 30000 40000 m et er DOY SECDOY 137 AC
Diff. -0.2 -0.1 0 0.1 0.2 25000 40000 55000 70000 85000 m et er DOY SECDOY 138 AC
Diff. -0.15 -0.1 -0.05 0 0.05 0.1 25000 27000 29000 31000 33000 35000 m et er DOY139 SECDOY139 AC
Diff.33
3.4.3 TRACK 內部精度評估
本研究將 4 組 GPS 接收儀安裝於船載重力測量的船上，取樣頻率為 1Hz。 解算時以連續航線為計算單位，以獲得每連續航線的 GPS 動態橢球高資訊及其 精度值，並輸出估計天頂方向的對流層附加參數值與訊號延遲誤差量，並藉此延 遲誤差判斷 TRACK 動態解算成果之合理性。 對流層延遲誤差由空氣中的乾延遲量與溼延遲量共同組成，乾延遲量由乾燥 的氣體所引起，其在天頂方向造成的延遲量約為 2.3 m，約占對流層總延遲量的 90%，乾空氣之密度會隨其表面壓力的改變而變化，由於其變化較為平緩穩定， 故可藉由量測大氣的溫度、溼度、壓力將此乾延遲量予以模式化並求解；而溼延 遲量則與大氣的水氣分壓有關，由於水氣會隨著時間與空間高度之改變而變化， 因此不易利用數學模式予以準確估計，在溼度較高的地區，溼延遲量可高達 35 cm(Bevis et al., 1992)，故若求解得的對流層總延遲量介於 2~3 m 之間、平面定位 精度介於 2~4 cm、高程定位精度介於 5~10 cm 之間，即可視該計算成果為合理 的結果(陳逸如, 2010)。此節以 2012 年 DOY136 之 TRACK 解算成果為例。 圖 3.14 TRACK 解算 2012 DOY136 之橢球高、天頂距對流層延遲量及軟體內部 精度34
3.5 改正模式
3.5.1 海潮改正模式
NAO.99b 及 NAO.99Jb 為 National Astronomical Observatory in Japan 研製短 週期海潮模型，分別為全球性以及區域性海潮模型。NAO.99b 涵蓋範圍由經度 0°-360°；北緯 90°-南緯 90°，解析度為 0.5 度，它結合了五年的 TOPEX/POSEIDON 測高儀之海潮資料(cycle 9-198)以及水文動力模型；NAO.99Jb 為一日本區域性海 潮模型，以 NAO.99b 結合 TOPEX/POSEIDON 資料，提供 16 個分潮：2N2, J1, K1, K2, L2, M1, M2, MU2, N2, NU2, O1, OO1, P1, Q1, S2, T2。在日本附近海域(東經 110°~165°；北緯 20°~65°，包含臺灣地區)結合日本沿岸驗潮站同化資料以發展 NAO.99Jb 模式，其空間解析度為 5 分，其中提供 16 個分潮調和常數，並以 Matsumoto et al.編撰的電腦程式完成 16 個分潮調和常數的載入(Matsumoto et al., 2000)。本研究使用 NAO.99Jb 模型進行海潮改正。然於不同的參考文獻當中發 現與岸邊驗潮站所蒐集之潮位資料結合，也是常見的海潮改正方法之一(Foster, 2009; Bonnefond, 2003)。
3.5.2 固體潮改正模式
地球表面會因太空中攝動天體(如月球、太陽)對彈性地球的引力作用而產生 週期性漲落的現象，稱為固體潮現象。此現象之概略描述為在地球地心與攝動天 體之連線方向上會被拉長，而連線垂直方向則趨於扁平，其對地表面之影響包含 長期性偏移以及週期性偏移，前者與緯度有關，後者由日週期和半日週期組成。 經由一天的靜態觀測，週期性偏移之影響大部分都可平滑消除，可觀測到的位移 量可達約 40 公分(Lambeck, 1988)，而長期性偏移方面，即使進行長時間觀測， 殘餘誤差最多仍可達 10 公分左右。在相對定位方法中，一百公里以內的基線二 測站固體潮影響可視為一致，可以差分方法消除。35
3.5.3 濾波
本研究所使用 Gaussian filter 對研究中所有計算出的海面高資料進行濾波， 而此濾波主要是針對船航行時的風浪，而非船速可能產生的誤差，原因是一般船 行速度較高的情況下，對海水面高訂定有 5~10 公分甚至更高的誤差可能產生， 但 當 船 速 低 於 3 m/s 情 況 時 ， 行 進 速 度對 海 面 高 訂 定 的 影 響小 於 一 公 分 (Bonnefond, 2003)，本研究航行船速控制在 3 m/s 左右，故船速對海面高影響忽 略不計。本研究使用之 Gaussian filter (Hwang et al., 2006)是以離選定資料越近而 權重越重的加權方法平滑化所有資料，以期平滑化後能夠消除船行進間姿態角因 風浪影響而擺幅過大等高頻誤差。由前人實驗發現 120 秒罩窗所得海水面高精度 比 60 秒罩窗於平差後 RMS 值改善約一公分(陳逸如, 2010)，而本次研究中亦由 無濾波、60 秒、120 秒和 360 秒罩窗的成果發現 120 秒罩窗成果在本研究中表現 最好，於是選定 120 秒罩窗為處理本次研究的罩窗大小。36
第四章 GPS 與 INS/GPS 資料計算
4.1 研究區域簡介
4.1.1 施測範圍、儀器、測量船隻簡介
本研究為臺灣船載重力計畫之一部分，此為一 2011 年起之 4 年期之計畫， 研究中所使用之研究數據為 2012 年 5 月至 6 月期間所蒐集，施測範圍為臺灣西 部近岸 12 海里(約 19.6 km)海域，涵蓋區域南起台中梧棲漁港、北迄新北市富基 漁港，規劃航線如圖 4.1 所示，實際施測航線以 TRACK 動態差分解算成果繪製， 如圖 4.2 所示。37
圖 4.15 2012 年近岸船載重力計畫航線規劃圖(摘自內政部 101 年度臺灣本島近岸 船載重力計畫書)
38 圖 4.16 2012 年近岸船載重力計畫實際施測航線圖 在考量計畫施測範圍、航線及航期規畫後，針對船隻續航力、連續作業能力、 船體穩定度及作業安全等因素，計畫所選用之船隻為 CT2 等級之尖再發號，船 隻規格表如表 4.1，船隻外觀如圖 4.3。多天線 GPS 系統方面，所使用之儀器為 大地測量型 GPS 接收儀 Trimble 5700，天線盤型號 TRM42149.00，架設位置如 圖 4.4 中藍色圓圈及標記 A、B、C、D 部分；GPS/INS 整合系統方面，慣性測量 儀(Inertial Measurements Unit, IMU)的型號為 LITEF LCI 戰術等級 IMU，並搭配 架設 NovAtel Propak 雙系統雙頻 GNSS 大地等級接收儀，架設位置如圖 4.4 兩個
39 綠色圓圈部分。 表 4.6 船載重力計畫選用船隻─尖再發號 項目 內容 船長(註冊尺度) 10.14 公尺 船寬(註冊尺度) 3.40 公尺 總噸位 15.75 公噸 淨噸位 4.73 公噸 油櫃容量 1600 公升 馬力 279 KW （380 PS） 船殼材質 玻璃纖維強化塑膠 圖 4.17 尖再發號船體外觀
40
圖 4.18 上為實驗儀器於船隻上架設位置圖、下為全測儀施測數據圖，紅色圓圈 及 A、B、C、D 為 GPS 接收天線位置，2 綠色圓圈為 IMU 及其 GNSS 天線位置。
41
