1

Wprowadzenie

W geodezji satelitarnej funkcjonuj ˛

a dwa układy współrz˛ednych zwi ˛

azane z ruchem obrotowym Ziemi Earth

Centered Earth Fixed (ECEF), oraz układ ekwinokcjalny. Ze wzgl˛edu na prost ˛

a form˛e matematyczn ˛

a a tak˙ze

niewielkie ilo´sci danych, orbity satelitów s ˛

a wysyłane w postaci elementów orbity keplerowskiej w układzie

ekwinokcjalnym. Orbity typu broadcast zostaj ˛

a zapisane w pliku RINEX n, w rozdzielczo´sci dwu godzinnej.

Dla wyznaczenia współrz˛ednych satelitów w układzie ECEF na dany moment nale˙zy dokona´c transformacji.
Poni˙zszy przykład ilustruje schemat przelicze´n mi˛edzy tymi dwoma układami współrz˛ednych.

2

Układy

2.1

Układ ekwinokcjalny

Współrz˛edne satelitów w układzie ekwinokcjalnym, składaj ˛

a si˛e z sze´sciu podstawowych parametrów or-

bity keplerowskiej:

Ω rektanscenzji w˛ezła wst˛epuj ˛

acego,

i inklinacji, ω argumentu perygeum, f anomali

prawdziwej,

a dłu˙zszej półosi orbity, e ekscentryczno´sci, W tym układzie o´s X skierowana jest w stron˛e

punktu równonocy, o´s

Z jest zgodna z osi ˛

a obrotu Ziemi, za to o´s

Y jest prostopadła do dwóch poprzednich

tworz ˛

ac układ prawoskr˛etny rys. 2a. ´Srodek układu znajduje si˛e w ´srodku mas Ziemi. Istotnymi punktami

zwi ˛

azanymi z obrit ˛

a s ˛

a: w˛ezeł wst˛epuj ˛

acy - jedno z miejsc przeci˛ecia si˛e płaszczyzny orbity z płaszczyzn ˛

a

równika, to w którym satelita porusza si˛e z południa na północ, oraz perygeum

P czyli punkt w którym

satelita znajduje si˛e najbli˙zej Ziemi.

2.2

Układ współrz˛ednych w płaszczy´znie orbity

Zgodnie z rys. 1 układ współrz˛ednych w płaszczy´znie orbity składa si˛e z osi

ξ skierowanej w stron˛e

perygeum, z osi

η skierowanej w stron˛e w˛ezła zst˛epuj ˛

acego (przeciwnie do w˛ezła wst˛epuj ˛

acego) oraz z

osi

ζ prostopadłej do płaszczyzny orbity. ´Srodek układu znajduje si˛e w centrum mas Ziemi.

2.3

Układ zwi ˛

azany z ruchem obrotowym Ziemi ECEF

W tym układzie o´s

X jest zwi ˛

azana z południkiem osiowym Greenwich który obraca si˛e wokół osi

Z układu

(osi obrotu Ziemi). Natomiast o´s

Y układu jest prostopadła do dwóch poprzednich i tworzy z nimi układ

prawoskr˛etny rys. 2a.

3

Współrz˛edne satelity w płaszczy´znie orbity

Zgodnie z rys. 1 współrz˛edne satelity w płaszczy´znie orbity (

ξ η ζ) uzyskujemy wi ˛

a˙z ˛

ac anomalie prawdziw ˛

a

z anomali ˛

a ekscentryczn ˛

a. Jak wynika z rys 1

ξ = r · cos(f) = a · cos(E) − a · e = a · (cos(E) − e)

(1)

η = r · sin(f) =

b

a

· sin(E) = b · sin(E) = a ·

p

1 − e

2

· sin(E)

(2)

ζ = 0;

(3)

Anomali˛e mimo´srodow ˛

a E otrzymamy, u˙zywaj ˛

ac równania Keplera (rozwi ˛

azanie iteracyjne):

E = µ + e · sin(E)

(4)

gdzie

µ to anomalia ´srednia czyli k ˛

at jaki ´srednio w czasie

t − t

0

zakre´sli promie´n wodz ˛

acy

satelity.Anomalia ´srednia

µ uzale˙zniona jest od ´sredniej pr˛edko´sci k ˛

atowej satelity:

n =

T

2 · π

=

pGM/a

3

(5)

1

x

h

u

0

x

S’

S

f

E

r

(-ae,a)

(-ae,b)

(-ae,0)

(0,0)

h

(a,0)

