ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZE FAKULTA STAVEBNÍ – OBOR GEODÉZIE A KARTOGRAFIE KATEDRA VYŠŠÍ GEODÉZIE název předmětu
Vyšší geodézie 1 úloha/zadání
název úlohy
GPS - Výpočet drah družic
2/3 školní rok
2010/11
semestr
1
skupina
NG1-88
zpracoval
datum
Jan Dolista
29. 10.
klasifikace
GPS - výpočet drah družic Zadání: K dispozici máte soubor vysílaných efemerid ve formátu RINEX (tzv. navigační zprávu, označeni koncovkou ’n’) pro stanici sítě EUREF EPN vytvořený dne 17.9.2009 a kartézské souřadnice stanice v systému WGS84. Pro družici danou PRN kódem spočítejte pomocí vysílaných efemerid její polohu v čase 𝑡1 dne 17.9.2009 (den 260). Souřadnice družice určete v systému WGS84, ve kterém jsou zadány i dráhové elementy. Pro výpočet Keplerových elementů v čase 𝑡1 použijte vysílané efemeridy v epoše nejbližší nižší zadanému času. Dále proveďte transformaci souřadnic družice do topocentrického pravotočivého systému (N W U). Jako zeměpisné souřadnice pozorovatele použijte zadané souřadnice stanice. Kontrolně vypočtěte délku družice-stanice v obou souřadnicových systémech. Pro výpočet Keplerových elementů v čase 𝑡1 použijte vysílané efemeridy v epoše nejbližší nižší zadanému času. Výsledkem tedy bude: ∙ Keplerovy elementy v čase 𝑡1 ∙ Souřadnice družice v systému WGS84 ∙ Souřadnice družice v topocentrickém systému (N W U) ∙ Délka stanice-družice vypočtená v obou souřadnicových systémech Jako geocentrickou gravitační konstantu použijte hodnotu 𝐺𝑀 = 398, 6005 · 1012 𝑚3 𝑠−2 . Jako úhlovou rychlost rotace Země 𝜔𝐸 = 7292115, 1467 · 10−11 𝑟𝑎𝑑 · 𝑠−1 .
Číselné zadání 3: číslo stanice EPN zadání 3 GOPE
PRN 12
hod 16
čas 𝑡1 min sec 12 45
X 3979316
souřadnice Y Z 1050312 4857067
Vypracování: Veškeré výpočty byly provedeny v programu Octave. Veškeré úhly jsou uváděny v radiánech a nejsou přičítány násobky 2𝜋, délky a souřadnice jsou uváděny v metrech.
1
Načtení parametrů z navigační zprávy
Ze zadaného navigačního souboru RINEX byly pomocí skriptu v programu Octave vybrány záznamy pro zadané PRN družice. Soubor navigační zprávy obsahuje nejprve hlavičku, za níž následují jednotlivé záznamy ve tvaru: PRN — 𝐶𝑢𝑐 𝑇 𝑜𝑒 𝐼0 𝐼˙ .. .
ROK 𝐶𝑟𝑠 𝑒 𝐶𝑖𝑐 𝐶𝑟𝑐 —
MĚSÍC ∆𝑛 𝐶𝑢𝑐 𝑙0 𝜔0 —
DEN 𝑀0 √ 𝑎 𝐶𝑖𝑠 ˙ Ω
HODINA
MINUTA
SEKUNDA
𝛿𝑐
𝜕𝛿𝑐 𝜕𝑡
𝜕 2 𝛿𝑐 𝜕𝑡2
—
Skriptem je tedy nejprve detekován konec hlavičky a poté vyhledáván vždy první řádek každého záznamu. Z tohoto řádku je načteno PRN a vyhodnoceno zda záznam přísluší zadané družici. Po vybrání všech záznamů pro dané PRN je vyhledán záznam s nejbližším nižším časem než je zadaný čas 𝑡1 . Při vyhledávání času je zohledněno i datum, neboť navigační zpráva může obahovat i záznamy z předchozího nebo následujícího dne.
Z vybraného záznamu (nejbližší nižší čas pro dané PRN) jsou do proměnných načteny oskulační elementy a další korekční členy.
2
Výpočet Keplerovských oskulačních elementů
Postupně byly vypočteny Keplerovské oskulační elementy
2.1
Střední anomálie 𝑀 = 𝑀0 + 𝑛(𝑡1 − 𝑡0 ) + ∆𝑛(𝑡1 − 𝑡0 ) = −1.520259187026𝑟𝑎𝑑
2.2
Excentrická anomálie 𝐸 = 𝑀 + 𝑒 sin 𝐸 = −1.523349928908𝑟𝑎𝑑
Excentrická anomálie byla určena iteračně, kdy v první iteraci je hodnota excentrické anomálie volena 𝐸0 = 𝑀 + 𝑒 sin 𝑀 . V dalších iteracích je hodnota anomálie dána vztahem 𝐸𝑖 = 𝑀 + 𝑒 sin 𝐸𝑖−1 . Výpočet je opakován, dokud rozdíl ve dvou po sobě jdoucích iteracích není menší než dvojnásobek přesnosti výpočtu (dáno možnostmi počítače).
2.3
Pravá anomálie
(︃ √
1+𝑒 𝐸 𝑣 = 2𝑎𝑡𝑎𝑛 √ 𝑡𝑔 1−𝑒 2
2.4
)︃
= −1.526440902535𝑟𝑎𝑑
Argument šířky 𝑢0 = 𝜔0 + 𝑣 = −2.113224826955𝑟𝑎𝑑 𝜔 = 𝜔0 + 𝐶𝑢𝑐 cos (2𝑢0 ) + 𝐶𝑢𝑠 sin (2𝑢0 ) = −0.586779172135𝑟𝑎𝑑 𝑢 = 𝜔 + 𝑣 = −2.113220074670𝑟𝑎𝑑
2.5
Průvodič 𝑟0 = 𝑎(1 − 𝑒 cos 𝐸) = 26555411.2𝑚 𝑟 = 𝑟0 + 𝐶𝑟𝑐 cos (2𝑢) + 𝐶𝑟𝑠 sin (2𝑢) = 26555207.4𝑚
Po výpočtu argumentu šířky a délky průvodiče byly určeny souřadnice družice v systému, kde oběžná dráha družice je v rovině xy, a osa x směřuje do výstupního uzlu (průsečík dráhy s rovinou rovníku). ⎛ ⎞ ⎛ ⎞ cos 𝑢 −13708152.5 ⎜ ⎟ ⎜ ⎟ XS = 𝑟 ⎝ sin 𝑢 ⎠ = ⎝ −22743473.6 ⎠ 0 0
2.6
Sklon ˙ 1 − 𝑡0 ) = 0.967734197596𝑟𝑎𝑑 𝐼 = 𝐼0 + 𝐶𝑖𝑐 cos (2𝑢) + 𝐶𝑖𝑠 sin (2𝑢) + 𝐼(𝑡
2.7
Délka výstupního uzlu ˙ − 𝜔𝐸 )(𝑡1 − 𝑡0 ) − 𝑇 𝑜𝑒𝜔𝐸 = −27.665143827554𝑟𝑎𝑑 𝑙 = 𝑙0 + (Ω
3
Rotace do WGS
Rotace do systému WGS probíhá ve dvou krocích. Nejprve sklopením kolem osy x o úhel 𝐼 v matematicky záporném smyslu. Tím je ztotožněna osa z s osou Z systému WGS. Druhým krokem je pak rotace kolem osy Z o úhel 𝑙 v matematicky záporném smyslu. XWGS = RZ (−𝑙)RX (−𝐼)XS ⎞
⎛
XWGS
3861164.8 ⎟ ⎜ = ⎝ 18422760.0 ⎠ −18731587.2 ⎞
⎛
4
⎞
⎛
1 0 0 ⎟ ⎜ RX (−𝐼) = ⎝ 0 cos (−𝐼) sin (−𝐼) ⎠ 0 − sin (−𝐼) cos (−𝐼)
cos (−𝑙) sin (−𝑙) 0 ⎟ ⎜ RZ (−𝑙) = ⎝ − sin (−𝑙) cos (−𝑙) 0 ⎠ 0 0 1
Transformace do lokálního (topocentrického) systému (NWU)
Rotace do systému NWU probíhá přes pomocný systém TOPO do kterého je transformace provedena ve třech krocích. Nejprve je v systému WGS provedem posun počátku ze středu elipsoidu do stanoviska stanice. Ve druhém kroku je provedena rotace kolem osy Z o úhel 𝜆 (geodetická délka stanoviska stanice) v matematicky kladném smyslu. Posledním krokem je pak rotace kolem nové osy Y o úhel 𝐵 (geodetická šířka stanoviska stanice). stanice XTOPO = RY (−𝐵)RZ (𝜆)(Xdruzice WGS − XWGS )
⎛
XTOPO ⎛
⎞
−15265751.8 ⎜ ⎟ = ⎝ 16827355.1 ⎠ −18494254.8 ⎞
⎛
cos (−𝐵) 0 − sin (−𝐵) ⎜ ⎟ 0 1 0 RY (−𝐵) = ⎝ ⎠ sin (−𝐵) 0 cos (−𝐵)
4.1
⎞
cos (𝜆) sin (𝜆) 0 ⎜ ⎟ RZ (𝜆) = ⎝ − sin (𝜆) cos (𝜆) 0 ⎠ 0 0 1
Výpočet geodetické šířky 𝐵 a geodetické délky 𝜆
K realizaci rotací bylo nutné určit zeměpisnou délku a geodetickou šířku stanoviska stanice. Ty lze přibližně, avšak dostatečně přesně, určit ze vztahů: (︂ )︂ 𝑍 ·𝑎 Θ = 𝑎𝑡𝑎𝑛 √ = 0.869503762462𝑟𝑎𝑑, 𝑋2 + 𝑌 2 · 𝑏 kde 𝑋, 𝑌, 𝑍 jsou souřadnice stanice v systému WGS a 𝑎 = 6378137𝑚, 𝑏 = 6356752.3142𝑚 jsou hlavní a vedlejší poloasa elipsoidu WGS-84. Geodetická šířka: (︃
𝑍 + 𝑒′2 𝑏 sin3 Θ 𝐵 = 𝑎𝑡𝑎𝑛 √ 𝑋 2 + 𝑌 2 − 𝑒2 𝑎 cos3 Θ kde 𝑒2 =
𝑎2 −𝑏2 𝑎2
a 𝑒′2 =
𝑎2 −𝑏2 𝑏2
)︃
= 0.871158507502𝑟𝑎𝑑,
jsou první a druhá excentricita elipsoidu WGS-84.
Geodetická délka:
𝑋 = 0.258057687805𝑟𝑎𝑑 𝑌 V pomocném systému TOPO směřuje kladná osa 𝑋 vzhůru, kladná osa 𝑌 k východu a kladná osa 𝑍 k severu. Transformace do systému NWU je pouze záměnou souřadnic. 𝜆 = 𝑎𝑡𝑎𝑛
⎛
XNWU
⎞
⎛
⎞
XTOPO (3) −18494254.8 ⎜ ⎟ ⎜ ⎟ = ⎝ −XTOPO (2) ⎠ = ⎝ −16827355.1 ⎠ XTOPO (1) −15265751.8
5
Výpočet délky družice - stanice v obou systémech
Na závěr byla vypočtena délka mezi družicí a stanicí v obou souřadnicových systémech: stanice 𝑟1 = |Xdruzice WGS − XWGS | = 29295742.3𝑚
𝑟2 = |XNWU | = 29295742.3𝑚 Porovnání obou délek bylo použito jako kontrola provedených rotací.
6
Číselné výsledky
Souřadnice družice v systému WGS:
XWGS
⎛
⎞
⎛
⎞
3861164.8 ⎟ ⎜ = ⎝ 18422760.0 ⎠ −18731587.2
Souřadnice družice v systému NWU:
XNWU
−18494254.8 ⎜ ⎟ = ⎝ −16827355.1 ⎠ −15265751.8
Délka stanice - družice: 𝑟 = 29295742.3𝑚
Závěr: Pro zadané PRN družice a čas 𝑡1 byly na základě navigační zprávy určeny souřadnice družice v systémech WGS-84 a místním topocentrickém systému NWU pro danou stanici. Ze souřadnic v obou systémech byla vypočtena délka mezi stanicí a družicí. Dvojím výpočtem délky byla částečně ověřena správnost výpočtu (pouze transformace mezi systémy, nikoliv výpočet Keplerových elementů). Výpočty byly provedeny v programu Octave. Zdrojový kód k výpočtům není přílohou technické zprávy (v případě potřeby bude zaslán).
V Kralupech nad Vltavou 29.10.2010
Jan Dolista (
[email protected])
