GLONASS衛星軌道演算及其定位方法

國立臺灣大學電機資訊學院電機工程學系 碩士論文

Department of Electrical Engineering

College of Electrical Engineering and Computer Science

National Taiwan University Master Thesis

GLONASS 衛星軌道演算及其定位方法 GLONASS Orbit Determination and

Positioning Method

許庭瑋 Ting-Wei Hsu

指導教授：張帆人 博士 王立昇 博士

Advisors：Fan-Ren Chang, Ph.D. Li-Sheng Wang, Ph.D.

中華民國 102 年 7 月

July 2013

i

誌 誌 誌

誌 謝 謝 謝 謝

兩年的碩士生涯，說長不長說短不短，期間的收穫與回憶不勝枚舉，然而完 成此碩士論文絕對是當中最難忘的里程碑。

首先，我要感謝指導教授張帆人博士與王立昇博士，從碩零暑假開始就賦予 我極具挑戰性的研究計畫。每兩週一次的富士通討論會議著實讓我學習與成長許 多，也大大提升我學理與實務應用雙方面之能力，同時，對於此篇論文的撰寫也 有莫大的幫助。除此之外，老師面對生活的態度與應對進退的能力更是身教勝於 言教，永遠不會忘記老師的幽默風趣與處事謹慎。

其次，我要感謝本實驗室與隔壁實驗室所有的學長姐、同學與學弟妹，以及 應力所的同學們。感謝弘唯學長於實驗上的幫助及生活休閒上的領導，感謝沛宏 學長不時的指點與建議，感謝逸群學長平時的閒談與人生方向的引導，感謝政昌 學長最後的畫龍點睛，感謝炳輝學長於修課時的協助，感謝政達學長每週的建議 與討論，感謝致豪學長的陪伴與建議，感謝宜璋學長與椀潔學姊處處幫忙。感謝 政憲、培耕、雨弘、澤民、敏琦與進業，我們一起努力一起奮鬥，一起成功走到 最後。感謝志陽學弟、勁璋學弟、兆慶學弟與張時中教授的學生們，因為你們使 兩年的實驗室生活過得更精彩更開心。能與大家一同度過這歡樂的碩士生涯是我 的榮幸，我將不會忘記你們。

最後，我要感謝我的父母、家人與朋友。在課餘時間你們總是陪伴著我，在 不順心的時候你們總是默默支持與鼓勵，這些都是讓我向前最大的動力，你們缺 一不可。謝謝有你們，讓我能夠度過難關並且順利完成學業。

摘 摘 摘

摘 要 要 要 要

全球導航衛星系統(Global Navigation Satellites System; GNSS)乃泛指所有以衛 星為基礎之導航系統。目前發展較為成熟之導航衛星系統包括美國的 GPS(Global Positioning System)以及俄羅斯的 GLONASS(Global Navigation Satellites System)。

隨著科技的演進，不同領域的應用對於導航與定位之精度需求越加精細。提 升定位精度最直接之方法即為提升可視衛星數目，而整合 GPS 及 GLONASS 可達 此目的。由於 GPS 及 GLONASS 兩系統不可完全兼容，在整合前需對各別系統有 相當程度之了解。因此，吾人對 GLONASS 進行研究，嘗試運用該系統作定位，

並觀察其結果。

在本文中，吾人從能量之觀點推導 GLONASS 衛星軌道模型，找出衛星運動 方程式，並將該系統之廣播星曆代入求解，藉龍格-庫塔數值積分法取得準確之衛 星位置。其後，利用衛星位置與偽距於偽距觀測量定位法進行定位。偽距觀測量 定位法中，需修正電離層延遲誤差，才能使定位精度提高，而 GLONASS 並無廣 播電離層模型之參數予以修正誤差，故吾人將 GPS 廣播之電離層模型參數應用於 GLONASS 定位之上。最後，引入高度限制條件式於定位方程中，而定位精度與穩 定度均有相當改善。

本文除顯示運用 GLONASS 定位之結果外，也詳細描述整個定位過程。從接 收到廣播星曆開始直至定位完成，期間計算所需注意之部分皆於本文當中。此舉 將便於未來相關應用。

關鍵詞：GLONASS 衛星軌道模型、龍格-庫塔數值積分法、偽距觀測量定位法、

電離層模型修正法、高度限制條件式

iii

Abstract

In general, GNSS (Global Navigation Satellites System) is described as the satellites-based navigation system. Presently, the available systems include GPS (Global Positioning system) of United States and GLONASS (Global Navigation Satellites System) of Russia.

With more available satellites, the accuracy of positioning could become more precise, and the combination of GPS and GLONASS could easily reach it. Although GPS and GLONASS are similar to each other, the both systems are not entirely compatible. Before combining these two systems, it is important to be familiar with each system in great details. Therefore, we made a study of GLONASS and attempted to specify the position by means of only GLONASS.

In this paper, we determine the satellite position through GLONASS ephemeris by using Runge-Kutta Method, and then obtain the receiver position with pseudorange measurement. To correct ionosphere delays on GLONASS satellites, we applied Klobuchar Model with little modification. Afterward, we also applied altitude-hold constraint to positioning equations, and consequently, the results of positioning point have become more precise and stable.

It shall be noted that we not only demonstrated the results but also elaborated on the whole process of positioning. This will benefit the future studies and some relative applications.

Key Words: GLONASS Orbit Model, Runge-Kutta Method, Pseudorange Measurement, Klobuchar Model, Altitude-Hold Constraint

目 目 目

目 錄 錄 錄 錄

誌謝 誌謝 誌謝

誌謝 ... i

摘要 摘要 摘要 摘要 ... ii

Abstract ... iii

目錄目錄 目錄目錄 ... iv

圖目錄 圖目錄 圖目錄 圖目錄 ... vi

表目錄 表目錄 表目錄 表目錄 ... ix

第一章第一章 第一章第一章 緒論緒論緒論緒論 ... 1

1.1 研究動機 ... 1

1.2 研究方向 ... 2

1.3 論文架構 ... 3

第二章 第二章 第二章 第二章 GLONASS 全球導航衛星系統概述全球導航衛星系統概述全球導航衛星系統概述全球導航衛星系統概述 ... 5

2.1 系統架構與衛星種類 ... 5

2.1.1 太空部分 ... 6

2.1.2 地面控制部分 ... 6

2.1.3 使用者部分 ... 7

2.1.4 衛星種類 ... 8

2.2 衛星訊號與導航訊息 ... 9

2.2.1 載波頻率之分配 ... 10

2.2.2 虛擬隨機測距碼 ... 10

2.2.3 導航訊息之產生 ... 12

2.2.4 導航訊息內容 ... 15

2.2.5 導航訊息架構 ... 19

2.3 座標系統與時間系統 ... 22

2.3.1 座標系統類別 ... 23

2.3.2 PZ-90.02 座標系統 ... 24

v

2.3.3 UTC(SU)時間系統 ... 25

第三章第三章 第三章第三章 GLONASS 衛星軌道模型衛星軌道模型衛星軌道模型衛星軌道模型 ... 27

3.1 衛星軌道方程式之推導 ... 27

3.2 數值積分法 ... 33

3.3 演算之結果 ... 36

3.3.1 與網站數據相較之誤差分析 ... 36

3.3.2 與接收機相較之誤差分析 ... 39

第四章 第四章 第四章 第四章 GLONASS 定位原理與方法定位原理與方法定位原理與方法定位原理與方法... 43