Rysunek 1: Układ współrz˛ednych w płaszczy´znie orbity, S - poło˙zenie satelity na orbicie

µ(t) = n · (t − t

0

)

(6)

4

Przeliczenia

4.1

Keplerowskie parametry ruchu satelitów

Aby wykona´c transformacj˛e

X, Y , Z na ξ, η, ζ nale˙zy wykona´c trzy obroty: wokół osi Z o k ˛

at rektanscenzji

w˛ezła wst˛epuj ˛

acego

Ω rys. 2a przy pomocy macierzy rotacji R

z

(Ω) w ten sposób przechodzimy do układu

przej´sciowego

x

′

,

y

′

,

z

′

, rys. 2b kolejno przy pomocy macierzy rotacji

R

x

(i) o k ˛

at inklinacji

i, przechodz ˛

ac

do układu

x”, y”, z”, a nast˛epnie wykonuj ˛

ac obrót przy pomocy macierzy rotacji

R

z

(ω), wokół osi z” rys.

2c. Na rys. 2, zaznaczono kolejne poło˙zenia transformowanego układu.

i

W

x

y

z

x’’

x’

y’

y’’

P

x

0

i

i

v

W

z’’

z’

rzu

t o

rb

ity

sa

te

lit

y

rzut

rów

nik

a

z

h

W

x

y

z

x’

y’

P

x

0

i

i

W

z’

rzu

t o

rb

ity

sa

te

lit

y

rzut

rów

nik

a

z

h

x

y

z

P

x

0

rzu

t o

rb

ity

sa

te

lit

y

rzut

rów

nik

a

z

h

Rysunek 2: Zestawienie kolejnych obrotów pomi˛edzy układami

ω, i, Ω

2

W notacji macierzowej zapiszemy:





X

s

(t

j

)

Y

s

(t

j

)

Z

s

(t

j

)





= R

z

(−Ω) · R

x

(−i) · R

z

(−ω) ·





r · cos(f)
r · sin(f)

0





(7)

gdzie, macierze rotacji

R

z

,

R

x

, dla dowolnego k ˛

ata

K to

R

x

=





1

0

0

0

cos(K)

sin(K)

0

−sin(K )

cos(K)





(8)

R

z

=





cos(K)

sin(K) 0

−sin(K )

cos(K)

0

0

0

1





(9)

Znaki − w (7) zwi ˛azane s ˛a z przej´sciem odwrotnym (tzn. ξ, η, ζ na X, Y , Z).

4.2

Poprawki do keplerowskich parametrów orbity

Poniewa˙z rzeczywisty ruch satelity odbiega od teorii keplerowskiej, w depeszy nawigacyjnej wyst˛epuj ˛

a tak˙ze

poprawki do sze´sciu podstawowych parametrów orbity. Poprawki te nale˙zy uwzgl˛edni´c 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˛edu nale˙zy najpierw obliczy´c ile czasu
min˛eło od czasu efemeryd

(t − t

oe

)

t

j

= t − t

oe

(10)

anomalia ´srednia w czasie

t

j

µ = µ

0

+ (

pGM/a

3

+ δn) · t

j

(11)

GM = 3.986005 ·10

14

m

3

/s

2

rozwi ˛

azanie iteracyjne równania Keplera (4),

nast˛epnie wyznaczy´c anomali˛e prawdziw ˛

a:

f

j

= arctan

 η

ξ



= arctan

√

1 − e

2

· sin(E

j

)

cos(E

j

) − e

!

(12)

rektascenzj˛e w˛ezła wst˛epuj ˛

acego:

Ω

j

= Ω

0

+ ( ˙

Ω

0

− ω

e

) · t

j

− ω

e

· t

oe

(13)

ω

e

= 7.292115147 · 10

−5

rad/s

Argument perygeum (tego parametru nie poprawiamy ze wzgl˛edu na cało´sciow ˛

a poprawk˛e do argumentu

szeroko´sci).

ω

j

= ω + f

i

+ C

ωC

· cos(2 · (ω + f

j

)) + C

ωS

· cos(2 · (ω + f

j

))

(14)

Odległo´s´c radialna(długo´s´c promienia wodz ˛

acego)

r

j

