cwiczenie nr1


1 Wprowadzenie
W geodezji satelitarnej funkcjonują dwa układy współrzędnych związane z ruchem obrotowym Ziemi Earth
Centered Earth Fixed (ECEF), oraz układ ekwinokcjalny. Ze względu na prostą formę matematyczną a także
niewielkie ilości danych, orbity satelitów są wysyłane w postaci elementów orbity keplerowskiej w układzie
ekwinokcjalnym. Orbity typu broadcast zostają zapisane w pliku RINEX n, w rozdzielczości dwu godzinnej.
Dla wyznaczenia współrzędnych satelitów w układzie ECEF na dany moment należy dokonać transformacji.
Poniższy przykład ilustruje schemat przeliczeń między tymi dwoma układami współrzędnych.
2 Układy
2.1 Układ ekwinokcjalny
Współrzędne satelitów w układzie ekwinokcjalnym, składają się z sześciu podstawowych parametrów or-
bity keplerowskiej: &! rektanscenzji węzła wstępującego, i inklinacji,  argumentu perygeum, f anomali
prawdziwej, a dłuższej półosi orbity, e ekscentryczności, W tym układzie oś X skierowana jest w stronę
punktu równonocy, oś Z jest zgodna z osią obrotu Ziemi, za to oś Y jest prostopadła do dwóch poprzednich
tworząc układ prawoskrętny rys. 2a. Środek układu znajduje się w środku mas Ziemi. Istotnymi punktami
związanymi z obritą są: węzeł wstępujący - jedno z miejsc przecięcia się płaszczyzny orbity z płaszczyzną
równika, to w którym satelita porusza się z południa na północ, oraz perygeum P czyli punkt w którym
satelita znajduje się najbliżej Ziemi.
2.2 Układ współrzędnych w płaszczyznie orbity
Zgodnie z rys. 1 układ współrzędnych w płaszczyznie orbity składa się z osi  skierowanej w stronę
perygeum, z osi  skierowanej w stronę węzła zstępującego (przeciwnie do węzła wstępującego) oraz z
osi ś prostopadłej do płaszczyzny orbity. Środek układu znajduje się w centrum mas Ziemi.
2.3 Układ związany z ruchem obrotowym Ziemi ECEF
W tym układzie oś X jest związana z południkiem osiowym Greenwich który obraca się wokół osi Z układu
(osi obrotu Ziemi). Natomiast oś Y układu jest prostopadła do dwóch poprzednich i tworzy z nimi układ
prawoskrętny rys. 2a.
3 Współrzędne satelity w płaszczyznie orbity
Zgodnie z rys. 1 współrzędne satelity w płaszczyznie orbity (  ś) uzyskujemy wiążąc anomalie prawdziwą
z anomalią ekscentryczną. Jak wynika z rys 1
 = r cos(f) = a cos(E) - a e = a (cos(E) - e) (1)

b
 = r sin(f) = sin(E) = b sin(E) = a 1 - e2 sin(E) (2)
a
ś = 0; (3)
Anomalię mimośrodową E otrzymamy, używając równania Keplera (rozwiązanie iteracyjne):
E = + e sin(E) (4)
gdzie to anomalia średnia czyli kąt jaki średnio w czasie t - t0 zakreśli promień wodzący
satelity.Anomalia średnia uzależniona jest od średniej prędkości kątowej satelity:

T
n = = GM/a3 (5)
2 Ą
1

(-ae,a)
S
(-ae,b)

S
r


f
E

(-ae,0) (0,0) (a,0)
Rysunek 1: Układ współrzędnych w płaszczyznie orbity, S - położenie satelity na orbicie
(t) = n (t - t0) (6)
4 Przeliczenia
4.1 Keplerowskie parametry ruchu satelitów
Aby wykonać transformację X, Y , Z na , , ś należy wykonać trzy obroty: wokół osi Z o kąt rektanscenzji
węzła wstępującego &! rys. 2a przy pomocy macierzy rotacji Rz(&!) w ten sposób przechodzimy do układu
przejściowego x2 , y2 , z2 , rys. 2b kolejno przy pomocy macierzy rotacji Rx(i) o kąt inklinacji i, przechodząc
do układu x , y , z , a następnie wykonując obrót przy pomocy macierzy rotacji Rz(), wokół osi z rys.
2c. Na rys. 2, zaznaczono kolejne położenia transformowanego układu.
z z z
z z
y 
z 