4.1 定位基本原理 ... 43

4.2 單機定位與偽距觀測量 ... 44

4.3 電離層延遲誤差 ... 50

4.3.1 訊號傳遞之特性與概念 ... 50

4.3.2 電離層模型修正法 ... 53

4.3.3 雙頻觀測量修正法 ... 56

4.4 高度限制條件式 ... 57

第五章 第五章 第五章 第五章 實驗結果與分析實驗結果與分析實驗結果與分析實驗結果與分析 ... 63

5.1 實驗設備簡介 ... 63

5.2 實驗方法介紹 ... 63

5.3 單機靜態定位實驗結果 ... 66

第六章第六章 第六章第六章 結論與未來展望結論與未來展望結論與未來展望結論與未來展望 ... 87

6.1 結論 ... 87

6.2 未來展望 ... 87

參考參考 參考參考文獻文獻文獻文獻 ... 89

圖目錄 圖目錄 圖目錄 圖目錄

圖 2.1 GLONASS 系統架構 ... 5

圖 2.2 GLONASS 衛星分佈圖 ... 7

圖 2.3 移位暫存器之結構 ... 11

圖 2.4 虛擬隨機距離碼與時脈產生之簡示方塊圖 ... 11

圖 2.5 資料序列生成之簡示方塊圖 ... 12

圖 2.6 導航訊息與虛擬隨機測距碼之二進位串列脈衝同步關係圖 ... 14

圖 2.7 導航訊息生成圖 ... 15

圖 2.8 超訊框架構 ... 20

圖 2.9(a) 第一至第四訊框架構 ... 21

圖 2.9(b) 第五訊框架構 ... 21

圖 2.10 資料串架構 ... 22

圖 2.11 地心慣性座標系統與地心地固座標系統關係示意圖 ... 23

圖 3.1 衛星 48 個時刻之位置誤差 ... 37

圖 3.2 衛星 48 個時刻之速度誤差 ... 37

圖 3.3 24 顆衛星第 2 時刻之位置誤差 ... 38

圖 3.4 衛星位置與傳遞訊號之時間描述 ... 40

圖 3.5 模擬之衛星位置與 SATXYZ data 之誤差 ... 41

圖 4.1 定位基本原理示意圖 ... 43

圖 4.2 衛星與接收機位置關係示意圖 ... 44

圖 4.3 訊號傳遞時間圖 ... 46

圖 4.4 電子總量(TEC)示意圖 ... 52

圖 4.5 經緯高座標與 ECEF 直角座標關係圖 ... 58

圖 4.6 三角形 OPQ 示意圖 ... 59

圖 5.1 定位演算法流程圖 ... 65

圖 5.2 各時刻之 GLONASS 衛星數目 ... 66

vii

圖 5.3(a) Calculate data/PSRPOS data 水平定位結果 ... 67

圖 5.3(b) Calculate data 各時刻之 ENU 座標 ... 68

圖 5.3(c) Calculate data 與參考點之水平與垂直誤差... 68

圖 5.3(d) Calculate data/PSRPOS data ENU 座標 E 方向 ... 69

圖 5.3(e) Calculate data/PSRPOS data ENU 座標 N 方向 ... 69

圖 5.3(f) Calculate data/PSRPOS data ENU 座標 U 方向 ... 70

圖 5.3(g) Calculate data 與 PSRPOS data 三方向之誤差 ... 70

圖 5.4(a) Calculate data^h/PSRPOS data 水平定位結果 ... 71

圖 5.4(b) Calculate data^h各時刻之 ENU 座標 ... 72

圖 5.4(c) Calculate data^h與參考點之水平與垂直誤差 ... 72

圖 5.4(d) Calculate data^h/PSRPOS data ENU 座標 E 方向 ... 73

圖 5.4(e) Calculate data^h/PSRPOS data ENU 座標 N 方向 ... 73

圖 5.4(f) Calculate datah/PSRPOS data ENU 座標 U 方向 ... 74

圖 5.4(g) Calculate datah 與 PSRPOS data 三方向之誤差 ... 74

圖 5.5 各時刻之 GLONASS 衛星數目 ... 75

圖 5.6(a) Calculate data/PSRPOS data 水平定位結果 ... 76

圖 5.6(b) Calculate data 各時刻之 ENU 座標 ... 77

圖 5.6(c) Calculate data 與參考點之水平與垂直誤差... 77

圖 5.6(d) Calculate data/PSRPOS data ENU 座標 E 方向 ... 78

圖 5.6(e) Calculate data/PSRPOS data ENU 座標 N 方向 ... 78

圖 5.6(f) Calculate data/PSRPOS data ENU 座標 U 方向 ... 79

圖 5.6(g) Calculate data 與 PSRPOS data 三方向之誤差 ... 79

圖 5.7(a) Calculate data^h/PSRPOS data 水平定位結果 ... 80

圖 5.7(b) Calculate data^h各時刻之 ENU 座標 ... 81

圖 5.7(c) Calculate data^h與參考點之水平與垂直誤差 ... 81

圖 5.7(d) Calculate data^h/PSRPOS data ENU 座標 E 方向 ... 82

圖 5.7(e) Calculate data^h/PSRPOS data ENU 座標 N 方向 ... 82 圖 5.7(f) Calculate datah/PSRPOS data ENU 座標 U 方向 ... 83 圖 5.7(g) Calculate datah 與 PSRPOS data 三方向之誤差 ... 83

ix

表 表 表 表目錄 目錄 目錄 目錄

表 2.1 GLONASS 與 GPS 之差異 ... 6

表 2.2 GLONASS 衛星之規格與特性 ... 8

表 2.3 GLONASS-M 衛星之規格與特性 ... 9

表 2.4 GLONASS-K 衛星之規格與特性 ... 9

表 2.5 GLONASS 衛星載波頻率分佈於 L1 與 L2 ... 13

表 2.6 GLONASS 衛星座標值與速度精度 ... 16

表 2.7 P1 旗標內容 ... 16

表 2.8 廣播星曆參數之字元特性 ... 17

表 2.9 FT 旗標內容 ... 18

表 2.10 GLONASS 非即時性資料於超訊框內之分配 ... 19

表 2.11 PZ-90.02 座標系統之常數與相關參數 ... 24

表 3.1 各擾動力對 GPS 衛星之影響 ... 29

表 3.2 SATXYZ data 內容參考範例 ... 39

表 5.1 單機靜態定位參考座標 ... 64

表 5.2 Calculate data/PSRPOS data 於 ENU 座標之平均值與標準差 ... 67

表 5.3 Calculate data^h/PSRPOS data 於 ENU 座標之平均值與標準差 ... 71

表 5.4 Calculate data/PSRPOS data 於 ENU 座標之平均值與標準差 ... 76

表 5.5 Calculate data^h/PSRPOS data 於 ENU 座標之平均值與標準差 ... 80

表 5.6 第一筆數據 Calculate data/Calculate data^h/PSRPOS data 之比較 ... 84

表 5.7 第二筆數據 Calculate data/Calculate data^h/PSRPOS data 之比較 ... 84

第一章 第一章 第一章

第一章 緒論 緒論 緒論 緒論

1.1 研究動機研究動機 研究動機研究動機

太空科技自問世以來即成為各國矚目之發展重點，而衛星導航定位技術就屬 其中極為重要之一環。衛星導航定位技術之發展，起源於軍事需求，舉凡飛彈及 戰鬥機、艦等之導航追蹤、地面部隊派遣定位及眾多軍事應用等皆倚賴於該技術，

