DOLOČITEV ABSOLUTNEGA POLOŽAJA GPS-SPREJEMNIKA IZ KODNIH OPAZOVAN) ABSOLUTE GPS POSITIONING FROM CODE OBSERVATIONS Polona Pavlovčič Prešeren, Bojan Stopar UDK: 528.28 POVZETEK V prispevku obravnavamo nekaj algoritmov določitve absolutnega položaja in popravka sprejemnikove ure iz podatkov opazovanj kodnih psevdorazdalj. Problem razdelimo v direktne in iterativne rešitve za primer, ko imamo na voljo podatke opazovanj štirih ali več satelitov. Dodatno predstavimo tudi rešitev določitve položaja s tremi psevdorazdaljami (rešitev z apriori znano višino), pogostega problema v avtomobilski navigaciji. Tekom posameznih izpeljav obravnavamo prednosti metod in geometrične pogoje, ki vodijo od več rešitvam k najbolj verjetni. Namen prispevka je izbrati optimalno metodo za določitev začetnih vrednosti neznank za rekurzivno obdelavo s Kalmanovim filtriranjem oziroma za določitev apriori vrednosti za nadaljnjo obdelavo s faznimi opazovanji. Metode izračuna primerjamo na izbranih trenutkih sprejema GPS-signala, z uporabo oddanih efemerid in ob modeliranju troposferske refrakcije in aberacije Zemlje v za ta namen izdelanem programu v okolju Matlab R14. Z nadaljnjimi analizami pokažemo, da na končne rezultate vpliva tudi število satelitov in čas trajanja opazovanj. KLJUČNE BESEDE GPS (Global Positioning System), absolutni položaj, enačba psevdorazdalje, pogrešek ure, rešitev z apriori znano višino, Kalmanovo filtriranje, oddane efemeride, troposferska Klasifikacija prispevka po COBISS-u: 1.02 ABSTRACT In the article we present the single point positioning procedures from pseudoranges for obtaining the absolute position of the GPS antenna and its clock offset. First we define the direct and iterative solutions from four or more pseudoranges, further we derive a solution from three pseudoranges used in situations with obstructions to the receiver-to-satellite line of sights, for example in vehicle navigation (altitude hold mode). We present the theoretic demonstration of preferences of each method and discuss the geometric conditions leading from several solutions to the most likely one. The main purpose of this study is finding the optimal computation method to attain a well-based initial solution guess for the recursive method of computation - Kalman filtering, and to obtain quality apriori values for further carrier-phase processing. In order to achieve a well-grounded method comparison in the Matlab R14 program, we apply all methods on the same set of observation data and in connection with broadcast ephemerides considering troposphere refraction models and Earth aberration. Additional analyses show that the number of satellites and observation time interval have a significant effect on the accuracy improvement. KEY WORDS GPS, single point positioning, pseudoranges, clock offsets, altitude hold mode, Kalman filtering, broadcast ephemerides, troposphere refraction models, Earth aberration Q 1 UVOD 55 Ob začetku obdelave GPS-opazovanj za potrebe geodezije je treba določiti približne vrednosti določenih parametrov, ki so izhodišče za nadaljnjo obdelavo. Ena takih nalog je določitev absolutnega položaja GPS-sprejemnika iz kodnih opazovanj, ki je lahko zanimiva tudi širšemu negeodetskemu krogu bralcev. Lahko bi rekli, da je določitev absolutnega položaja GPS-sprejemnika iz kodnih opazovanj splošno najbolj uporabljena metoda GPS-izmere v realnem času, saj je dovolj, da imamo na voljo najcenejše GPS-sprejemnike. V prispevku se bomo podrobneje seznanili z različnimi metodami izračuna, ki nas privedejo od opazovanih psevdorazdalj do položaja antene GPS-sprejemnika. Metode določitve absolutnega položaja GPS-sprejemnika delimo v direktne, iterativne in rekurzivne. Z direktnimi metodami algebraično rešujemo problem neznanega položaja GPS-sprejemnika in sprejemnikove ure iz psevdorazdalj in položajev GPS-satelitov, ne da bi vnaprej poznali približni položaj opazovališča.1 V iterativnih metodah problem rešujemo z več koraki; znane so rešitve z oziroma brez linearizacije enačb psevdorazdalj, do rezultata pa pridemo po večkratni obdelavi opazovanj istega trenutka. V primeru iterativnih metod, ki slonijo na linearizaciji psevdorazdalj, moramo poznati približne vrednosti položaja opazovališča. Nadgradnja teh metod so rekurzivne metode, kamor uvrščamo postopek Kalmanovega filtriranja. Tu na osnovi zaporednega niza opazovanj poizkušamo najti optimalno rešitev z minimiziranjem variance napake opazovanj. Tudi v tem primeru moramo imeti na voljo približne koordinate opazovališča. Zaradi obsežnosti rešitve s Kalmanovim filtriranjem se bomo omejili le na določitev položaja opazovališča iz opazovanj enega trenutka. V določenih primerih, na primer v avtomobilski navigaciji, nimamo vedno na voljo opazovanj do štirih satelitov. Zato omenjamo tudi metode določitve položaja sprejemnika na osnovi treh psevdorazdalj. V danem primeru tekom obdelave izhajamo iz predpostavke, da se višinska komponenta položaja h ne spreminja tako hitro kot horizontalna (9, X)in jo zato obravnavamo kot znano, vrednost pa privzamemo iz prejšnjega koraka. Rešitev poznamo pod imenom metoda z znano apriori višino (angl. altitude hold mode) oziroma določitev 2D-položaja s tehnologijo GPS. Direktne metode, kjer v enem koraku pridemo do rešitve, imajo prednosti tudi zato, ker z njimi pospešimo hitrost obdelave GPS-opazovanj oziroma hitrost konvergiranja Kalmanovega filtra. Direktne algebraične rešitve določitve absolutnega položaja GPS-sprejemnika so znane že dolgo (Bancroft (1985), Fang (1986), Krause (1987), Kleusberg (1994)), v zadnjem času pa so jih izboljšali Abel in Chafee (1991, 1994), Grafarend in Shan (1996), Crocetto et al. (1998). Znani so tudi pristopi reševanja določitve absolutnega položaja z umetno nevronsko mrežo (Chansarkar, 2001; Jwo, 2004). 1.1 Določitev absolutnega položaja Metodo lahko uporabimo tako pri kodnih kot pri faznih opazovanjih. V prvem primeru bomo pridobili položaj sprejemnika z nekaj desetmetrsko točnostjo, v drugem primeru pa z nekaj 1 Opazovališče v širšem pomenu besede predstavlja točka, kjer izvajamo GP^-opaz^^vanja. Vedeti je treba, da tekom obdelave OPS-opazo^vanj razlikujemo točko, na katero se nanašajo GPS-opazovanja (fazni center GPS-antene), in trajno stabilizirano točko, za katero dejansko določamo položaj. Da lahko to točko določimo, merimo višino antene. Če je ne upošt^^amo, je opazo^vališče končna točka GPS-opaz^^anj -GPS-antena. centimetrsko - govorimo o tehniki natančne določitve absolutnega položaja PPP (angl. Precise Point Positioning). V obeh primerih je način izvedbe izmere enak: z enim GPS-sprejemnikom izvajamo opazovanja signala vsaj štirih GPS-satelitov (slika 1), le opazovane količine so različnega tipa. V prvem primeru bomo obdelovali kodne psevdorazdalje, v drugem primeru pa dodatno še opazovanja faze nosilnega valovanja. Princip določanja absolutnega položaja z GPS-opazovanji je podoben kot pri terestrični trilateraciji. Ko z opazovanji pridobimo razdalje do GPS-satelitov, nadalje položaj določimo s presekom krogel, katerih polmer je enak opazovanim razdaljam satelit-sprejemnik. Za neznano določitev položaja v 3D-prostoru bi morali imeti na voljo razdalje od najmanj treh satelitov, vendar šele četrta razdalja omogoča določitev razlike urinega stanja sprejemnikove ure (glede na nominalni GPS-čas). GPS-sprejemnik je tekom izmere lahko v gibanju, vendar je natančnost določitve položaja odvisna tudi od trajanja opazovanj na posameznem stojišču. Poleg tega na natančnost položaja vplivajo še drugi faktorji: kakovost GPS-sprejemnika, geometrijska razporeditev satelitov in sprejemnika (faktor GDOP), način obdelave (za boljšo določitev lahko dodatno vključimo tudi Kalmanovo filtriranje). Slika 1: Določitev absolutnega položaja na osnovi opazovanj GPS. 1.2 Določitev absolutnega položaja na osnovi kodnih psevdorazdalj V prispevku se omejujemo na določitev absolutnega položaja za potrebe navigacije (merska točnost položaja). V sprejemnikovi kodni zanki poteka spremljanje kode satelitskega signala. S primerjavo s satelita oddane kode in v sprejemniku generirane kode pridobimo navigacijsko sporočilo in kot stranski produkt vrednosti opazovanj psevdorazdalj. Ta opazovanja se nanašajo na urino stanje sprejemnika (sprejemnikova in satelitova ura nista usklajeni). Razdalja, ki jo je signal prepotoval od satelita do sprejemnika, je enaka produktu hitrosti svetlobe (c = 299792458 m/s) in časa potovanja signala: Pi = c(Ti - ti), (1.2.1) kjer je T. trenutek sprejema in t^ trenutek oddaje signala. Na hitrost potovanja signala vplivajo tudi troposferski in ionosferski pogoji. Zato so opazovanja, ki jih obdelujemo, le približne razdalje CN med sateliti do sprejemnikov - imenujemo jih psevdorazdalje: Pi =-,!{xj - xf' + {yi - y )2 + {zj - z)2 + c ■ dT - c ■ dti + Tropoi + lonoj + e,, kjer so: (1.2.2) x i = yi x = dT dt Tropo. lono. x y z znan vektor položaja satelita i neznan vektor položaja GPS-sprejemnika pogrešek sprejemnikove ure (neznana količina) pogrešek ure satelita i (dano v navigacijskem sporočilu ali pa v SP3-datotekah) vpliv troposferske refrakcije vpliv ionosferske refrakcije drugi neodstranjeni vplivi. Enačbo psevdorazdalje poenostavimo tako, da ne upoštevamo vpliva ionosfere in drugih neodstranjenih vplivov: P - x)2 +{yi - y)2 +{zi - z)2 + c ■ dT - c ■ dti + Tropoi. (1.2.3) Tretji člen izraza (1.2.3), to je c ■ dti, bomo določili iz podatkov efemerid2, pogrešek troposferske refrakcije Tropo pa določimo s pomočjo ustreznega modela. Za pogrešek satelitovih ur in troposferski vpliv popravljeno psevdorazdaljo matematično zapišemo kot: I^i - x)2 + (yi - y)2 +{zi - z)2 + c ■ dT = p, + b, (1.2.4) kjer razdaljo med i-tim satelitom in sprejemnikom označimo z p,, produkt c ■ dT pa označimo z b. Enačba (1.2.4), kjer nastopajo štiri neznanke (x, y, z, b), nam predstavlja izhodišče za nadaljnje določanje položaja GPS-opazovališča iz popravljenih psevdorazdalj z različnimi metodami. 2 MATEMATIČNI MODELI DOLOČITVE ABSOLUTNEGA POLOŽAJA IZ TREH KODNIH PSEVDORAZDALJ •o tig 55 2.1 Direktna rešitev - brez pogreška sprejemnikove ure Tu predpostavljamo, da nimamo pogreška urinega stanja sprejemnikove ure3 (b = 0), zato je metoda zgolj približna. Služi za izhodišče določanja položaja s štirimi ali več sateliti z iteracijskimi metodami. Za tri satelite tako zapišemo enačbo psevdorazdalje v najbolj poenostavljeni obliki, ^ Pogrešek satelitovih ur rekonstruiramo nazaj s pomočjo koeficientov polinoma, ki so podani v podatkih broadcast efemerid, ali pa vrednost direktno dobimo datotek s preciznimi efemeridami. 3 Pogrešek sprejemnikove ure 1 ns vodi do pogreška v izmerjeni psevdorazdalji približno 0,29 m. e. tj. brez vplivov: - +{yi - y +{zi - z i = 1,2,3 ter enačbe kvadriramo: -z - 2XiX - 2 yiy - 2ZiZ- + y2 + z2 - P,^ )= O, i- =1,2,3 Če odštejemo tretjo enačbo iz sistema (2.1.2) od prve in druge, dobimo: aix + by + Ciz + = O «2 X + b2 y + C2 z + d2 = O Sistem enačb (2.1.3) zapišemo matrično: (2.1.1) (2.1.2) (2.1.3) x a1 b1" -1 - d1 - C1z 'a ß~ - d1 - C1z y a2 b2 _ - d2 - C2 z J S_ - d2 - C2 z (2.1.4) kjer so: ai = x3 -xi, bj = y3 -yj, cj = z3 -zj , dj = I 2 2 2 r,2\ { 2 2 2 n2 - Ix2 + y3 + z3 - P3 j + X2 + y2 + z2 - Pj 2 j = 1,2. (2.1.5) Od tod izrazimo parametrično enačbo za položaj: X = x0 +Ä-1, y = y0 1, z = t, kjer so: x0 =-«■ d1 -ß- d2, y0 =-/■ d1 -S- d2 Ä =-a ■ C1 -ß- C2, ^ = C1 -S- C2. (2.1.6) (2.1.7) Izraze (2.1.5) do (2.1.7) vključimo v sistem enačb (2.1.2) in jih preuredimo, da dobimo kvadratno enačbo (neznanka t): g1 ■ t2 + 2 • g2 • t + g3 = O (2.1.8) kjer so: g1 =12 + +1, g2 = -[(x:1 -x0)^ + (y1 -yo)ß + zl], g3 =(x1 -x0)2 + (y1 -yo)2 + z2 -P12 (2.1.9) Dobimo 3 enačbe s tremi neznankami, kjer najprej določimo t, potem pa še X in ju ter naprej neznani položaj sprejemnika X z izrazi (2.1.6). 0 i? 1 -V o IT 2.2 Direktna rešitev - z apriori znano elipsoidno višino Iz sistema treh enačb psevdorazdalj (1.2.4) določamo štiri neznanke ob predpostavki, da se višina opazovališča h ne spreminja oziroma je vnaprej znana (angl. altitude hold). To zapišemo s pomočjo kvadratne forme (Phatak et al., 1999): X ■ r ■ rT ■ xT = 1 Pogoj o znani višini vključuje diagonalna matrika: 1 1 1 r = diag aWGS-84 + h aWGS-84 + h WGS-84 + h (2.2.1) (2.2.2) K !?: o ;s pri čemer sta awGs-84 in bWGS _84 polosi elipsoida WGS-84. Problem rešimo tako, da neznani vektor položaja sprejemnika X zapišemo kot funkcijo količine b, naprej pa s formiranjem polinoma četrte stopnje za b s substitucijo za X določimo b. Ko poznamo b, določimo še X. Modifikacija nekaterih metod (npr. Krause, 1987) omogoča direktno rešitev problema in direktno računanje v koordinatnem sistemu ECEF (angl. Earth Centered Earth Fixed). Pri izpeljavi modificirane rešitve Krause izhajamo iz treh enačb oblike (1.2.4). Definirajmo dva vektorja V1, V 2 ter matriko O: = xT - x^, V2 = xT - xj, n = V1 v 2 (2x3) Naprej določimo r1 = • n + u • l1 , kjer je: ri = xT - xT za i = 1,2,3, (2.2.3) (2.2.4) pri čemer je u enotski vektor, pravokoten na ravnino, ki jo določajo trije sateliti. Rešitev za podajajo izrazi: ^T = \ Mr1 rešitev za l1 pa: ,2 S11 Soo M1 = n • nT, Sii = v i • vT, Ji -P12, i =1,2, (2.2.5) i12 =P12 -1^1 ■ ^2. (2.2.6) Ko določimo in l1, določimo r1 (2.2.4), naprej pa velja:x = x1 + r1. V primeru psevdorazdalj imamo dodatno neznanko b, zato enačbo za l1 (2.2.6) zapišemo drugače: 2 i12 =(/1 - b)2 -1^1 ■ n-(^1 - )■ n kjer sta: (2.2.7) to iT = T Mf1 I-Ji , JJi = _PI2, i = 1,2. Velja: Ji - Ji = 2 • b ■ K in i; T - ^ T ^ b. M r-1 Ki K 2 , kjer je: Ki = Pi- Pi Zgornje izraze vključimo v (2.2.7), da določimo l1 kot funkcijo b: 't • 1 - bf - gl - b ■« (2.2.8) (2.2.9) (2.2.10) (2.2.11) kjer sta: gl = Ii ■ n, « = [^1 KtiMf1«. l1 izračunamo iz kvadratne enačbe: 'i2 = a1 ■ b2 + a2 ■ b + a3, (2.2.12) a1 = 1 - M ■ a2 = 2 ■(gg1 ■ M T -1~1 a3 = .P12 - g1 ■ g1 . V osnovi so metodo izpeljali tako, da so b določili iz enačb štirih psevdorazdalj (rešujemo kvadratno enačbo za b). V danem primeru pa do kvadratne enačbe pridemo tako, da X določimo kot funkcijo b in l1 in naprej vključimo v izraz (2.2.2). Iskani vektor položaja sprejemnika zapišemo kot: X = X1 + • n + u • l1 = X1 +• n- )• n + u • l1 = X1 + g1 - m • b + u • l1 • X (2.2.13) Upoštevamo še: X-r = a + ß- b + 'y-/1 a = (x1+ gl)-r, ß = -®■ r, y = u■ r. (2.2.14) Zgornji izraz vstavimo v pogoj (2.2.1) ob upoštevanju kvadratne enačbe za l1 (2.2.13), da dobimo: A + ^-1} + C-b2 = l1 -{d + E ■ b),kjer so: (2.2.15) A = a aT +y-yT ■ a3 -1, B =y yT ■ a2 + 2-a-ßT C = ß ßT +y yT ■ a1, D = -2■a yT, E = -2■ß yT (2.2.16) izračunamo iz polinoma četrte stopnje, ki ga dobimo tako, da kvadriramo obe strani enačbe (prvi izraz 2.2.15) in naredimo ponovno substitucijo za l1 (2.2.12) ■I Sf to ^^ _^ (c2 - a^E2 + {lBC - la^DE - a-^E2 )b3 + _ + (b2 + 2AC - a^D^ - la2DE - a^E"^-)b2 + (lAB - a2D^ - la3DE)b + (a2 - a^D"^ )= 0 (2.2.16) Od tu naprej rešujemo enačbo (2.2.16) analitično - določamo ničle polinoma četrte stopnje. Ko določimo pravo vrednost za b, preko /i izračunamo vektor neznanega položaja sprejemnika X. 55 .t Pravo vrednost izmed štirih izračunanih za b in pravo vrednost izmed dveh izračunanih za /1 (torej osem vrednosti) dobimo po premisleku. Predznak za /1 izberemo tako, da zadostimo (2.2.15), ostanejo še 4 možne rešitve za X. Izmed teh sta navadno dve, ki ne zadostita enačbi (1.2.4), od preostalih možnih rešitev pa kot verjetno rešitev izberemo tisto, ki zadosti geometrijskemu pogoju in zastavljenemu pogoju pogreška ur sprejemnika. Geometrijski pogoj se nanaša na predpostavko, da morajo biti vsi koti med sprejemnikom in sateliti ostri - drugače povedano: produkt (x,- - x) • X mora biti pozitiven. Pri uporabi drugega pogoja moramo vedeti, ali situacijo določitve položaja sprejemnika iz treh psevdorazdalj rešujemo na začetku ali v sredini niza zaporednih opazovanj. V drugem primeru si pomagamo z izračunanimi parametri ur iz prejšnjega koraka in jih ekstrapoliramo za izbrani trenutek. Če situacijo s tremi opazovanji rešujemo na začetku opazovanj, so parametri ur v znanem časovnem intervalu.4 Dodatno zanemarimo še nesoglasje trenutka oddaje glede na GPS-čas in upoštevamo, da so pogreški sprejemnikove ure v definicijskem območju [-10 ms, 10 ms], kar ustreza razdaljam v območju [-3000000 m, 3000000 m]. Velja pa vedeti, da opisana metoda ne deluje vedno; takrat moramo počakati na opazovanja četrtega satelita. 3 MATEMATIČNI MODELI DOLOČITVE ABSOLUTNEGA POLOŽAJA IZ ŠTIRIH ALI VEČ PSEVDORAZDALJ Znanih je več matematičnih modelov določitve absolutnega položaja iz štirih ali več psevdorazdalj. Glede na način rešitve jih lahko razdelimo v direktne ali iterativne, slednje pa glede na to, ali tekom obdelave izvedemo postopek linearizacije ali ne. Pri linearizaciji je dobro poznati približne vrednosti neznank, da postopek izračuna poteka hitreje. Hkrati je delovanje metod odvisno tudi od geometrijske razporeditve satelitov in opazovališča (faktor GDOP). Najbolj elegantno rešitev problema (tudi časovno najmanj potratno) za situacijo z natanko štirimi opazovanji, ki so jo kasneje izpeljali tudi za večje število opazovanj, je podal Bancroft. Na podobni osnovi so direktno rešitev predstavili tudi drugi avtorji (Crocetto et al., 1998). Slednji so predstavili tudi iterativno rešitev brez linearizacije enačb psevdorazdalj. V smislu geometričnega razumevanja situacije je zelo zanimiva analitična rešitev, ki jo je podal Kleusberg (1994). Vsaka izmed metod ima svoje prednosti in slabosti, ki jih podrobneje opisujemo v nadaljevanju. V smislu razumevanja problema in rešitev predstavljamo naštete metode z vmesnimi izpeljavami in komentarji. ^ GPS-sprejemniki namreč delujejo tako, da za začetni trenutekprivzamejo trenutek oddaje (po prvem sprejemu signala) in ga dodatno povprečijo za čas potovanja signala, kar znaša približno 70 ms. 3.1 Algorite m tJi^n croft Leta 1985 je Bancroft predstavil direktno rešitev določitve absolutnega položaja za primer enostranskih opazovanj dolžin. Temelji na matričnem reševanju sistema linearnih enačb in nadaljnji rešitvi skalarne kvadratne enačbe. Ta sicer vodi do dveh rešitev, vendar je ena izmed njih povsem nesmiselna. Algoritem lahko uporabimo tudi na drugem tipu opazovanj, znana pa je tudi uporaba algoritma na primeru kombinacije eno- in dvostranskih opazovanj, kjer določamo ničle polinoma četrte stopnje (Anderson in Tran, 2003). Izhajamo iz izhodišča, da imamo v danem trenutku na voljo položaj vsaj štirih satelitov v geocentričnem koordinatnem sistemu (/-ti satelit: x, y, z) ter smo opazovani psevdorazdalji Pj prišteli vpliv popravka satelitovih ur c ■ dti, v drugi iteraciji pa odšteli tudi popravek troposferske refrakcije (Tropo)5 in upoštevali aberacijo Zemlje. V vektor s,- uvrstimo znane količine vezane na /-ti satelit, v vektor u pa neznanke (položaj opazovališča: y, z ter pogrešek sprejemnikove ure dT, ki ga zaradi numerične stabilnosti modela pomnožimo s hitrostjo svetlobe c): s,- = Xj y, z, iT , = [x y z Za satelit / ter za določene vplive popravljeno psevdorazdaljo ter za neznane količine iz vektorja u zapišemo funkcionalni model za izbrani trenutek (Yang in Chen, 2001): Pj - b = - X)2 + - y)2 ^{z^ - z)2 . Funkcijo kvadriramo: - 2Pjb + b2 = ixj - x)2 + (y, - y)2 + (z, - z)2 2 ^ 2 2^ 22^ 2 = xi - 2xix + x + yi - 2yjy + y + z,- - 2ziz + z in preuredimo: xf + yf + zf -)- 2(x,x + y,y + z,z - I'ib)+[x^ + y2 + z2 - b2) = 0. Nadalje enačbo (3.1.3) preuredimo tako, da si pomagamo z Lorentzovim produktom6: (3.1.1) (3.1.2) (3.1.3) {si,s) -{si, ^ + \{u,u = 0 (3.1.4) 5 Troposfersko refrakcijo model/ramo v odvisnost/ od višinskega kota satelita /n meteoroloških pogojev na opazoval/šču, k/ so odv/^n/ od v/š/ne opazovališča h. V prvi iteraciji predpostavljamo, da se nahajamo v središču Zemlje (x0, y0, Z0) = (O, O, 0), zato nimamo dovolj podatkov za določitev parametrov, s katerimi določimo vpliv troposferske refrakcije na izmerjeno dolžino. 6 Lorentzov produktje definiran kot: ^u, V = uTMv - u^Vi + u^v^ + U3V3 — U4V4, kjer sta: u,v e Ä4 in matrika M = 13x3 0 0 -1 , M"' = M = Mt . Velja tudi: (Mu, M^ = (u, v) . i? I o IT Izhajamo iz oznak: X1 J1 Z1 P1 ( s1,s0 "1' x2 y 2 z 2 P~2 , a = 1. 2 ( 1 B = X3 y 3 Z 3 ( 83. , e = 1 X4 y 4 Z 4 P4 { s4,s^ 1 A = U'^ • (3.1.5) I P V primeru, da imamo na voljo n psevdorazdalj, je matrika B dimenzije n x 4, vektorja a, e sta velikosti n x 1, medtem ko je A skalar. Sistem n enačb z oznakami (3.1.5) preuredimo v obliko (3.1.4): a - BMu + Ae = 0 (3.1.6) Če imamo na voljo opazovanja od štirih satelitov, lahko rešimo enačbo (3.1.6): u = M ■ B"1 -(Ae + a) (3.1.7) Ker količina A tudi vsebuje vektor u, enačbo (3.1.7) vstavimo v enačbo (3.1.6), uporabimo zvezo (Mu,M^ = (u,\) in enačbo kvadriramo, da dobimo: (3.1.8) Obstajata dve rešitvi kvadratne enačbe (3.1.8), ena izmed teh je pravilna (primerljiva velikost s polmerom Zemlje). Če imamo na voljo več kot štiri psevdorazdalje, enačbo (3.1.6) v skladu z metodo najmanjših kvadratov preoblikujemo v niz normalnih enačb tako, da levo stran (3.1.6) pomnožimo z BT: BTa - BtBMu + B T Ae = 0. (3.1.9) Nadalje postopamo enako kot prej, da dobimo kvadratno enačbo, kjer je: ^B*e, B*^A2 + B*^ - B*a, B*^ = 0 (3.1.10) in izmed dveh rešitev kot pravilno vzamemo bolj verjetno. Neznani vektor u naprej določimo s pomočjo izraza (3.1.7). •iS Sfj 3.2 Direktna rešitev - Crocetto et al. Izpeljava temelji na podobnem izhodišču kot algoritem Bancroft, tj. na kvadriranju enačbe psevdorazdalje (3.1.1). S postopkom eliminacije neznanke iz sistema štirih enačb psevdorazdalj dobimo: (x4 - xi)■ x + (j4 - y + (24 - Zj)-z - -1\)■ b + fi = 0 (x4 - X2 )■ x + (>^4 - y2 )' y + (24 - Z2 )' z-(^4 - Pi h b + fi = 0 (x4 - X3)■ x + (y4 - y3)■ y + (24 - Z3)■ z - (P4 - /3)■b + f, = 0, kjer so: fi y2 + zi -p^)-(x| + y| + zj -Pj i = 1,2,3. (3.2.1) (3.2.2) Sistem enačb lahko zapišemo kot: a1x + b1y + c1d +d1b + f1 = 0 aj x + bj y + Cjd +d jb + fi = 0 a3x + b3 y + c3d +d 3b + f3 = 0 (3.2.3) Neznani položaj antene GPS-sprejemnika izračunamo iz zgornjega sistema enačb, zapisanega matrično: x a1 b1 C1 -1 d1 ■ b - /1 " «1 71 --d1 ■ b - /l " xq +1 ■ b y = a2 b2 C2 - d2 ■ b - /2 = «2 ß2 72 - d2 ■ b - f = y0 + m ■ b z _a3 b3 c3 _ - d2 ■ b - /3 «3 ß3 73 _ d2 ■ b - /3 _ zq + n ■ b _ (3.2.4) kjer so: x0 = y0 = z0 = -«1 ■ f1 -ß\ ■ f2-71 ■ f3 I -«2 ■ f1 -ßi ■ fl-72 ■ f3 m -«3 ■ f1 -ß3 ■ f2-73 ■ f3 n = -«1 ■ d1 -ßi- d2-71 ■ d3 = - «2 ■ d1 - ß2 ■ d2 -72 ■ d3 = -«3 ■ d1 -ß3 ■ d2-73 ■ d3 (3.2.5) Pogreške sprejemnikove ure (b = c ■ dT), določimo tako, da v sistem enačb (3.2.4) vstavimo v (3.1.3). Z oznakami: 72 2 2 1 e1 = l + m + n -1 e2 = (x0 -x1)-l + (y0 -y1 )■m + (z0 -z1 )■ n + 1x0 - x1)2 +(y0 - y1)2 +(z0 - z1 )2 - p2 (3.2.6) rešujemo kvadratno enačbo: e1 ■ b2 + 2 ■ e2 ■ b + e3 = 0. (3.2.7) Dobimo dve rešitvi za b, ki ju naprej uporabimo za določitev neznanega položaja sprejemnika (3.2.4). Kot pravo rešitev (položaja sprejemnika) izberemo tisto, katere dolžina vektorja bolj ustreza polmeru Zemlje. Opisana rešitev je zanimiva zato, ker najbolj nazorno opisuje situacijo, da je rešitev kvadratne ■o o IT e = enačbe (3.2.7) odvisna od geometrije sateliti-sprejemnik (kar vključujejo količine a, ,ßi, , fj) in od pogreškov izmerjenih psevdorazdalj (vključujejo dj in fi). Metodo izračuna položaja opazovališča lahko uporabimo tudi za situacijo, ko imamo na voljo več kot štiri psevdorazdalje. V tem primeru izračunamo vse kombinacije rešitev in iz teh izračunamo povprečje. Število kombinacij štirih satelitov pri večjem številu satelitov na obzorju n strmo naraste: ko imamo na voljo opazovanja petih satelitov, je kombinacij 5, pri sedmih 35, pri devetih pa že 126.7 Obstaja pa tudi hitrejša rešitev: približni položaj sprejemnika določimo le za eno kombinacijo štirih satelitov, drugi korak naredimo s postopkom linearizacije enačb psevdorazdalj (podpoglavje 3.4). 3.3 Rešitev z vektorsko algebro - algoritem Kleusberg Kleusbergova rešitev temelji na vektorski algebri, tj. rešitvi geometrije tridimenzionalnega hiperboličnega preseka, in na predpostavki, da imamo na voljo opazovanja štirih satelitov. Vzemimo za referenčni satelit (Sat 0) enega izmed satelitov na obzorju, npr. satelit z največjim višinskim kotom. Naprej izračunamo razdalje b. med referenčnim in drugimi sateliti in enotske vektorje e, med referenčnim satelitom in drugimi sateliti (Sat i). Tako dobljene količine opisujejo dejansko geometrijo razporeditve satelitov glede na naše opazovališče. Neznani položaj antene sprejemnika X opišemo s pomočjo neznanega enotskega vektorja e, tj. normiranega vektorja med referenčnim satelitom in sprejemnikom. Z df označimo razliko razdalje satelit-sprejemnik od razdalje referenčni satelit-sprejemnik (Kleusberg, 1999): di - x )2 - y )2 ^{z, - z )2 ^(xo - x)2 +(yo - y )2 +(zo - z)2 . (3.3.1) Izračunamo razdalje od referenčnega do ostalih satelitov: bi -xo)2 +{yi -yo)2 +{zi -zo)2, i = 1,2,3 od tod pa enotske vektorje: (3.3.2) xi - xo yi - yo zi - zo bi bi bi (3.3.3) •o tig 55 Naprej rešujemo geometrijo trikotnikov Sat i - Sat 0 - P, pri čemer iz skalarnega produkta med enotskima vektorjema e in e i določimo vmesni kot med vektorjema: s^2 = bf + s2 - 2bjS0 (e ■ ei) (3.3.4) 7 Izhajamo iz enačbe za izračun št^-vila kombinacij: ' 4! {n - 4)! ,, kjer je n št^-vilo vseh satelit^^, ki j^h uporabimo. T e = n Slika 3: Geometrija tridimenzionalnega hiperboličnega preseka. Enačbo (3.3.1) lahko zapišemo kot: s, = d^ + s0, (s0 predstavlja razdaljo med referenčnim satelitom in sprejemnikom, s, pa razdalje med satelitom i in sprejemnikom), da v vsakem od trikotnikov dobimo povezavo med razlikami psevdorazdalj in razdaljami med sateliti in sprejemnikom: sf = df + S02 + 2diS0 (3.3.5) Enačbi (3.3.4) in (3.3.5) izenačimo, da dobimo: , 0 di + bi (e ■ e,-), i = 1,2,3. (3.3.6) V sistemu treh enačb (3.3.6) so tri neznanke: dve neodvisni komponenti enotskega vektorja e in razdalja s0. Neznano razdaljo iz sistema enačb (3.3.6) eliminiramo: bi - dj bi+1 - di+1 i = 1,2 . di + bi (e ■ ei) di+1 + bi+1 (e ■ ei+1 K Na osnovi distributivnostnega zakona vektorske algebre enačbi (3.3.7) preuredimo: i = 1,2 (3.3.7) bi bi 'i+1 b2 - d} bi+1 - di+1 -ei+1 e = di i+1 di bh - df,1 b2 - d} in v krajši obliki zapišemo: F1 ■ e =U,, F2 ■ e =U2. (3.3.8) (3.3.9) Vektor F1 leži v ravnini, ki jo določajo sateliti 1, 2 in 3, vektor F2 pa v ravnini, ki jo določajo sateliti 1, 3 in 4. Oba vektorja in skalarja U^ in U2 lahko izračunamo iz znanih koordinat satelitov in iz razlik dolžin med njimi. Vendar želimo geometrično rešitev predstaviti na osnovi vektorske algebre. Enačbi (3.3.9) sta skalarna produkta vektorjev F1 oziroma F2 in enotskega vektorja e. V T o ' -cž ^^ H ^ n •i? o ^^ .!3 splošnem ima problem dve rešitvi: prva je nad in druga pod ravnino, ki jo določata vektorja F1 in F2, le v primeru, ko sta vektorja vzporedna, rešitev ne obstaja. Algebraično rešitev enačbe (3.3.9) lahko pridobimo z vektorskim produktom: ex(f1 xf2) = F1 -(e■ F2)-F2 -(e■ F1 ) = U2 ■ F1 -U1 ■ F2 (3.3.10) Z okrajšavama: G = F1 X F2, H = U2F1 - U;F2 (3.3.11) enačbo (3.3.10) zapišemo v bolj enostavni obliki: ex G = H. (3.3.12) Obe strani enačbe (3.3.13) množimo z G, za levo stran pa uporabimo razčlenitev vektorskega produkta: e-(g • G)" G-(G • e) = G x H. (3.3.13) Skalarni produkt drugega člena leve strani enačbe (3.3.13) zapišemo kot produkt dolžin vektorjev in kosinusa vmesnega kota 5: G ■ e = ^ G ■ G ■ cos((y). Kot 5 se nahaja tudi v izrazu za izračun dolžine vektorja H iz enačbe (3.3.13): 4h ■ H ^(e X G )-(e x G) = V G ■ G ■ sin(