y
y


P P
P
i
i
i


0 0
0
y y
y


i
i
x
x
x 
x x
x
Rysunek 2: Zestawienie kolejnych obrotów pomiędzy układami , i, &!
2
a
k
a
a
i
k
k
n
i
i
n
w
n
ó
w
r
w
ó
t
ó
r
r
u
t
z
t
u
r
u
z
r
z
r
y
t
i
y
l
y
t
i
t
e
l
i
t
l
e
a
t
e
s
t
a
a
s
y
s
t
i
y
t
y
b
i
t
r
i
b
o
r
b
r
t
o
o
u
t
z
t
u
r
u
z
r
z
r
W notacji macierzowej zapiszemy:
ł łł ł łł
Xs(tj) r cos(f)
ł ł ł ł
Ys(tj) = Rz(-&!) Rx(-i) Rz(-) r sin(f) (7)
Zs(tj) 0
gdzie, macierze rotacji Rz, Rx, dla dowolnego kąta K to
ł łł
1 0 0
ł ł
Rx = 0 cos(K) sin(K) (8)
0 -sin(K) cos(K)
ł łł
cos(K) sin(K) 0
ł ł
Rz = -sin(K) cos(K) 0 (9)
0 0 1
Znaki - w (7) związane są z przejściem odwrotnym (tzn. , , ś na X, Y , Z).
4.2 Poprawki do keplerowskich parametrów orbity
Ponieważ rzeczywisty ruch satelity odbiega od teorii keplerowskiej, w depeszy nawigacyjnej występują także
poprawki do sześciu podstawowych parametrów orbity. Poprawki te należy uwzględnić przed wstawieniem
ostatecznych parametrów do równania (7). Efemerydy satelitów zostały wyznaczone na pewien konkretny
moment (o pełnej godzinie z interwałem co dwie godziny), z tego względu należy najpierw obliczyć ile czasu
minęło od czasu efemeryd (t - toe)
tj = t - toe (10)
anomalia średnia w czasie tj

= 0 + ( GM/a3 + n) tj (11)
GM = 3.986005 1014m3/s2
rozwiązanie iteracyjne równania Keplera (4),
następnie wyznaczyć anomalię prawdziwą:

"

 1 - e2 sin(Ej)
fj = arctan = arctan (12)
 cos(Ej) - e