衛星導航定位儼然成為世界各強國必備之現代化科技。

冷戰期間，蘇聯與美國兩大政府於太空探索領域的競爭十分激烈。西元 1957 年 10 月 4 日，蘇聯成功發射世上第一顆人造衛星至太空，其與美國太空軍事競賽 的序幕就此展開。期間蘇聯與美國各自發展 GLONASS 與 GPS 衛星導航系統，雖 然蘇聯較美國早先一步於航太領域方面取得突破，但其後蘇聯解體，種種因素以 致 GLONASS 之發展遲滯，美國最終成為太空科技的領銜者，GPS 因而廣為眾人 所知也較多為採用。

西元 1991 年，俄羅斯聯邦政府接續 GLONASS 計劃之開發，直至西元 1995 年終將該系統之衛星數目提升至 24 顆，實現其完整佈署[1]。往後幾年，俄羅斯爆 發經濟危機，航太領域之預算因而遭大幅刪減，政府無力維持 GLONASS 之完整 佈署，故該系統之發展再度受阻。爾後，隨著俄羅斯的經濟好轉，以及於該政府 之政策規劃下，GLONASS 成為重點發展項目，新型的導航衛星 GLONASS-M 與 GLONASS-K 相繼誕生。發展至今，GLONASS 之衛星佈署達運作衛星 24 顆、備 用衛星 4 顆，其導航範圍已可覆蓋全球表面及近地高空，其精度可達數米左右。

事實上，衛星導航科技與定位技術之發展已相當成熟，以 GPS 為基礎架構之 技術發展尤是。雖然目前已可單獨使用 GLONASS 達到導航及定位之目的，但由 於 GPS 發展較早且已完備，單以 GLONASS 為基礎架構之技術應用難以進入廣大 市場，取而代之的是 GPS 與 GLONASS 兩者整合之系統。隨著科技的演進，不同 領域的應用對於導航與定位之精度需求越加精確，以往車輛船舶之導航精度可以

2

是公尺等級，而各種建築基線與地形之測量則需公分或是公厘等級之精度，不同 系統之整合應用對於提升定位之性能是有幫助的。

定位系統之整合實屬不易，時間系統之別、空間座標系統之別以及其餘系統 特性之差異等，於整合時皆需審慎考量。進行系統整合前，須對各個系統有完整 且深入之研究，如此才能熟悉各個系統之基本特性，進而由基礎出發，找出整合 系統之優良且合適的技術發展。有關 GPS 之研究已屬廣博，而 GLONASS 相關之 研究相較之下則為鮮少許多，故於整合 GPS 與 GLONASS 兩系統前，吾人對 GLONASS 進行相當程度之研究，對該系統特性有所掌握並在此予以分享。

1.2 研究方向研究方向 研究方向研究方向

基於實現單用 GLONASS 定位之理念，本論文以 GLONASS 衛星所傳遞之廣 播星曆資料(Ephemeris Data)，透過數值積分龍格-庫塔法(Runge Kutta Method)計算 出各個時刻 GLONASS 衛星的位置，配合接收機所算出之偽距(Pseudorange)進行定 位。本論文之研究主要分為兩大階段，第一個階段係由廣播星曆資料計算出各個 衛星於各個時刻之精確位置，第二階段接續由計算出之精確衛星位置配合偽距進 行定位，並引入高度限制條件式，求得接收機所在位置之最佳解。

從接收衛星訊息至最終定位完成的過程中，大部分皆使用 GLONASS 衛星傳 遞之訊息即可，唯獨電離層延遲量(Ionospheric Delay)。電離層延遲量之估測可分 為兩種方式，一種屬單頻接收機所使用，另一種則為雙頻接收機可用。單頻接收 機估測電離層延遲量需使用電離層模型，而電離層模型之相關參數 GLONASS 並 無廣播，需由 GPS 廣播之訊息取得；雙頻接收機估測電離層延遲量則可直接由雙 頻之偽距求出，不需 GPS 廣播之訊息，故以 GLONASS 單獨定位較適合使用雙頻 偽距觀測量以修正電離層之影響。由於考量整合 GPS 與 GLONASS 雙系統為主要 趨勢，本論文於 GLONASS 定位時使用電離層模型之修正法(依賴 GPS 廣播訊息) 取得電離層延遲量。

1.3 論文架構論文架構 論文架構論文架構

本論文共分為六章，茲將各章之要點簡述如下：

第一章為緒論，簡述本論文之研究動機、研究方向與論文架構。

第二章為 GLONASS 全球導航衛星系統概述，簡介 GLONASS 全球導航衛星 系統，並比較其與 GPS 之異同。首先介紹 GLONASS 系統架構，再者會提及其導 航訊息內容與載波頻率之分配，最後說明 GLONASS 所使用之座標系統與時間系 統。

第三章為 GLONASS 衛星軌道模型，為求得衛星各個時刻之精確位置，必須 先行建立衛星軌道模型。GLONASS 之衛星位置係由直角坐標軌道方程式經數值積 分龍格-庫塔法(Runge Kutta Method)積分求得。因未來 GLONASS 將與 GPS 作結合，

於求得衛星位置後需經座標系統轉換，轉換至 GPS 所使用之座標系統。最後將求 得之衛星位置與官方公布之衛星位置作對照與分析。

第四章為 GLONASS 定位原理與方法，介紹定位基本原理，進而推廣至衛星 系統定位方式，並且探討訊號傳遞之誤差的修正方式。訊號傳遞之誤差絕大部分 來自電離層延遲量的影響，故吾人在此深入探討電離層延遲量及其修正方法。最 後，提出高度限制條件式，並將其應用於定位方程式中。

第五章為實驗結果與分析，使用合適之接收機取得衛星廣播星曆訊息，接著 利用第三章之衛星軌道模型求得衛星位置，再配合接收機相關計算出之偽距，透 過第四章之定位方法求得接收機位置之最佳解，最後呈現定位結果並與接收機提 供之值及實際位置作比較。

第六章為結論與未來展望，總結本論文之研究成果並探討未來可進行之研究 方向。

4

第二章 第二章 第二章

第二章 GLONASS 全球導航衛星系統概述 全球導航衛星系統概述 全球導航衛星系統概述 全球導航衛星系統概述

GLONASS(Global Navigation Satellite System)是一個全球性、全天候之三維空 間即時定位系統，提供定位及時間校準功能予陸海空甚至是近地太空之使用者，

對於使用者的數量並無限制。起初，GLONASS 是以國防軍事為主要目的，大多用 於飛機、飛彈、船艦等之導航與時間校準，爾後俄羅斯政府漸漸開放該系統給民 間使用，民用信號於國內外皆可免費且無限制地使用。

GLONASS 功能架構與 GPS 相似，但略有相異之處。俄羅斯位於高緯度且該 系統原先設計目的為提供國防軍事方面之服務，因此 GLONASS 衛星軌道仰角偏 高，衛星星群之設置主要涵蓋高緯度地區，與美國 GPS 有所不同。在訊號傳遞方 面，GLONASS 採用分頻多工(Frequency Division Multiple Access; FDMA)，而 GPS 則採用分碼多工(Code Division Multiple Access; CDMA) ，詳情如 2.2 節所述。在 座標系統與時間系統方面，GLONASS 採用 PZ-90.02 座標系統與 UTC(SU)時間系 統，兩者皆與 GPS 相異。GLONASS 與 GPS 之差異，見表 2.1。