= a · (1 − e · cos(E

j

) + C

rC

· cos(2 · (ω + f

j

)) + C

rS

· cos(2 · (ω + f

j

))

(15)

inklinacj˛e:

i

j

= i

0

+ idot

j

+ C

iC

· cos(2 · (ω + f

j

)) + C

iS

· cos(2 · (ω + f

j

))

(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 ,

´Sledzi´nski J., Geodezja satelitarna

6

Przykład obliczeniowy

W kolejnej cz˛e´sci zrealizowane zostanie przeliczenie parametrów orbit kepplerowskich na współrz˛edne
ECEF. Do zadania zostały u˙zyte dane z 11 11 2008 roku ze stacji WROC. Zaprezentowano mo˙zliwo´s´c wyz-
naczenia orbity w dowolnym momencie pomi˛edzy kolejnymi epokami

t

oe

.

6.1

Pobieranie danych

Sie´c stacji permanentnych obejmuje swoim zasi˛egiem cał ˛

a kul˛e ziemsk ˛

a, jedna ze stacji IGS znajduje si˛e w

budynku Wydziału In˙zynierii Kształtowania ´Srodowiska i Geodezji UP we Wrocławiu. Stacja ma symbol
kodowy WROC. Dane pochodz ˛

ace z tej stacji mo˙zna pobra´c ze strony:

http://igs.bkg.bund.de/, zakładka: Download: IGS/obs/2008/doy/wroc[doy]0.08n.Z, otworzy´c, klikn ˛

a´c

showfile, skopiowa´c zawarto´s´c do notatnika.

Produkty sieci IGS, w tym precyzyjne orbity,

mo˙zna sci ˛

agn ˛

a´c z tego samego serwera:

http://igs.bkg.bund.de/,zakładka:

Download:

IGS/products/Tydzie´n GPS/Dzie´n tygodnia/igs[Tydzie´n

GPS][dzie´n tygodnia].sp3.Z, otworzy´c, klikn ˛

a´c showfile, skopiowa´c zawarto´s´c do notatnika.

6.2

Dane

Dane znajduj ˛

ace si˛e w tabeli 1.

pochodz ˛

a z pliku WROC3150.08n z fragmentu zamieszczonego

poni˙zej.Cz˛e´s´c 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˛e´s´c 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

10

n

Satelita nr

7

Rekord w pliku *.08n

0

Czas zegara satelity toc

[08 11 11 14 00 0.0]

Bł ˛

ad ´sredni zegara satelity a0

2.311961725354D-05

Dryft zegara satelity a1

3.410605131648D-13

Pr˛edko´s´c dryftu zegara satelity a2

0.000000000000D+00

Rekord w pliku *.08n

1

IOED

4.300000000000D+01

Poprawka do długo´sci promienia wodz ˛

acego

r − C

rs

-1.465625000000D+01

Poprawka do ruchu ´sredniego

δn

4.187317159676D-09

Anomalia ´srednia

µ

0

w momencie przej´scia satelity przez w˛ezeł wst˛epuj ˛

acy

-1.114328016966D+00

Rekord w pliku *.08n

2

Poprawka do anomali ´sredniej

C

uc

-9.424984455109D-07

Ekscentryczno´s´c

e

2.313868375495D-03

Poprawka do anomali ´sredniej

C

us

1.226924359798D-05

Pierwiastek z dłu˙zszej półosi orbity

√a

5.153689096451D+03

Rekord w pliku *.08n

3

Czas efemerydy pokładowej

t

oe

2.232000000000D+05

Poprawka do k ˛

ata inklinacji

C

ic

-9.313225746155D-08

Rektascenzja w˛ezła wst˛epuj ˛

acego

Ω

0

-2.831047113047D-02

Poprawka do k ˛

ata inklinacji

C

is

-3.539025783539D-08

Rekord w pliku *.08n

4

Inklinacja orbity

i

0

9.650716011630D-01

Poprawka do długo´sci promienia wodz ˛

acego

C

rc

1.459062500000D+02

Argument perygeum

ω

2.829229163768D+00

Pr˛edko´s´c zmiany rektascencji w˛ezła wst˛epuj ˛

acego ˙

Ω

-7.744251462327D-09

Rekord w pliku *.08n

5

Pr˛edko´s´c zmiany k ˛

ata nachylenia orbity

(idot) ˙i

5.432368999081D-10

Do kontroli oblicze´n wykorzystano informacje o poło˙zeniu satelitów z pliku igs15052.sp3. Pozycje

satelitów w układzie ECEF s ˛

a wyznaczone z interwałem 15 minutowym. Plik jest dost˛eny po około tygodniu

od momentu obserwacji (dane z postprocessingu). Cz˛e´s´c 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˛e´s´c zawieraj ˛

aca dane o poło˙zeniu 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˛edne 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˛edko´s´c obrotowa Ziemi

w

e

= 7.2921151467D − 005

(18)

Dłu˙zsza póło´s orbity satelity

a = (

√a)

2

= 26560511.302

(19)

´srednia pr˛edko´s´c k ˛

atowa satelity

n

0

=

 GM

(a

3

)



1

/2

= 1.4585D − 004

(20)

Czas na który chcemy policzy´c współrz˛edne satelity (w sekundach tygodnia GPS)
t = [08 11 11 16 00 0.0]
t =

2

d

16

h

00

m

0

s

t = 230400 s
Czas zegara satelity (w sekundach tygodnia GPS)
t

0

c = [08111114000.0]

t =

2

d

14

h

00

m

0

s

t = 223200 s
Poprawka do czasu zegata satelity

dt = a0 + a1 · (t − t

oc

) + a2 · (t − t

oc

)

2

= 2.3122D − 005

(21)

Czas na jaki chcemy wyznaczy´c współrz˛edne satelity z uwzgl˛ednieniem bł˛edów zegara satelity

t

j

= t − t

oe

− dt = 7.200 − 2.3122D − 005

(22)

Czas

t

oc

mo˙zna uzyska´c z pliku *.n w sekundach tygodnia GPS, mo˙ze zdarzy´c si˛e, ˙ze

(t − t

oc

) czas

odbiega od zało˙zonego wzorca (´srodek tygodnia 302400), wypadaj ˛

ac w chwili po ko ´ncu tygodnia lub przed

pocz ˛

atkiem. Nale˙zy poprawi´c wtedy: je´sli

(t − t

oc

) > 302400 to (t − t

oc

) − 604800; natomiast je´sli

(t − t

oc

) < −302400 to (t − t

oc

) + 604800;

pr˛edko´s´c k ˛

atowa satelity

n = n

0

+ δn = 1.4586D − 004

(23)

Anomalia ´srednia

µ

j

= µ

0

+ n · t

j

= −0.064158915

(24)

Anomalia mimo´srodowa pierwsza iteracja

E

0

j

= µ

j

+ e · sin(µ

j

)

(25)

Kolejne iteracje

E

1

j

= µ

j

+ e · sin(E

0

j

)

(26)

Warto´sci kolejnych przybli˙ze´n znajduj ˛

a si˛e w tabeli.

Tablica 2: Warto´sci

E anomalii mimo´srodowej z kolejnych iteracji

E

0

1.000000000

E

1

-0.062211862

E

2

-0.064302773

E

3

-0.064307601

E

4

-0.064307612

E

5

-0.064307612

7

Anomali˛e prawdziw ˛

a

f (12)- nale˙zy znale´z´c stosuj ˛

ac sposób oblicze´n podobny do wyznaczania azy-

mutów ze współrz˛ednych (uwzgl˛edniaj ˛

ac czwartaki).

f = 6.218728826

(27)

Suma k ˛

atów argumentu perygeum

ω oraz anomalii prawdziwej f daje w efekcie argument szeroko´sci u.

u = f + ω = 9.047957989

(28)

Wyznaczenie poprawek do argumentu szeroko´sci, długo´sci wektora wodz ˛

acego oraz inklinacji

du

j

= C

us

· sin(2 · u) + C

uc

· cos(2 · u) = −9.083085863D − 006

(29)

dr

j

= C

rs

· sin(2 · u) + C

rc

· cos(2 · u) = 1.164244966D + 002

(30)

di

j

= C

is

· sin(2 · u) + C

ic

· cos(2 · u) + idot · t

j

= 3.867610821D − 006

(31)

Uwzgl˛ednienie poprawek do argumentu szeroko´sci, długo´sci wektora wodz ˛

acego, inklinacji oraz rektascen-

zji.

u

j

= u + du

j

= 9.047948906

(32)

r

j

= a · (1 − e · cos(Ek)) + dr

j

= 2.649929723D + 007

(33)

i

j

= i

o

+ di

j

= 0.965075468

(34)

Ω

j

= Ω

0

+ ( ˙

Ω − ω

e

) · t

j

− ω

e

· t

oe

= −16.829399526

(35)

Ze wzoru (7) - opuszczaj ˛

ac ostatni ˛

a rotacj˛e, poniewa˙z k ˛

at

ω został przeliczony na k ˛

at

u (28). Ostatecznie w

wyniku otrzymano:

X = 5702699.51
Y = -24605519.25
Z = 8016258.12
Wyniki skontrolowano z plikiem SP3 (igs15052.sp3) otrzymuj ˛

ac odchyłki.

dX = -1.97
dY = -0.33
dZ = 1.18

8