rektascenzję węzła wstępującego:
Ł
&!j = &!0 + (&!0 - e) tj - e toe (13)
e = 7.292115147 10-5rad/s
Argument perygeum (tego parametru nie poprawiamy ze względu na całościową poprawkę do argumentu
szerokości).
j =  + fi + CC cos(2 ( + fj)) + CS cos(2 ( + fj)) (14)
Odległość radialna(długość promienia wodzącego)
rj = a (1 - e cos(Ej) + CrC cos(2 ( + fj)) + CrS cos(2 ( + fj)) (15)
inklinację:
ij = i0 + idotj + CiC cos(2 ( + fj)) + CiS cos(2 ( + fj)) (16)
3
5 Materiały
Szewczyk J., Góral W., Zastosowanie technologii GPS w precyzyjnych pomiarach deformacji ,
Strang G., Borre K., Linear Algebra, Geodesy, and GPS ,
Lamparski J., NAVSTAR GPS od teorii do praktyki ,
Śledziński J., Geodezja satelitarna
6 Przykład obliczeniowy
W kolejnej części zrealizowane zostanie przeliczenie parametrów orbit kepplerowskich na współrzędne
ECEF. Do zadania zostały użyte dane z 11 11 2008 roku ze stacji WROC. Zaprezentowano możliwość wyz-
naczenia orbity w dowolnym momencie pomiędzy kolejnymi epokami toe.
6.1 Pobieranie danych
Sieć stacji permanentnych obejmuje swoim zasięgiem całą kulę ziemską, jedna ze stacji IGS znajduje się w
budynku Wydziału Inżynierii Kształtowania Środowiska i Geodezji UP we Wrocławiu. Stacja ma symbol
kodowy WROC. Dane pochodzące z tej stacji można pobrać ze strony:
http://igs.bkg.bund.de/, zakładka: Download: IGS/obs/2008/doy/wroc[doy]0.08n.Z, otworzyć, kliknąć
showfile, skopiować zawartość do notatnika.
Produkty sieci IGS, w tym precyzyjne orbity, można sciągnąć z tego samego serwera:
http://igs.bkg.bund.de/,zakładka: Download: IGS/products/Tydzień GPS/Dzień tygodnia/igs[Tydzień
GPS][dzień tygodnia].sp3.Z, otworzyć, kliknąć showfile, skopiować zawartość do notatnika.
6.2 Dane
Dane znajdujące się w tabeli 1. pochodzą z pliku WROC3150.08n z fragmentu zamieszczonego
poniżej.Część nagłówkowa pliku:
2.11 NAVIGATION DATA RINEX VERSION / TYPE
SPIDER V2,2,0,2479 IGG 2008 11 11 22:54 PGM / RUN BY / DATE
1.2107D-08 -7.4506D-09 -1.1921D-07 5.9605D-08 ION ALPHA
9.8304D+04 -8.1920D+04 -1.9661D+05 4.5875D+05 ION BETA
-5.587935447693D-09-2.575717417130D-14 405504 1505 DELTA-UTC: A0,A1,T,W
14 LEAP SECONDS
END OF HEADER
Część z danymi:
7 08 11 11 14 00 0.0 2.311961725354D-05 3.410605131648D-13 0.000000000000D+00
4.300000000000D+01-1.465625000000D+01 4.187317159676D-09-1.114328016966D+00
-9.424984455109D-07 2.313868375495D-03 1.226924359798D-05 5.153689096451D+03
2.232000000000D+05-9.313225746155D-08-2.831047113047D-02-3.539025783539D-08
9.650716011630D-01 1.459062500000D+02 2.829229163768D+00-7.744251462327D-09
5.432368999081D-10 1.000000000000D+00 1.505000000000D+03 0.000000000000D+00
2.000000000000D+00 0.000000000000D+00-1.071020960808D-08 4.300000000000D+01
2.232000000000D+05 0.000000000000D+00
4
Tablica 1: Zestawienie danych z pliku WROC3150.08n. Litera D przy liczbie oznacza 10n
Satelita nr 7
Rekord w pliku *.08n 0
Czas zegara satelity toc [08 11 11 14 00 0.0]
Błąd średni zegara satelity a0 2.311961725354D-05
Dryft zegara satelity a1 3.410605131648D-13
Prędkość dryftu zegara satelity a2 0.000000000000D+00
Rekord w pliku *.08n 1
IOED 4.300000000000D+01
Poprawka do długości promienia wodzącego r - Crs -1.465625000000D+01
Poprawka do ruchu średniego n 4.187317159676D-09
Anomalia średnia 0 w momencie przejścia satelity przez węzeł wstępujący -1.114328016966D+00
Rekord w pliku *.08n 2
Poprawka do anomali średniej Cuc -9.424984455109D-07
Ekscentryczność e 2.313868375495D-03
Poprawka do anomali średniej Cus 1.226924359798D-05
"
Pierwiastek z dłuższej półosi orbity a 5.153689096451D+03
Rekord w pliku *.08n 3
Czas efemerydy pokładowej toe 2.232000000000D+05
Poprawka do kąta inklinacji Cic -9.313225746155D-08
Rektascenzja węzła wstępującego &!0 -2.831047113047D-02
Poprawka do kąta inklinacji Cis -3.539025783539D-08
Rekord w pliku *.08n 4
Inklinacja orbity i0 9.650716011630D-01
Poprawka do długości promienia wodzącego Crc 1.459062500000D+02
Argument perygeum  2.829229163768D+00
Ł
Prędkość zmiany rektascencji węzła wstępującego &! -7.744251462327D-09
Rekord w pliku *.08n 5
Ł
Prędkość zmiany kąta nachylenia orbity (idot) i 5.432368999081D-10
Do kontroli obliczeń wykorzystano informacje o położeniu satelitów z pliku igs15052.sp3. Pozycje
satelitów w układzie ECEF są wyznaczone z interwałem 15 minutowym. Plik jest dostęny po około tygodniu
od momentu obserwacji (dane z postprocessingu). Część nagłówkowa pliku:
#cP2008 11 11 0 0 0.00000000 96 ORBIT IGS05 HLM IGS
## 1505 172800.00000000 900.00000000 54781 0.0000000000000
+ 32 G01G02G03G04G05G06G07G08G09G10G11G12G13G14G15G16G17
+ G18G19G20G21G22G23G24G25G26G27G28G29G30G31G32 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
++ 0 3 2 2 2 2 3 3 3 2 2 2 3 3 2 3 3
++ 3 3 2 2 3 3 3 3 2 2 3 3 2 3 3 0 0
++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
%c G cc GPS ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc
%c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc
%f 1.2500000 1.025000000 0.00000000000 0.000000000000000
%f 0.0000000 0.000000000 0.00000000000 0.000000000000000
5
%i 0 0 0 0 0 0 0 0 0
%i 0 0 0 0 0 0 0 0 0
/* FINAL ORBIT COMBINATION FROM WEIGHTED AVERAGE OF:
/* cod emr esa gfz grg jpl mit ngs sio
/* REFERENCED TO IGS TIME (IGST) AND TO WEIGHTED MEAN POLE:
/* PCV:IGS05_1502 OL/AL:FES2004 NONE Y ORB:CMB CLK:CMB
Część zawierająca dane o położeniu satelitów:
2008 11 11 16 0 0.00000000
*
PG01 6819.100318 22684.542613 -11726.615294 999999.999999
PG02 -17329.253956 2856.130178 -19947.396290 173.908276 6 7 5 143
PG03 20102.586173 4117.986849 16755.586569 309.673456 5 4 4 118
PG04 -11279.751100 -11738.202242 -21045.153370 -229.059453 6 7 6 168
PG05 -13724.580949 15036.800475 -17288.041562 -382.023204 6 7 5 167
PG06 20292.734384 8594.496315 15018.444743 145.315970 5 5 6 97
PG07 5702.701494 -24605.518923 8016.256940 23.140425 7 10 6 136
PG08 -3754.712097 -18439.057161 18461.527494 -177.581258 4 7 4 161
PG09 -14336.241697 20967.943635 6906.222068 24.380977 9 9 9 76
PG10 -26725.180394 931.949808 -565.870759 -7.101661 7 7 4 155
PG11 11615.427277 -22551.724720 7073.884357 15.874061 8 8 6 138
PG12 -18863.802451 11492.080442 -14608.599782 -341.232918 8 10 8 97
PG13 -2207.463435 -22060.025269 -14740.545862 276.836701 8 10 7 121
PG14 16085.267891 20868.063299 -3705.901104 -195.869876 6 7 6 93
PG15 -14494.281590 5746.985824 21537.437156 -209.170250 7 7 6 161
PG16 25725.167399 2974.040717 -6330.842531 95.982310 8 7 7 133
PG17 -15646.303743 -21510.484459 -1729.366793 40.383722 5 8 8 109
PG18 2321.089161 15546.041710 21633.147203 -114.009287 7 9 5 118
PG19 14056.482313 -6262.646171 21723.703717 35.585619 5 5 2 167
PG20 15592.381753 -13630.642791 -16722.241393 99.393020 7 9 6 131
PG21 -1294.818127 23109.759466 12670.360164 43.009188 5 6 7 175
PG22 15591.952708 12368.989530 17780.047243 209.478767 2 5 5 152
PG23 5637.964924 -16307.046109 -20163.708828 379.693574 9 10 7 115
PG24 -13017.135109 23320.427687 -341.221443 140.238800 8 9 6 181
PG25 8558.692736 -25236.036504 1066.432266 15.431173 5 8 4 168
PG26 -16582.840398 -1768.208949 20310.699333 14.048462 6 6 7 145
PG27 1429.315846 -24149.088523 11518.417269 232.011155 7 7 4 168
PG28 -14040.989291 -12665.283190 19183.217733 -23.480632 3 5 5 130
PG29 634.799620 21301.354863 -15773.708139 -30.471641 8 9 7 169
PG30 -4601.989205 15798.916457 -21194.796410 107.264033 6 7 5 138
PG31 12315.062801 8637.400382 -21674.224103 -37.838442 7 9 5 146
PG32 22283.915124 -7747.817278 -11563.079862 314.328065 6 7 8 125
Współrzędne precyzyjne satelity nr 7 z SP3 to:
X = 5702.701494
Y = -24605.518923
Z = 8016.256940
6.3 Obliczenia
Stałe: potencjał Ziemski
GM = 3.986005D + 14 (17)
6
prędkość obrotowa Ziemi
we = 7.2921151467D - 005 (18)
Dłuższa półoś orbity satelity
"
a = ( a)2 = 26560511.302 (19)
średnia prędkość kątowa satelity
1/2
GM
n0 = = 1.4585D - 004 (20)
(a3)
Czas na który chcemy policzyć współrzędne satelity (w sekundach tygodnia GPS)
t = [08 11 11 16 00 0.0]
t = 2d 16h 00m 0s
t = 230400 s
Czas zegara satelity (w sekundach tygodnia GPS)
t0c = [08111114000.0]
t = 2d 14h 00m 0s
t = 223200 s
Poprawka do czasu zegata satelity
dt = a0 + a1 (t - toc) + a2 (t - toc)2 = 2.3122D - 005 (21)
Czas na jaki chcemy wyznaczyć współrzędne satelity z uwzględnieniem błędów zegara satelity
tj = t - toe - dt = 7.200 - 2.3122D - 005 (22)
Czas toc można uzyskać z pliku *.n w sekundach tygodnia GPS, może zdarzyć się, że (t - toc) czas
odbiega od założonego wzorca (środek tygodnia 302400), wypadając w chwili po końcu tygodnia lub przed
początkiem. Należy poprawić wtedy: jeśli (t - toc) > 302400 to (t - toc) - 604800; natomiast jeśli
(t - toc) < -302400 to (t - toc) + 604800;
prędkość kątowa satelity
n = n0 + n = 1.4586D - 004 (23)
Anomalia średnia
j = 0 + n tj = -0.064158915 (24)
Anomalia mimośrodowa pierwsza iteracja
0
Ej = j + e sin(j) (25)
Kolejne iteracje
1 0
Ej = j + e sin(Ej ) (26)
Wartości kolejnych przybliżeń znajdują się w tabeli.
Tablica 2: Wartości E anomalii mimośrodowej z kolejnych iteracji
E0 1.000000000
E1 -0.062211862
E2 -0.064302773
E3 -0.064307601
E4 -0.064307612
E5 -0.064307612
7
Anomalię prawdziwą f (12)- należy znalezć stosując sposób obliczeń podobny do wyznaczania azy-
mutów ze współrzędnych (uwzględniając czwartaki).
f = 6.218728826 (27)
Suma kątów argumentu perygeum  oraz anomalii prawdziwej f daje w efekcie argument szerokości u.
u = f +  = 9.047957989 (28)
Wyznaczenie poprawek do argumentu szerokości, długości wektora wodzącego oraz inklinacji
duj = Cus sin(2 u) + Cuc cos(2 u) = -9.083085863D - 006 (29)
drj = Crs sin(2 u) + Crc cos(2 u) = 1.164244966D + 002 (30)
dij = Cis sin(2 u) + Cic cos(2 u) + idot tj = 3.867610821D - 006 (31)
Uwzględnienie poprawek do argumentu szerokości, długości wektora wodzącego, inklinacji oraz rektascen-
zji.
uj = u + duj = 9.047948906 (32)
rj = a (1 - e cos(Ek)) + drj = 2.649929723D + 007 (33)
ij = io + dij = 0.965075468 (34)
Ł
&!j = &!0 + (&! - e) tj - e toe = -16.829399526 (35)
Ze wzoru (7) - opuszczając ostatnią rotację, ponieważ kąt  został przeliczony na kąt u (28). Ostatecznie w
wyniku otrzymano:
X = 5702699.51
Y = -24605519.25
Z = 8016258.12
Wyniki skontrolowano z plikiem SP3 (igs15052.sp3) otrzymując odchyłki.
dX = -1.97
dY = -0.33
dZ = 1.18
8


Wyszukiwarka