2.1 系統架構與衛星種類系統架構與衛星種類 系統架構與衛星種類系統架構與衛星種類

GLONASS 之系統架構，主要由太空部分、地面控制部分與使用者部分等三

大結構組成[2]，參閱圖 2.1。

圖 2.1 GLONASS 系統架構[2]

6

表 2.1 GLONASS 與 GPS 之差異

系統 GLONASS GPS

衛星數目 24+ 24+

軌道面數目 3 6

每個軌道上的衛星數 8 4

軌道面傾角 64.8° 55°

軌道半徑 25510 公里 26560 公里

週期 約 11 小時 15 分 約 11 小時 58 分

訊號調變技術 FDMA CDMA

載波頻率

L1 約 1602 MHz L2 約 1246 MHz

L1 = 1575.42 MHz L2 = 1227.60 MHz

時間系統 UTC(SU) UTC(USNO)

座標系統 PZ-90.02 WGS84

2.1.1 太空部分

GLONASS 太空部分目前係由 24 顆運作衛星與 4 顆備用衛星所組成之衛星星 群，衛星分佈如圖 2.2 所示。24 顆運作衛星平均分布於三個軌道面上，各個軌道 面之間相互間隔 120°，軌道面相對於地球赤道呈現 64.8°之傾斜角。每一個軌道面 上有 8 顆衛星，彼此以間隔 45°均勻分佈。每顆衛星運行於距地球表面 19100 公里 高之橢圓軌道，其運行週期約為 11 小時又 15 分鐘。衛星的空間分佈如此，能保 證全球表面與近地太空無時無刻皆涵蓋於此系統之下。

2.1.2 地面控制部分

地面控制部分包含系統控制中心(System Control Center)，以及遍佈於俄羅斯境

內之指揮與追蹤站(Command and Tracking Stations)間彼此連結成的網域。其主要功 能在於監控 GLONASS 衛星星群狀態，並且校正衛星軌道參數，以及上載導航相 關資料。

2.1.3 使用者部分

使用者部分包含接收機與處理器，負責接收與處理 GLONASS 之導航訊號，

並且計算使用者之位置座標，進而提供速度與時間等相關資訊。

圖 2.2 GLONASS 衛星分佈圖[3]

8

2.1.4 衛星種類

GLONASS 衛星共有三種，分別為 GLONASS、GLONASS-M 與 GLONASS-K，

目前 24 顆運作衛星皆為 GLONASS-M。各種衛星之簡介及規格如下[3]：

(1) GLONASS：

提供導航訊息及精確的時間訊號，使軍民之陸海空各種用戶皆可使用，其規 格與特性如表 2.2。截至目前為止已發射共 81 顆。

(2) GLONASS-M：

同樣提供導航訊息及精確的時間訊號，使軍民之陸海空各種用戶皆可使用，

值得注意的是 GLONASS-M 衛星傳送兩種民用訊號，使定位能夠更加精確，其規 格與特性如表 2.3。截至目前為止已發射共 26 顆。

(3) GLONASS-K：

同樣提供導航訊息及精確的時間訊號，使軍民之陸海空各種用戶皆可使用。

新的訊號傳遞技術誕生，且新增分碼多工(Code Division Multiple Access; CDMA) 訊號調變技術之訊號，可與 GPS/Compass/Galileo 相容。GLONASS-K 衛星完全由 俄羅斯打造，無使用任何進口設備，其規格與特性如表 2.4。

表 2.2 GLONASS 衛星之規格與特性[3]

壽命 3 年

質量 1415 公斤

發射訊號種類

L1：標準精度(民用)、高精度(軍用) L2：高精度(軍用)

耗電量 1000 瓦

方向定位精度 0.5°

太陽能板指向精度 5°

表 2.3 GLONASS-M 衛星之規格與特性[3]

壽命 7 年

質量 1415 公斤

發射訊號種類

L1：標準精度(民用)、高精度(軍用) L2：標準精度(民用)、高精度(軍用)

耗電量 1450 瓦

表 2.4 GLONASS-K 衛星之規格與特性[3]

壽命 10 年

質量 875 公斤

發射訊號種類

L1：標準精度(民用)、高精度(軍用)，FDMA L2：標準精度(民用)、高精度(軍用)，FDMA L3：標準精度(民用)、高精度(軍用)，FDMA

L1CR (1575.42 MHz)：CDMA L5R (1176.45 MHz)：CDMA

耗電量 1600 瓦

2.2 衛星訊號與導航訊息衛星訊號與導航訊息 衛星訊號與導航訊息衛星訊號與導航訊息

GLONASS 衛星訊號係藉由射頻(Radio Frequency; RF)之無線電波傳遞，所使

用之頻段於 L 頻帶範圍內，與 GPS 相似。目前有兩個主要載波頻段，即 1.6GHz 之 L1 頻段與 1.2GHz 之 L2 頻段。傳遞之導航訊號有兩種精度，標準精度(Standard

10

Precision)與高精度(High Precision)，其中標準精度導航訊號之時鐘振盪頻率為 0.511MHz，可供國內外商業用戶及一般民眾使用，而高精度導航訊號之時鐘振盪 頻率為 5.11MHz，僅提供予俄羅斯官方授權之單位使用[4]。

GLONASS 衛 星 訊 號 與 GPS 衛 星 訊 號 最 主 要 之 不 同 是 訊 號 調 變 方 式，

GLONASS 採用分頻多工(Frequency Division Multiple Access; FDMA)之訊號調變 方式，即不同 GLONASS 衛星傳遞之訊號頻率皆相異，各個衛星皆有專屬之頻道，

藉以分辨訊號來源，而 GPS 則採用分碼多工(Code Division Multiple Access; CDMA)，

即不同 GPS 衛星傳遞之訊號編碼皆相異，藉以分辨訊號來源。

2.2.1 載波頻率之分配

GLONASS 於 L1 與 L2 之載波頻率可由以下公式表示[2]：

1 01 1

fK = f +K f∆ (2.1)

2 02 2

fK = f +K f∆ (2.2) 其中 f01=1602MHz，∆f1=562.5kHz 為 L1 頻帶所使用；f02=1246MHz，∆f2=437.5kHz 為 L2 頻帶所使用。而 K 為衛星傳遞訊號之頻率號碼，K=-7,-6,…,0,+1,…,+6。

GLONASS 衛星載波頻率分佈如表 2.5 所示。

2.2.2 虛擬隨機測距碼

虛擬隨機測距碼(Pseudo Range Ranging Code)為一最大長度之序列，其係由移 位週期為 1 毫秒(Millisecond; ms)且位元率(Bit Rate)為 511kbps 之移位暫存器(Shift Register)所產生[2]。虛擬隨機測距碼取至 9 級移位暫存器中之第 7 級輸出，如圖 2.3。暫存器之初始值為 111111111，其週期為 1 毫秒。產生該序列之多項式為[2]：

5 9

( ) 1

G x = +x +x (2.3) 虛擬隨機距離碼與時脈產生之簡示方塊圖如圖 2.4 所示。

圖 2.3 移位暫存器之結構[2]

圖 2.4 虛擬隨機距離碼與時脈產生之簡示方塊圖[2]

12

2.2.3 導航訊息之產生

導航訊息之基礎單位為持續 2 秒不間斷之重複性資料串(String)，其前 1.7 秒為 85 位元之導航資料，後 0.3 秒為傳輸之時間標記。導航訊息係由以下兩二進位串 列之和再除以 2 取之餘數(Modulo 2 Addition)而得[2]：

導航訊息數位資料之相關碼所組成的位元序列，其 1 位元佔 20 毫秒

干擾序列(Meander Sequence)，其 1 位元佔 10 毫秒

時間標記之二進位碼為一 30 位元之虛擬隨機序列，其每一位元需 10 毫秒，

產生該序列之多項式為[2]：

3 5

( ) 1

G x = +x +x (2.4) 產生結果為 111110001101110101000010010110。資料序列產生之簡示方塊圖，如 圖 2.5。

圖 2.5 資料序列生成之簡示方塊圖[2]

表 2.5 GLONASS 衛星載波頻率分佈於 L1 與 L2[1]

衛星序號 軌道平面 頻率號碼 於 L1 之頻率(MHz) 於 L2 之頻率(MHz)

1 1 +1 1602.5625 1246.4375

2 1 -4 1599.75 1244.25

3 1 +5 1604.8125 1248.1875

4 1 +6 1605.375 1248.625

5 1 +1 1602.5625 1246.4375

6 1 -4 1599.75 1244.25

7 1 +5 1604.8125 1248.1875

8 1 +6 1605.375 1248.625

9 2 -2 1600.875 1245.125

10 2 -7 1598.0625 1242.9375

11 2 0 1602 1246

12 2 -1 1601.4375 1245.5625

13 2 -2 1600.875 1245.125

14 2 -7 1598.0625 1242.9375

15 2 0 1602 1246

16 2 -1 1601.4375 1245.5625

17 3 -2 1600.875 1245.125

18 3 -7 1598.0625 1242.9375

19 3 0 1602 1246

20 3 -1 1601.4375 1245.5625

21 3 +4 1604.25 1247.75

22 3 -3 1600.3125 1244.6875

23 3 +3 1603.6875 1247.3125

24 3 +2 1603.125 1246.875

14

每一數位資料串之第一個位元為 0，此為一閒置字元作為前一資料串時間標記 之增補字元。資料串串列中，包含資料位元、干擾位元、時間標記位元與測距碼 位元，各種類位元彼此同步於各傳輸導航訊號之範圍中。干擾位元與資料位元之 邊界同時位於測距碼初始位元之前緣，而時間標記最新位元之後緣對應到之時刻 當日時刻相異，存在整數秒之誤差。

導航訊息與虛擬隨機測距碼之二進位串列彼此脈衝同步關係，如圖 2.6。導航 訊息生成顯示於圖 2.7。導航訊息之內容與格式會詳述於下一節。

圖 2.6 導航訊息與虛擬隨機測距碼之二進位串列脈衝同步關係圖[2]

圖 2.7 導航訊息生成圖[2]

2.2.4 導航訊息內容

導航訊息主要目的為提供使用者必要之數據，以達定位、時間校準及導航等 功能。導航訊息之內容主要分為兩大部分，即時性資料與非即時性資料。即時性 資料屬各獨立衛星之資訊，包含列舉之衛星時間標記、衛星時間與 GLONASS 系 統時間之差值、衛星載波頻率與其標稱值之差值以及其他軌道參數等廣播星曆資 料(Ephemeris Data)；非即時性資料屬衛星群體資訊，包含衛星之粗略狀態(Status Almanac)、衛星時間與 GLONASS 系統時間之粗略修正量(Phase Almanac)、衛星粗 略之軌道參數(Orbit Almanac)以及 GLONASS 系統時間與 UTC(SU)時間標準之修 正量等粗略星曆資料(Almanac Data)[2]。

廣播星曆(Ephemeris)及粗略星曆(Almanac)簡介如下[2]：

(1)廣播星曆(Ephemeris)：

系統控制中心週期性地計算及上載廣播星曆參數，GLONASS 衛星座標值與速 度之方均根誤差列於表 2.6 中。

16

表 2.6 GLONASS 衛星座標值與速度精度[2]

誤差分量 方均根誤差值

Predicted coordinates (m) Velocity (cm/s) 衛星種類 GLONASS GLONASS-M GLONASS GLONASS-M

Along track 20 7 0.05 0.03

Cross track 10 7 0.1 0.03

Radial 5 1.5 0.3 0.2

廣播星曆資料各字元之特性列於表 2.8 中，舉其中幾項重要之星曆資料加以解 釋。

tb：當日根據 UTC(SU)+03hours00min 之時間區間指標。

P1：即時性資料上載之旗標，其顯示兩相鄰時間指標之時間間隔，內容見表 2.7。

n：衛星之指標，其與衛星之序號相符合。

FT：偽距(pseudorange)測量值之精度，其內容見表 2.9。

γn(tb)：載波頻率預測值與標稱值之相對偏差，定義如下：

( ) ⁿ( )^{b Nn}

n b

Nn

f t f

t f

γ = ⁻ (2.5)

其中 fn(tb)為第 n 顆衛星之載波預測值，fNn為第 n 顆衛星之載波標稱值。

τn(tb)：第 n 顆衛星時間 tn與 GLONASS 系統時間 tc之修正量，其計算式如下：

( ) ( ) ( )

n tb t tc b t tn b

τ = − (2.6)

表 2.7 P1 旗標內容[2]

字元 P1 兩相鄰時間指標 tb之時間間隔 (min)

00 0

01 30

10 45

11 60

表 2.8 廣播星曆參數之字元特性[2]

字元 位元數 比例因子(LSB) 有效範圍 單位

M 4 1 0…15 無單位

tk

5 1 0…23 時

6 1 0…59 分

1 30 0;30 秒

tb 7 15 15…1425 分

M ⁽¹⁾ 2 1 0;1 無單位

γn(tb) ⁽²⁾ 11 2^-40 ±2^-30 無單位 τn(tb) ⁽²⁾ 22 2^-30 ±2^-9 秒 xn(tb), yn(tb), zn(tb) ⁽²⁾ 27 2^-11 ±2.7·10⁴ 公里 xn(tb), yn(tb), zn(tb) ⁽²⁾ 24 2^-20 ±4.3 公里/秒 xn(tb), yn(tb), zn(tb) ⁽²⁾ 5 2^-30 ±6.2·10^-9 公里/秒²

B_n 3 1 0…7 無單位

P ⁽¹⁾ 2 1 00,01,10,11 無單位

NT (1) 11 1 0…1461 天

FT (1)

4 見表 2.9

n ⁽¹⁾ 5 1 0…31 無單位

∆τn (2)

5 2^-30 ±13.97·10^-9 秒

En 5 1 0…31 天

P1 2 見表 2.7

P2 1 1 0;1 無單位

P3 1 1 0;1 無單位

P4 ⁽¹⁾ 1 1 0;1 無單位

ln

(1) 1 1 0;1 無單位

註：標記(1)：此類字元為 GLONASS-M 中之導航訊息。

標記(2)：此類字元有正負符號之分，其第一位元表示正負號，0 表示正號 1 表示負號。

18

表 2.9 FT旗標內容[2]

字元 FT 測量精度 (m)

0 1

1 2

2 2.5

3 4

4 5

5 7

6 10

7 12

8 14

9 16

10 32

11 64

12 128

13 256

14 512

15 不使用

(2) 粗略星曆(Almanac)：

粗略星曆即為非即時性資訊，其中包含 GLONASS 時間上之資料、各個 GLONASS 衛星之機上時間尺度的資訊以及各個 GLONASS 衛星之軌道參數與健 康狀態等，其資訊皆與整體 GLONASS 衛星星群相關。由於解定位時以廣播星曆 為主，故於此對粗略星曆資料不多作介紹，僅舉相關之時間校正參數。

τc：GLONASS 時間與 UTC(SU)時間之時間尺度校正量。

τGPS：GLONASS 時間與 GPS 時間之校正量，關係式如下：

GPS GLONASS GPS

T −T = ∆ +T τ (2.7) 其中 ∆T 為整數部分，來自於 GPS 導航訊息，而 τGPS為小數部分。

2.2.5 導航訊息架構

導航訊息是以數位資料之型式傳遞，其透過漢明碼(Hamming Code)編碼而成。

結構上，這些資料為一串連續且重複之超訊框(Super Frame)，超訊框內含多個訊框 (Frame)，而訊框內含多個資料串(String)。不同衛星之資料串、訊框及超訊框彼此 同步於 2 毫秒內[2]。

超訊框完整傳送需持續 2.5 分鐘，其內含 5 個訊框，任一訊框之傳輸時間為 30 秒，其內含 15 資料串，任一資料串之傳輸時間為 2 秒。不同衛星發送之完整超 訊框內，含整體衛星星群之非即時性資料，卻僅含該衛星之即時性資料，超訊框 完整架構如圖 2.8。每一訊框內含該衛星之完整即時性資料，卻僅含部分衛星之非 即時性資料，各訊框內含非即時性資料之衛星序號見表 2.10。

每一衛星發送之超訊框內，第一訊框至第四訊框結構相同，各訊框之第一至 第四資料串為該衛星之即時性資料，而第六至第十五資料串為該訊框對應之衛星 序號的非即時性資料。每一衛星之非即時性資料佔據 2 個資料串，第一至第四訊 框皆含 5 顆衛星之非即時性資料，其非即時性資料共佔 10 個資料串，而第五訊框 僅含 4 顆衛星之非即時性資料，故其結構與前四訊框相異。訊框架構詳見圖 2.9(a)(b)。

表 2.10 GLONASS 非即時性資料於超訊框內之分配[2]

Frame number within superframe Satellite numbers, for which almanac is transmitted within given superframe

1 1-5

2 6-10

3 11-15

4 16-20

5 21-24

20

圖 2.8 超訊框架構[2]

圖 2.9(a) 第一至第四訊框架構[2]

圖 2.9(b) 第五訊框架構[2]

22

每一資料串之傳遞皆需持續 2 秒，後 1.3 秒為時間標記， 其共含 30 位元，每 一位元占 10 毫秒，前 1.7 秒為 85 位元之導航資料。資料串之結構詳見圖 2.10。

圖 2.10 資料串架構[2]

2.3 座標系統與時間系統座標系統與時間系統座標系統與時間系統座標系統與時間系統

在天文學與太空科學領域中，時間與空間之描述成為研究載具運動之重要部 分。座標系統是精確描述天體與人造衛星相互關係之重要基準，而時間系統同時 也是進行導航及定位時之重要元素，故瞭解空間系統與時間系統之基本概念，有 助於導航衛星系統之應用。

GLONASS 與 GPS 兩者所使用之座標系統與時間系統皆不同，GLONASS 之 座標系統為 PZ-90.02 而 GPS 則為 WGS84 座標系統；GLONASS 之時間系統為 UTC(SU)，而 GPS 為 UTC(USNO)時間系統。本節將介紹 PZ-90.02 座標系統與 UTC(SU)時間系統，以及整合兩系統時需特別注意之系統轉換方式。

2.3.1 座標系統類別

GLONASS 衛星量測中，主要之參考座標系統有兩種：地心慣性( Earth-Centered Inertial; ECI)座標系統與地心地固(Earth-Centered Earth Fixed; ECEF)座標系統[5]。

(1)地心慣性座標系統：

地心慣性座標系統為一慣性座標系，不隨地球自轉而轉動，但會隨地球公轉 而移動，即其三軸方向不變但其原點會隨地球位置而移動。以公元 2000 年 1 月 1 日正午 12 時為該座標系統之基準(J2000)，其原點為地心，X 軸指向 J2000 之平均 春分點，Z 軸指向平均北極點，Y 軸與 X 軸及 Z 軸構成一右旋直角座標系統。春 分點係指春分時太陽直射地球赤道，此時赤道面與黃道面於地球上之交點。

(2)地心地固座標系統：

地心地固座標系統為一隨地球自轉而旋轉之座標系統，其原點位於地心，X 軸指向本初子午線與赤道交點，Z 軸指向平均北極，Y 軸與 X 軸及 Z 軸構成一右 旋直角座標系統。

地心慣性座標系統與地心地固座標系統關係如圖 2.11。

圖 2.11 地心慣性座標系統與地心地固座標系統關係示意圖[6]

24

2.3.2 PZ-90.02 座標系統

GLONASS 傳遞之廣播星曆資料，其中位置資訊是採用 PZ-90.02 座標系統所 描述，該座標系統為一地心地固座標系統。PZ-90.02 座標原點位於地心，其 Z 軸 指向由國際地球自轉服務機構(International Earth Rotation Service; IERS)所示之地 球北極，X 軸指向本初子午線與地球赤道之交點，Y 軸與 X 軸及 Z 軸構成一右旋 直角座標系統[2]。基於 PZ-90.02 座標系統定義之地球常數與相關重要參數，皆列 舉於表 2.11 中。

表 2.11 PZ-90.02 座標系統之常數與相關參數[2]

項 目 參 數 值 地球自轉率 7.292115×10^-5 radians/s 重力常數 398600.4418×10⁹ m³/s² 大氣重力常數 (fMa) 0.35×10⁹ m³/s²

光速 299792458 m/s

地球半長軸 6378136 m

地球扁率 1/298.25784

地球赤道之重力加速度 978032.84 mGal 海平面重力加速度之大氣修正量 0.87 mGal 第二階帶諧係數 (J20

) 1082625.75×10^-9 第四階帶諧係數 (J40

) -2370.89×10^-9 橢球表面之正向位能 (U0) 62636861.4 m²/s²

GLONASS 採用 PZ-90.02 座標系統之定義與 GPS 之 WGS84 不盡相同，整合 兩系統之座標系統必須經過座標轉換，將 PZ-90.02 座標系統下之座標轉換至

WGS84 座標系統，轉換公式如下[7]：

( )

84 90.02

1

1 1

WGS 1 PZ

x x u

y y s v

z z w

δω δψ

δ δω δε

δψ δε ₋

∆ −

       

       

= ∆ + + ⋅ − ⋅

       

∆ −

       

       

(2.8)

其 中 各 參 數 之 值 有 眾 多 相 關 研 究 ， 國 際 GLONASS 實 驗 組 織 (International GLONASS Experiment 1998; IGEX-98)經長期觀測與分析而提出以下參數值[8]：

參數 ∆x ∆y ∆z δε δψ δω δs

值 0.06m 0.07m -0.57m 1.697·10^-7 -1.018·10^-7 -1.736·10^-6 -1·10^-8

2.3.3 UTC(SU)時間系統

GLONASS 之時間系統是以 UTC(SU)時間系統為參考基準，該系統時間與國 際參考時間 UTC(SU)相差約 1 毫秒，透過導航訊息中之參數修正，可將其時間差 調整至 1 微秒內，每日會有兩次之時間修正。GLONASS 系統時間是依據 GLONASS 中央同步裝置之時間(Central Synchronizer time)而得，其每日不穩定性優於 2×10^-15。 GLONASS 時間係參照國際時間局(Bureau International de l'Heure; BIH)公告之潤秒 修正而跟進修正，每年或每一年半作一次潤秒修正[2]。

GLONASS 衛星會持續不斷地將時間刻度調整至 GLONASS 系統之標準時間，

使各時刻 GLONASS 衛星時間與 GLONASS 系統標準時間差小於 10 奈秒。然而，

該系統時間與 UTC(SU)固定存在 3 小時之差異，其關係如下[2]：

( ) 03hours00min

GLONASS UTC SU

t =t + (2.9) 而 GLONASS 衛星時間與 UTC(SU)之關係可利用廣播星曆資料修正而得，其關係 如下[2]：

( )

( ) 03hours00min ( ) ( )

UTC SU c n b n b b

t + = +t τ +τ t −γ t ⋅ t−t (2.10) 其中 τc、τ_n、tb、γ_n為星曆參數。

26

第 第 第

第三 三 三章 三 章 章 章 GLONASS 衛星軌道模型 衛星軌道模型 衛星軌道模型 衛星軌道模型

在整個定位過程中，吾人將其分為兩大階段，第一階段係由廣播星曆計算出 各個衛星於各個時刻之精確位置，得到各時刻各個衛星之精確位置後，即進入第 二階段利用相關定位方法求得各個時刻之接收機位置，相關定位方法將於第四章 詳細介紹。

本章將先建立 GLONASS 衛星軌道模型，推得衛星運動之微分方程式，再透 過數值積分龍格-庫塔法，積分得到衛星之速度與位置，最後將其結果與國際 GNSS 服務網之衛星位置及 Novatel 接收機算出之衛星位置比較，求得誤差並分析之。

3.1 衛星軌道方程式之推導衛星軌道方程式之推導 衛星軌道方程式之推導衛星軌道方程式之推導

GLONASS 衛星位置之計算與 GPS 相當不同，GPS 衛星位置主要是依據軌道

六元素計算而得，而 GLONASS 廣播星曆則提供該系統衛星於 PZ-90.02 之位置、

速度與太陽與月球引力對衛星造成之加速度，該系統之廣播星曆一般約 30 分鐘上 載更新一次，故於 30 分鐘內 GLONASS 之衛星位置與速度，可由廣播星曆提供的 數據當作初始值，對衛星運動之微分方程式進行積分求得。

根據牛頓運動定律，衛星繞地球運行是由於地球與衛星之間有作用力存在，

而主要之作用力為地球重力場造成，其重力位能之關係可以表示為[9]

V U

r

= µ+ (3.1)

其中 V 為總重力位能 µ 為地球重力常數 r 為衛星與地心之距離

U 為地球重力位能非球體部分

28

若只考慮地球重力位能球體部分，且忽略其他非重力影響，衛星之運動即為一純 克卜勒運動。事實上，並不能忽略其他力之影響，因此不能用一般之克卜勒橢圓 來描述衛星運動。GPS 為了補償此影響，該系統於廣播星曆中加入了一些修正參 數用於一假想之克卜勒軌道，而 GLONASS 則是利用不同方法解決此問題。

根據 Heiskanen 與 Moritz 於 1967 年提出之展開法，地球重力位能非球體部分 U 可展開成球諧函數(Spherical Harmonics)，其表示式為[10]

( ) ( )

2 0

cos cos sin

n n E

nm nm nm

n m

U a P c m s m

r r

µ θ λ λ

∞

= =

 

=   ⋅ +

 

∑∑

其中 aE 為地球赤道半徑

r, λ, θ, 為地球固定極座標(半徑、經度、餘緯度) n, m 為球諧函數展開之階層(Degree and Order) Pnm(cosθ) 為相關雷建德函數(Legendre Functions) cnm, snm 為球諧函數之係數

由於地球重力位能於一階近似為旋轉對稱(即與 λ 獨立)，其帶諧部分(Zonal Harmonics, 即 (3.2) 式 中 m=0) 對 衛 星 運 動 之 影 響 ， 遠 比 其 田 諧 部 分 (Tesseral Harmonics, 即(3.2)式中 0<m<n)與扇諧部分(Sectorial Harmonics, 即(3.2)式中 m=n) 重要。因此，於短時間軌道積分時，可以忽略田諧與扇諧部分之影響，而(3.2)式可 簡化為

( )

1

0 0

2

cos

n E

n n

E n

U a c P

a r

µ θ

∞ +

=

 

=   ⋅

 

∑

1996 年，Spilker 整理出各擾動力對 GPS 衛星之影響，其大約影響列於表 3.1 中。由表 3.1 中可發現，僅太陽與月球重力之影響比二階帶諧少一個級數(Order)，

其餘影響皆比二階帶諧少兩個級數以上，因此可忽略其其餘擾動影響。由於各系

統衛星所受之力應大致相同，固可將此結果應用於 GLONASS 衛星之上，即(3.3) 式僅保留二次項，可將(3.3)式表示為

( )

2

20 20

3^E cos

U a c P

r

µ θ

= ⋅ (3.4)

因雷建德多項式 P20(cosθ) = 3/2 cos²θ – 1/2，可將(3.4)式改寫為

2

2 3 20

3 1

2cos 2 aE

U c

r

µ ^ θ ^

= ⋅ − 

  (3.5) 總重力位能即可表示為

( )

2

2 3 20

1 3cos 1

2 aE

V c

r r

µ µ

θ

= + ⋅ − (3.6)

表 3.1 各擾動力對 GPS 衛星之影響[9]

擾動力來源 最大擾動加速度 (m/s²) 一小時內之最大偏移量 (m)

Spherical Earth 5.65·10^-1 -

Second zonal harmonic 5.3·10^-5 300

Fourth zonal harmonic 10^-7 0.6

Gravity anomalies 10^-8 0.06

Lunar gravity 5.5·10^-6 40

Solar gravity 3·10^-6 20

Solar radiation pressure 10^-7 0.6

All other forces 10^-8 0.06

重力位能造成之衛星加速度於卡式座定義為

xr= ∇V

&& (3.7)

30 其分量型式可寫為

i

i i i i

dV V r V V

x dx r x x x

λ θ

λ θ

∂ ∂ ∂ ∂ ∂ ∂

= = + +

∂ ∂ ∂ ∂ ∂ ∂

&&r

(3.8)

其中 xi為 x, y, z。

上式各項偏微分結果為

^{2 3 42 20}

(

)

2

V aE

r r r c

µ µ

∂ θ

= − − ⋅ −

∂

V 0 λ

∂ =

∂

2 3 20

3 ^E cos sin

V a

r c

µ θ θ

θ

∂ = − ⋅

∂

且

r x x r

∂ =

∂ , r y y r

∂ =

∂ , r z z r

∂ =

∂

1 ₃ sin

xz

x r

θ

θ

∂ = ⋅

∂ ^{, 3}

1 sin

yz

y r

θ

θ

∂ = ⋅

∂ ^,

2 2

3

1 sin

r z

z r

θ

θ

∂ −

= − ⋅

∂

將上述結果代入(3.8)式，可得衛星加速度為

2 2

3 20 5 2

3 1 5

2

aE z

x x c x

r r r

µ µ ^{ }

= − ⋅ + ⋅ ⋅ − 

 

&& (3.9)

2 2

3 20 5 2

3 1 5

2

aE z

y y c y

r r r

µ µ ^{ }

= − ⋅ + ⋅ ⋅ − 

 

&& ( 3 . 8 )

2 2

3 20 5 2

3 3 5

2

aE z

z z c z

r r r

µ µ ^{ }

= − ⋅ + ⋅ ⋅ − 

 

&& ( 3 . 8 )

太陽與月球引力對衛星造成之加速度，於短時間積分內可視為常數&&xLS_, &&y_LS, z&&LS ，且忽略先前所述之其他擾動力，則完整 GLONASS 衛星之加速度可表示為

2 2

3 20 5 2

3 1 5

2

E

LS

a z

x x c x x

r r r

µ µ ^{ }

= − ⋅ + ⋅ ⋅ − +

 

&& && (3.10)

2 2

3 20 5 2

3 1 5

2

E

LS

a z

y y c y y

r r r

µ µ ^{ }

= − ⋅ + ⋅ ⋅ − +

 

&& && ( 3 . 1 8 )

2 2

3 20 5 2

3 3 5

2

E

LS

a z

z z c z z

r r r

µ µ ^{ }

= − ⋅ + ⋅ ⋅ − +

 

&& && ( 3 1 . 8 )

由於牛頓運動定律只滿足於慣性座標系統，上述推導過程皆於慣性座標系統，

故(3.10)式只適用於慣性座標系統。而 GLONASS 之參考座標系統 PZ-90.02 為一地 心地固(ECEF)座標系統，且定位時，衛星運動描述於地心地固座標系統，對使用 上較為方便且直接，因此需要透過一系列之座標轉換，將衛星運動方程式(3.10)轉 換至地心地固座標系統。

慣性座標系統轉換至地心地固座標系統，包含地球自轉、進動(Precession)、

章動(Nutation)及極動(Polar motion)等影響，其轉換式如下[6]

M R N P

ECEF INS

xr = R ⋅R ⋅R ⋅R ⋅xr

(3.11) 其中 xrECEF

為地心地固座標系統之座標 xrINS

為地心慣性座標系統之座標 R ^M 為極動之旋轉矩陣

R ^R 為地球自轉之旋轉矩陣 R ^N 為章動之旋轉矩陣 R^P 為進動之旋轉矩陣

32

地球之極動、章動及進動是非常緩慢的過程，故在短時間積分內其對座標轉 換之影響不大，而地球自轉為主要影響。地球自轉導致柯氏力之影響，衛星運動 方程式最終可表示為[2][9]

dx x

dt = & (3.12)

dy y

dt = & ( 3 . 1 8 )

dz z

dt = & ( 3 1 . 8 )

2 2

2

3 20 5 2

3 1 5 2

2

E

LS E E

dx a z

x c x x x y

dt r r r

µ µ

ω ω

 

= − ⋅ + ⋅ ⋅ − + + ⋅ + ⋅

 

&

&& &_{0 )}

2 2

2

3 20 5 2

3 1 5 2

2

E

LS E E

dy a z

y c y y y x

dt r r r

µ µ

ω ω

 

= − ⋅ + ⋅ ⋅ − + + ⋅ − ⋅

 

&

&& & _{8 )}

2 2

3 20 5 2

3 3 5

2

E

LS

dz a z

z c z z

dt r r r

µ µ ^{ }

= − ⋅ + ⋅ ⋅ − +

 

&

&& ( 3 1 . 8 )

其中 x,y,z 為衛星座標位置

x&, y& , z& 為衛星速度 xLS

&& ,&&y_LS,z&& LS 為太陽與月球引力對衛星造成之加速度 r =^+ ^+ ^ 為衛星與地心距離

µ = 398600.4418×10⁹ m³/s² 為地球重力常數

c₂₀ = −1082625.75×10^-9 為二階諧帶係數(由 J20得) a_E = 6378136 m 為地球赤道半徑

ωE = 7.292115×10^-5 radians/s 為地球自轉率

(3.12)式為 GLONASS 衛星運動方程式，其描述之衛星座標於地心地固座標系 統(ECEF)，積分該式即可得到衛星於地心地固座標系統之位置。GLONASS 之廣 播星曆資料中包含 tb參考時刻之衛星位置 x ,y,z、速度x&, y& , z& 及太陽與月球引力 對衛星所造成之加速度&& ,xLS &&yLS,z&& ，這些數據可為(3.12)微分方程式中之初始值。LS

而 µ, c20,aE,ωE則採用 PZ-90.02 座標系統之常數與相關參數，見表 2.11。

3.2 數值積分法數值積分法 數值積分法數值積分法

前一節中介紹了衛星運動方程式，若要得知衛星運動情形，則必需對(3.12)

式作積分以求出衛星各個時刻之位置與速度。本論文使用四階龍格-庫塔數值積分 法求解衛星之位置與速度，以下將詳細介紹該方法[6][11]。

GLONASS 廣播星曆資料提供之衛星位置、速度及太陽與月球引力對衛星所造 成之加速度可當作初始值，如此求解各時刻衛星位置與速度即可視為初始值問題 (Initial Value Problem)，而透過尤拉方法(Euler’s Method)可簡單解決此類問題。

由微積分之定義，y 對 x 之導數可表示為(本節之 x 及 y 與前一節無關)

1 1

( i ) ( )i

i i

dy y x y x

y dx x x

+ +

′ = ≈ −

− (3.13) 若xi+1−xi =h，h 為步階，可將y′視為(x,y)之函數，其初始值(x0,y₀)代入上式得

1 0

0 0 0 0

( ) ( ) ( , ) ( , ( )) y x y x y f x y f x y x

h

′ = = ≈ − (3.14)

經移項整理後可得

1 0 0 0

( ) ( ) ( , ( ))

y x ≈ y x +hf x y x (3.15) 從上式得知，利用(x y x0, ( 0))可求得y x ，且( )1 x1=x0+h，故可類推於y x( 2)之後的 計算，即

( i 1) ( )i ( , ( ))i i

y x₊ ≈ y x +hf x y x , i = 1,2,3,… (3.16)

34 (3.16)式稱之為尤拉方法。

承接上式，於x0 ≤xi ≤ xn時，若對(x y₀, ₀= y x( ₀))作泰勒展開可得

2

( 1)

1 ( , ) '( , ) ... ( , )

2! !

m m

i i i i i i i i

h h

y y hf x y f x y f x y

m

−

+ = + + + + (3.17)

m 為展開之階數，(3.17)式即為 m 階泰勒法(mth-order Taylor Method)，因此尤拉方 法相當於一階泰勒法。

相較於尤拉方法是以單一導數 ( ,f x y 來計算，龍格-庫塔法之原理則是取多i i) 筆導數函式加權平均，其方法表示如下

1

1 n

i i j j

j

y+ y h ω k

=

= +

∑

二階龍格-庫塔法表示如下

1 ( 1 1 2 2)

i i

y+ = y +h ωk +ω k (3.19) 其中

1

2 1

( , )

( , )

i i

i i

k f x y

k f x ah y bk

 =

 = + +

 (3.20) 而ω ω₁, ₂, ,a b為待定參數，且

( ,_{i i}) dy f x y y

′ = dx (3.21) 再來推算 y 之二階導數

( , )_{i i} ( )_i

y′ = f x y = f (3.22)

( ,_{i i}) ( _{x y} )_i

y′′= f x y′ = f + f f (3.22) 將k 泰勒展開至一階可得 ₂