Doktorska disertacija UPORABA KALMANOVEGA FILTRA PRI KINEMATIČNIH GEODETSKIH MERITVAH Sonja Gamse, univ.dipl.inž.geod. Julij, 2010 IZJAVA O AVTORSTVU Podpisana Sonja Gamse, univ. dipl. inž. geod., izjavljam, da sem avtor doktorske disertacije z naslovom: UPORABA KALMANOVEGA FILTRA PRI KINEMATIČNIH GEODETSKIH MERITVAH. Ljubljana. (podpis avtorja) STRAN ZA POPRAVKE Stran z napako Vrstica z napako Namesto Naj bo Ne moreš obstajati, če ne delaš, ne moreš obstajati, če ne misliš. (Dr. France Bučar) Posvečeno mami ZAHVALA Ob koncu enega dela prehojene poti - nočem reči končanega dela, saj že na tem mestu vidim nove izzive - je tudi priložnost, da se ozrem nazaj v čas in naredim analizo prehojene poti. Kaj mi je prinesla, kaj sem si znala vzeti, kateri so bili tisti ključni trenutki, vrhovi in padci, kdo me je spremljal na tej poti ali se ji priključil ... Težko bi bilo vse občutke in misli strniti na eno stran, zato le kratka, a iskrena hvala vsem, ki ste mi na tej poti tako ali drugače pomagali. Hvala mentorju izr. prof. dr. Dušanu Kogoju za dano priložnost, zaupanje in strokovne napotke. Ich bedanke prof. dr. Thomas A. Wunderlich fur das feste Vertrauen, dass alle Tiire am Lehrstuhl fur Geodasie an der Technischen Universitat Miinchen mir immer offen waren. Ich bedanke Ihnen und Ihren Mitarbeiter - dr. Werner W. Stempfhuber, dr. Karl Foppe und bedonders dr. Peter Wasmeier fur die weitergegebene Kenntnisse und die Hinweise bei meiner Arbeit. Ich bedanke mich fur die schone und fruchtbare Zeit in Miinchen. Hvala Jaki in podjetju Geograd za izposojo instrumenta na začetku moje poti. Hvala podjetju Geoservis in Gregorju za vso posredovano znanje. Hvala Skladu dr. Otta Likarja, dipl. inž., in Karle Likar, prof. dr. Mihi Boltežarju, in skladu Deutscher Akademischer Austausch Dienst za zaupanje in finančno podporo, brez katere bi bila izbrana pot izjemno težko uresničljiva. Hvala prof. dr. Goranu Turku in izr. prof. dr. Mihaelu Permanu za temeljit pregled doktorske disertacije in plodne konzultacije, ki so delu dale dokončno obliko in vsebino. Hvala tistim, ki ste me znali poslušati. Z vami mi je bilo lažje, ob vas nisem pozabila še na druge plati življenja. Hvala Bojan, Gašper, Iztok, Dejan, Barbi, Primož, Tina, Robert, Nataša, Primož, Oskar, Dejan, Lucija. Ich bedanke mich Walkiria und Frank fur die Hilfe und unvergefiliche Zeit und Beziehung in Miinchen. Hvala Ivi, Vesni, Speli, Maji, Maidi; hvala vam za čudež, ki ste ga s svojim strokovnim delom in osebnim pristopom naredile pri Veroniki in Aleksandru. Hvala vam, da sem lahko bila vsak trenutek sproščena pri svojem delu, ko sta bila Veronika in Alekdander v vašem varstvu. Car trdega dela je tudi v tem, da ti v življenje lahko prinese posebne ljudi, ki ostanejo s tabo tudi, ko delo dobi zaključeno obliko. Hvala staršem za svobodo. Od kar se zavedam sebe, sta mi pustila svet okoli mene v lastno presojanje in samostojno izbiro poti. Hvala Robi za ... premalo imam prostora; nama bi lahko posvetila knjigo. BIBLIOGRAFSKO - DOKUMENTACIJSKA STRAN IN IZVLEČEK UDK: Avtor: Mentor: Komentor: Naslov: Obseg in oprema: Ključne besede: Izvleček 519.22:528(043.3) Sonja Gamse, univ. dipl. inž. geod. izr. prof. dr. Dušan Kogoj univ. prof. dr. - ing. habil. Thomas A. Wunderlich Uporaba Kalmanovega filtra pri kinematičnih geodetskih meritvah. 96 str., 5 pregl., 59 si., 81 en. diskretni model Wienerjevega procesa s pospeškom, elektronski tahimeter, izboljšava meritev, Kalmanov filter, kinematični proces, stanje sistema Zajem deformacij in premikov, tako v smislu merskega zajema kot analize rezultatov, danes zahteva interdisciplinarno delo strokovnjakov iz različnih področij. Prednost geodetskih postopkov zajema deformacij in premikov je zajem le-teh v prostoru in času. V primeru, da obravnavamo deformacije in premike kot kinematične, lahko za obdelavo kinematičnih terestričnih opazovanj uporabimo Kalmanov filter. V okviru doktorske disertacije so bile pridobljene konkretne meritve kinematičnega procesa. Meritve so bile izvedene z enim izmed najsodobnejših elektronskih tahimetrov z možnostjo samodejnega sledenja gibajočega reflektorja. Kinematičen proces je bil izveden na osnovi gibanja vozička z reflektorjem na tirnici, podani v lokalnem koordinatnem sistemu. Za zajem meritev je bil pripravljen program za vzpostavitev komunikacije z instrumentom, samodejno sledenje reflektorja in pridobivanje ter shranjevanje meritev. Za obdelavo meritev je bil uporabljen diskretni model Wienerjevega procesa s pospeškom, s predhodno uporabo zakona o prenosu varianc in kovarianc. Poudarek doktorske disertacije je na statistični oceni razvitega modela, ki nam omogoča odkrivanje prisotnosti grobih pogreškov v meritvah in podaja ustreznost izbire vhodnih količin modela. Izvirnost doktorske disertacije predstavlja neodvisna kontrola zanesljivosti matematično stohastičnega modela, ki temelji na neodvisni referenčni tirnici. Prednost referenčne tirnice je možnost istočasnega vrednotenja modela in zmogljivosti instrumenta. Prispevek dela podaja ključna izhodišča, na katera je potrebno biti pozoren pri kinematičnih geodetskih terestričnih postopkih. BIBLIOGRAPHIC - DOCUMENTALISTIC INFORMATION UDC: Author: Supervisor: Supervisor: Title: 519.22:528(043.3) Sonja Gamse, univ. dipl. ing. geod. Assoc. Prof. Dr. Dušan Kogoj Univ. Prof. Dr. - Ing. Habil. Thomas A. Wunderlich The Use of Kalman Filer for Kinematic Geodetic Observations Notes: Key words: 96 p., 5 tab., 59 fig., 81 eq. discrete Wiener process acceleration model, electronic tacheometer, measurement innovation, Kalman filter, kinematic process, system state Abstract Surveying of deformations and movements, which means measuring of processes and analyzing of observations, today requires interdisciplinary work of experts from different disciplines. The advantage of geodetic procedures in surveying of deformations and movements is their definition in space and time. In the case of deformations and movements regarded as kinematic, Kalman filter can be used for the processing of kinematic terrestrial observations. For the PhD thesis the concrete measurements of the kinematic process have been performed. Measurements were made with one of the most modern electronic tacheometers capable of automatic tracking of a moving reflector. Kinematic process was carried out based on the motion of the trolley with a reflector moving along a straight trajectory, given in a local coordinate system. To capture measurements a software was prepared to establish communication with the instrument, to automatically track the reflector and to acquire and store the measurements. For the processing of measurements the discrete Wiener process acceleration model was used, with prior application of the law of propagation of variances and covariances. The dissertation focuses on the statistical evaluation of the developed model, which allows us to detect the presence of gross errors in measurements and provides adequacy of the given input values. Originality of the thesis is represented by independent control of mathematical stochastic model reliability based on an independent reference trajectory. The advantage of the reference trajectory is the possibility of simultaneous evaluation of the model and capabilities of the instrument. The contribution of the work is in giving several demands that need extra attention in the case of kinematic geodetic terrestrial processes. Kazalo vsebine 1 uvod i 1.1 Razvrstitev deformacijskih modelov............................................1 1.2 Kalmanov filter..................................................................2 1.3 Opis problema....................................................................4 1.4 Vsebina doktorske disertacije....................................................6 2 kalmanov filter 9 2.1 Linearni Kalmanov filter........................................................11 2.1.1 Kalmanova matrika izboljšave Kk......................................14 2.2 Diskretizacija: Prehod iz zveznega problema na diskretni linearni problem . 16 2.2.1 Definicija kovariančne matrike šuma sistema Qk......................18 2.3 Razširitve Kalmanovega filtra..................................................19 2.3.1 Ocena stanja sistema v primeru nelinearnih sistemov................19 2.3.2 Postopek linearizacije: Razvoj v Taylorjevo vrsto ....................20 3 statistična analiza časovno odvisnih procesov 21 3.1 Ocena stanja sistema v diskretnem linearnem časovno odvisnem modelu . . 23 3.2 Preizkušanje konsistenčnosti postopka za ocenjevanje stanja sistema .... 23 3.2.1 Testna statistika v domeni meritev....................................23 3.2.2 Testna statistika v domeni vektorja stanja sistema....................25 3.3 Zmožnost nadzora in opazovanosti sistema....................................25 3.3.1 Zmožnost nadzora sistema..............................................25 3.3.2 Zmožnost opazovanosti sistema........................................26 4 kalmanov model meritev in sistema (računski primer) 27 4.1 Opis problema....................................................................27 4.2 Vzpostavitev praktičnega primera..............................................27 4.2.1 Referenčni okvir..........................................................27 4.2.2 Instrumentarij............................................................30 4.3 Funkcionalni in stohastični model: Diskretni model Wiener-jevega procesa s pospeškom......................................................................37 5 vrednotenje diskretnega modela wienerjevega procesa s pospeškom 45 5.1 Računske zahteve................................................................46 5.2 Kazalci notranjega zaupanja....................................................46 5.2.1 Zmožnost opazovanosti DMWPP......................................46 5.2.2 Zmožnost nadzora DMWPP............................................46 5.2.3 Lastnosti matrike P ....................................................46 5.2.4 Lastnosti matrike K....................................................51 5.3 Statistični testi..................................................................58 5.3.1 Položaj reflektorja ......................................................58 5.3.2 Skladnost DMWPP v domeni meritev ................................66 5.3.3 Skladnost DMWPP v domeni stanja sistema..........................69 5.4 Referenčna tirnica kot parameter ocenjevanja ................................72 5.5 Povzetek rezultatov..............................................................76 6 ZAKLJUČKI 81 6.1 Povzetek prispevka dela ........................................................82 6.2 Nadaljnje raziskave..............................................................84 6.3 Zaključek ........................................................................86 7 POVZETEK 89 8 SUMMARY 91 VIRI 93 PRILOGE 97 A Rezultati programa CAPLAN 99 B Razlike v meritvah 103 Kazalo preglednic 1.1 Razvrstitev in lastnosti deformacijskih modelov [Welsch, Heunecke (2001)] 2 4.1 Koordinate stojišča instrumenta s standardnimi deviacijami, Priloga A . . 29 4.2 Tehnični podatki instrumenta TCRA1201 Leica Geosystems [Tech. Data] . 32 5.1 Vsota absolutnih razlik med napovedanimi oz. filtriranimi in merjenimi vrednostmi........................................................................63 5.2 Odstotek filtriranih vrednosti znotraj intervala \a in intervala 2a ..........65 Kazalo slik 1.1 Diagram sistem-meritve-Kalmanov filter, povzeto po [Gelb (1974)]..... 3 2.1 Zanka KF, [Brown, Hwang (1992)] ...................... 10 2.2 Zanka KF.................................... 13 4.1 Geodetski laboratorij Katedre za geodezijo, TUM, z referenčno tirnico ... 28 4.2 Reflektor GRZ4 Leica Geosystems 360° na vozičku.............. 29 4.3 Položaj stojišča instrumenta (StandPunkt_S) in danih točk v referenčnem (y, x) sistemu laboratorija........................... 30 4.4 Komunikacijski protokoli Leica Geosystems.................. 33 4.5 Pregled aplikacije odjemalec - strežnik [GeoCOM Reference Manual V 1.20] 33 4.6 Razvoj geodetskih instrumentov........................ 35 4.7 Celotni računski postopek KF, z zakonom o prenosu varianc in kovarianc . 43 5.1 Konvergenca sledi a posteriori kovariančne matrike stanja sistema P+ ... 47 5.2 Konvergenca lastnih vrednosti 1,2,3 a posteriori kovariančne matrike P+ . 48 5.3 Konvergenca lastnih vrednosti 4,5,6 a posteriori kovariančne matrike P+ . 48 5.4 Konvergenca lastnih vrednosti 7,8,9 a posteriori kovariančne matrike P+ . 49 5.5 Pogojno število k(P+), brez grobih pogreškov ................ 50 5.6 Pogojno število k(P+), z grobimi pogreški.................. 50 5.7 Pogojno število k(P+) za različne vrednosti aw................ 51 5.8 Elementi Kalmanove matrike izboljšave K; rdeče - elementi matrike K za komponento x merskih popravkov, zeleno - eleme nt i matrike K za komponento y merskih popravkov, modro - eleme nt i matrike K za kompone nto z merskih popravkov: Položaj.......................... 52 KK xK nento y merskih popravkov, modro - eleme nt i matrike K za kompone nto z merskih popravkov: Hitrost .......................... 53 KK xK nento y merskih popravkov, modro - eleme nt i matrike K za kompone nto z merskih popravkov: Pospešek......................... 54 5.11 Elementi Kalmanove matrike izboljšave K; rdeče - elementi matrike K za komponento x merskih popravkov, zeleno - elementi matrike K za komponento y merskih popravkov, modro - elementi matrike K za kompone nto z merskih popravkov: Položaj, aw = 0.1 . . 55 KK x merskih popravkov, zeleno - elementi matrike K za komponento y merskih popravkov, modro - elementi matrike K za komponento z merskih popravkov: Hitrost, aw = 0.1 . . 55 KK £ merskih popravkov, zeleno - elementi matrike K za komponento y merskih popravkov, modro - elementi matrike K za komponento z merskih popravkov: Pospešek, aw = 0.1 . 56 KK £ merskih popravkov, zeleno - elementi matrike K za komponento y merskih popravkov, modro - elementi matrike K za kompone nto z merskih popravkov: Položaj, aw = 0.000001 56 5.15 Elementi Kalmanove matrike izboljšave K; rdeče - elementi matrike K za komponento x merskih popravkov, zeleno - elementi matrike K za komponento y merskih popravkov, modro - elementi matrike K za komponento 2 merskih popravkov: Hitrost, aw = 0.000001 57 KK x merskih popravkov, zeleno - elementi matrike K za komponento y merskih popravkov, modro - elementi matrike K za komponento 2 merskih popravkov: Pospešek, aw = 0.000001 57 5.17 Položaj x reflektorja: rdeče - x položaj izračunan iz meritev, zeleno - napovedane vrednosti položajev x, modro - filtrirane vrednosti položajev x . 59 xx xx 5.19 Položaj y reflektorja: rdeče - y položaj izračunan iz meritev, zeleno - napovedane vrednosti položajev y, modro - filtrirane vrednosti položajev y . . 60 yy yy 5.21 Položaj z reflektorja: rdeče - z položaj izračunan iz meritev, zeleno - napovedane vrednosti položajev z, modro - filtrirane vrednosti položajev z . . 61 yy yy 5.23 Horizontalna ravnina............................................................62 5.24 Standardne deviacije položajnih komponent ..................................62 5.25 Napaka položaja x (polna rdeča črta), interval 1a (zelena polna črta) in interval 2a (prekinjena zelena črta), brez prisotnosti grobih pogreškov, aw = 0.1 63 5.26 Napaka položaja y (polna rdeča črta), interval 1a (zelena polna črta) in interval 2a (prekinjena zelena črta), brez prisotnosti grobih pogreškov, aw = 0.1 64 5.27 Napaka položaja z (polna rdeča črta), interval 1a (zelena polna črta) in interval 2a (prekinjena zelena črta), brez prisotnosti grobih pogreškov, aw = 0.1 64 5.28 Normirana kvadratna izboljšava................................................66 5.29 Izboljšave meritev KF ..........................................................67 x y z 5.33 Normirana napaka stanja sistema....................... 70 5.34 Histogram popravkov komponent stanja sistema - položajne komponente . 70 5.35 Histogram popravkov komponent stanja sistema - komponente hitrosti ... 71 5.36 Histogram popravkov komponent stanja sistema - komponente pospeška . . 71 5.37 Pravokotna oddaljenost filtriranih (rdeče) in merjenih (modra) vrednosti od referenčne tirnice v horizontalni ravnini in interval 1 cm......... 72 5.38 Pravokotna oddaljenost filtriranih (rdeče) in merjenih (modro) vrednosti z 1 cm 5.39 Hitrost v vseh treh smereh (modro-filtrirane vrednosti, rdeče-referenčne vrednosti) in interval 2a (modro)....................... 74 5.40 Pospešek v vseh treh smereh (modro-filtrirane vrednosti, rdeče-referenčne vrednosti) in interval 2a (modro)....................... 74 5.41 Napaka hitrosti v vseh treh smereh ter interval a (zelena polna črta) in interval 2a (zelena prekinjena črta)...................... 75 5.42 Napaka pospeška v vseh treh smereh ter interval a (zelena polna črta) in interval 2a (zelena prekinjena črta)............................................75 5.43 Standardne deviacije hitrosti v vseh treh smereh..............................76 5.44 Standardne deviacije pospeška v vseh treh smereh............................76 5.45 Direktna opazovanja............................................................78 5.46 Indirektna opazovanja ali meritve KF..........................................78 B.l Razlike v merjenih dolžinah, d(i+l)-d(i) ...................104 B.2 Razlike v opazovanih horizontalnih smereh, hz(i+l)-hz(i)..........105 B.3 Razlike v opazovanih zenitnih razdaljah, zr(i+l)-zr(i)............106 SIMBOLI "fc -k a △am A B(t ) C d dk Dfc em f (•) Fk G Rnxn Gk hz Ho Ha Hk G Rmxn I KF Kk G Rnxm LKF lk m n nr Nk O Pfc" P+ P k Qk g nnxn Rn Rl Rk g nmxm △t u v k G Rmx1 zr pospešek največja sprememba pospeška znotraj časovnega intervala časovno zvezna matrika prehoda stanja časovno zvezna matrika vhoda matrika nadzora poševna dolžina izboljšava meritev kovariančna matrika izboljšave meritev a priori napaka a posteriori napaka pričakovana vrednost nelinearna funkcija sistema matrika prehoda stanja diskretna matrika vhoda horizontalna smer ničelna domneva alternativna domneva nelinearna funkcija meritev merska matrika enotska matrika Kalmanov filter/Kalmanovo filtriranje Kalmanova matrika izboljšave linearni Kalmanov filter vektor direktnih meritev tk število neznank / število komponent stanja sistema število matematičnih operacij matrika parcialnih odvodov matrika opazovanosti a priori kovariančna matrika napake a posteriori kovariančna matrika napake kovariančna matrika popravkov stanja sistema kovariančna matrika šuma sistema ostanek Taylorjeve vrste kovariančna matrika direktnih meritev kovariančna matrika meritev časovni interval vektor nadzora vnosa merski šum vektor popravkov stanja sistema zenitna razdalja zk G nmx1 wk G nnx1 Wk xk G Rnx1 x- x+ xk rk Q2 Qd,k QV 2*'k Xr,1-a k(P) a ad ahz azr am vektor meritev za tk šum sistema sprememba pospeška vektor neznank/stanja sistema za tk a priori ocena stanja sistema a posteriori ocena stanja sistema vektor izboljšave testna statistika v domeni meritev testna statistika v domeni stanja sistema X2-porazdelitev z r prostostnimi stopnjami in stopnjo značilnosti a pogojno število standardna deviacija standardna deviacija dolžine standardna deviacija horizontalne smeri standardna deviacija zenitne razdalje skalar intenzitete šuma sistema 1 UVOD Zajem deformacij grajenih ali naravnih objektov, tako v smislu merskega zajema kot analize rezultatov, danes zahteva interdisciplinarno delo strokovnjakov iz različnih področij, kot so gradbeništvo, geodezija, geologija, hidrologija, seizmologija, ... ter disciplin, ki se ukvarjajo s teorijo sistemov (matematika, statistika) in obdelavo signala (elektrotehnika). Interdisciplinarno delo posledično pogojuje podrobne zahteve glede opredelitev postopkov in analiz deformacij. Geodetske meritve so le del postopkov spremljanja deformacij in premikov. Ena od glavnih geodetskih nalog je zajem in analiza deformacij in premikov naravnih in grajenih objektov, ki se lahko obravnavajo na različne načine glede na upoštevanje časa in delujočih sil. Danes najpogosteje uporabljeni modeli so dinamični in kinematični modeli. Sodobni geodetski instrumenti omogočajo zajem kinematičnih procesov z visoko frekvenco meritev. Ta lastnost geodetskih instrumentov daje geodetskemu zajemu deformacij in premikov ključno prednost, t.j. opredelitev deformacij in premikov v prostoru in času. Z ustreznim deformacijskim modelom je vektor neznank (običajno so to prostorske koordinate referenčnih točk objekta), s pripadajočo oceno natančnosti, mogoče določiti v realnem času z upoštevanjem vseh časovnih zakasnitev. Pri obravnavi objekta kot dinamični ali kinematični proces, se predpostavlja, daje objekt stalno v gibanju. Ob predpostavki, da so spremembe objekta zajete z merskim sistemom, ki omogoča enolično določitev položaja (npr. TPS, GPS, laserski interferometer), za posamezni časovni trenutek nimamo na voljo nadštevilnih meritev. Klasični geodetski postopki naknadne obdelave, ki temeljijo na izravnavi nadštevilnih geodetskih opazovanj, v tem primeru niso uporabni. Potrebno je vključiti druge tehnike, ki omogočajo oceno opazovanj v realnem času. Kalmanov filter (KF) predstavlja tehniko naprednih geodetskih postopkov kinematičnih procesov. Ocena temelji na napovedovanju prihodnjega stanja sistema na osnovi znanega preteklega vedenja. 1.1 Razvrstitev deformacijskih modelov Geodetska deformaeijska analiza danes pomeni geodetsko analizo dinamičnih in kinematičnih procesov [Welsch, Heunecke (2001)], kar pomeni vključevanje dejavnikov vpliva in časovne odvisnosti. V tabeli 1.1 je podana razvrstitev in lastnosti deformacij skih modelov. V klasični geodeziji je najpogosteje uporabljen kongruenčni model. Deformacije in premiki objekta so definirani s spremembami položajev referenčnih točk objekta med dvema izmerama. V modelu niso vključene vzročne sile in časovna odvisnost. Ker se predpostavlja, daje objekt v času meritev v stanju mirovanja, je mogoče izvesti nadštevilne meritve. Nadštevilne meritve zagotavljajo klasično izravnavo geodetskih opazovanj [Bakker in sod. (1995)], katere rezultat je ocena neznank in pripadajoča natančnost. Izravnava je izvedena na osnovi metode najmanjših kvadratov in omogoča objektivno kontrolo geodetskih Deformacijski Kongruenčni Kinematični Statični Dinamični model model model model model Časovna brez modeliranja premiki so brez modeliranja premiki so komponenta funkcija časa funkcija časa Vzročne brez modeliranja brez modeliranja premiki so premiki so sile funkcija sil funkcija sil Stanje objekta zadostno stanje mirovanja v gibanju zadostno stanje mirovanja v gibanju Preglednica 1.1: Razvrstitev in lastnosti deformacijskih modelov [Welsch, Heunecke (2001)] opazovanj, ki je zahteva pri vseh geodetskih delih ne glede na zahtevano natančnost, [Le-ick (1990)]. Novejša geodetska deformacijska analiza vse bolj vključuje geodetsko analizo dinamičnih in kinematičih procesov, z vključitvijo sil in časovnega poteka nastalih deformacij in/ali premikov. Modeli, ki opisujejo deformacijo ali premik kot funkcijo sil in časa, se imenujejo dinamični modeli. Prva poenostavitev dinamičnih modelov so kinematični modeli, ki opisujejo deformacijo zgolj kot funkcijo časa. V preteklosti se je kinematika deformacijskih procesov obravnavala s transformacijo kinematičnega problema v zaporedje statičnih meritev. Moderni geodetski sistemi, kot sta globalni sistem določanja položaja, angl. Global Positioning System (GPS), in terestrični sistem določanja položaja, angl. Terrestrial Positioning System (TPS), omogočajo samodejno, ponavljajoče spremljanje objektov v gibanju in kot taki omogočajo zajem premikov in deformacij skoraj v realnem času. Geodetski instrumenti ponujajo možnost izpolnitve smernic v inženirski geodeziji, ki poleg geometričnih sprememb predvidevajo tudi in predvsem opis sprememb objekta v prostoru in času kot posledici sil. Tak postopek zahteva uporabo dinamičnih modelov, ki obravnavajo objekt, vzročne sile in nastale deformacije kot celoto. S poenostavitvijo dinamičnih modelov so izvedeni kinematični modeli. V kinematičnem pristopu obravnave koordinat točk kot funkcije časa je deformiranje opisano s koordinatami, hitrostjo in pospeškom. Prav tako so področja kot so računalništvo, nove merske tehnike in novi algoritmi vrednotenja bistveno pripomogli h kvalitetnejši deformacijski analizi. 1.2 Kalmanov filter V primeru, da opazujemo objekt, ki je v nenehnem gibanju, in pod predpostavko, da je kinematičen proces opazovan samo z elektronskim tahimetrom kot terestričnim sistemom za določanje položaja, ni mogoče izvesti nadštevilnih opazovanj, ki bi omogočala izračun natančnosti določitve neznank. Klasična geodetska izravnava z večjim številom opazovanj, kot jih je nujno potrebnih za enolično določitev neznank sistema, ne more biti izvedena. Za ovrednotenje neznanega vektorja stanja sistema in pripadajoče statistične ocene v realnem času je potrebno uporabiti druge metode. V kinematičnih procesih imamo opravka s časovnimi vrstami in standardnimi tehnikami odstranitve šuma iz časovne vrste, kot so filtriranje in glajenje. V smislu metode najmanjših kvadratov je odstranitev šuma optimalna z uporabo IvF. IvF reši problem enolične določitve neznank. Filter računsko zagotovi nadštevilna opazovanja, ki so potrebna za odstranitev vplivov merskih napak [So-renson (1970)] in izračun natančnosti ocene neznank. Teorija IvF je bila izpeljana iz del Gaussa. Ivolmogorova in izsledkov Wienerjevega filtra v začetku 60. let prejšnjega stoletja. Wienerjev pristop filtriranju po metodi najmanjših kvadratov je uporabljen tudi v IvF. Oba modela sta osnovana na definiciji utežne funkcije: kako utežiti vhodne podatke, da bodo zagotavljali pridobitev najboljše ocene neznank v danem trenutku. IvF je optimalni rekurzivni linearni algoritem za obdelavo meritev, v katerih je prisoten šum. Za izračun ocene vektorja neznank/stanja sistema algoritem vključi vsako opazovanje. ne glede na natančnost. Algoritem vključi v model sistema in model meritev vse razpoložljive informacije o. slika 1.1: - dinamiki sistema (časovnem spreminjanju in silah). - merskih sistemih. - statističnih karakteristikah merskih napak in napak sistema ter - začetnih vhodnih podatkih. Izhodna količina Kalmanovega filtra je ocena stanja sistema v danem trenutku, slika 1.1. Napaka sistema Merska napaka Ocena vektorja stanja sistema Slika 1.1: Diagram sistem-meritve-Ivalmanov filter, povzeto po [Gelb (1974)]. Figure 1.1: Block diagram depicting system-measurement-Ivalman filter, derived from [Gelb (1974)] Začetki praktične uporabe KF segajo na področje navigacije, to je določanja tirnice gibanja objekta v prostoru in času. Algoritem KF se danes uporablja na področjih geodezije, elektrotehnike, medicine, ..., povsod v nalogah, kjer so neznane količine določene kot funkcija pogrešenih opazovanj. Vhodne količine so določene z deterministično komponento in stohastičnimi motnjami. Namen metode KF je filtriranje podatkov, odkrivanje sistematičnih in grobih pogreškov ter določanje vpliva slučajnih pogreškov. Algoritem upošteva vse razpoložljive podatke o meritvah, predhodno znanje o sistemu in merskem instru-mentariju, za pridobitev ocene stanja sistema na način, ki statistično minimizira napake. Učinkovitost filtra se je izkazala uspešna zlasti pri nalogah, kjer ni na voljo nadštevilnih meritev, t.j. v geodetskih nalogah navigacije ter stalnem zajemanju deformacij in premikov. Pri obdelavi podatkov rekurzivni algoritem KF ne zahteva shranjevanje vseh meritev in stanj sistema. Filter namreč izvede izračun neznanih količin v danem trenutku samo na osnovi podatkov predhodnega stanja sistema. Ta lastnost je ena temeljnih prednosti algoritma, kadar se uporablja za velike sisteme z veliko količino podatkov. 1.3 Opis problema Glede na hiter razvoj merskih sistemov, ne le geodetskih, in potrebe po interdisciplinarni obravnavi deformacij in premikov morajo biti tudi sodobni geodetski merski sistemi in postopki poznani v vseh podrobnostih. Zmogljivosti kinematičnih postopkov GPS so dandanes že dobro poznane. Po drugi strani pa tudi terestrična geodezija s sodobnimi instrumenti ponuja možnosti za obravnavo deformacij in premikov naravnih in grajenih objektov kot kinematičnih procesov. Vključitev časovne komponente neposredno v obdelavo terestričnih kinematičnih opazovanj zahteva dobro poznavanje zmogljivosti sodobnih elektronskih tahimetrov in postopkov obdelave terestričnih kinematičnih opazovanj. Na voljo je veliko literature, ki opisuje teoretične osnove ocenjevanja dinamičnih in kinematičnih sistemov in podajajo numerične primere filtriranja. V tem delu sta bili kot teoretična osnova uporabljeni predvsem deli: [Bar-Shalom in sod. (2001)] in [Simon (2006)]. Poleg navedenih temeljnih del so bila za teoretično podlago uporabljena naslednja dela: [Gelb (1974)], [Brown, Hwang (1992)], [Blackman, Popoli (1999)], [Grewal, Andrews (2001)], [Chui, Chen (1999)], [Welch, Bishop (2004)] in [Eibeiro (2004)]. Omenjena dela vključujejo številne kratke, vendar za razumevanje algoritma izjemno uporabne, numerične primere uporabe KF. Razvoj geodetskih postopkov in metod spremljanja in vrednotenja deformacij je opisan v dveh člankih: [Welsch (2002)] in [Welsch, Heunecke (2001)] ter knjigi [Welsch in sod. (2000)]. V doktorski disertaciji [Giilal (1997)] je podrobno opisana uporaba KF na primeru geodetskih opazovanj grajenega objekta - jezu. Nadaljnje praktične izvedbe KF so podane v [Busse, How (2002)], [Eichhorn (2005)], [Gao in sod. (2005)], [Negenborn (2003)], [Lippitsch, Lasseur (2006)] in [Pelzer (1987)]. Matematične osnove, potrebne za razumevanje KF, so podane v [Golub, Van Loan (1991)] in [Bronštejn, Semendjajev (1967)]. V knjigah [Ristic in sod. (2004)] in [Haykin (1999)] so predstavljene adaptivne tehnike filtriranja. Zmogljivosti sodobnih elektronskih tahimetrov so opisane v [Stempfhuber (2004)], [Bayoud (2006)] in [Scholderle (2006)]. Glavna naloga dela je strnjena v dveh delih: 1. Izvedba kinematičnega procesa - simulacija kinematičnega procesa in - zajem le-tega z elektronskim tahimetrom; 2. Zapis modela vrednotenja - razviti ustrezen model vrednotenja in - preveriti ustreznost modela z elementi notranjega zaupanja, statističnimi testi in neodvisnim referenčnim okvirjem. Glavni poudarek dela je na opisu terestričnega kinematičnega procesa. Predmet doktorske disertacije je razviti primeren model vrednotenja opazovanj, zajetih s sistemi TPS, še posebej v nalogah, kjer je potrebno izpolniti visoke zahteve po natančnosti. Pri delu s praktičnimi problemi, ki vključujejo diskretne podatke, zlasti pri vrednotenju le-teh v realnem času, je pomembno, da je uporabljena metoda računsko izvedljiva, kot tudi matematično pravilna. Za preverjanje uporabljene metode so običajno izvedene predhodne simulacije in testi meritev. Druga možnost je vzpostavitev neodvisnega sistema kontrole. Simulacija kinematičnega procesa V okviru doktorske disertacije je bil simuliran kinematični proces na poznani referenčni tirnici v geodetskem laboratoriju Tehniške univerze Miinchen. Referenčna tirnica, na kateri je nameščen voziček, je definirana v lokalnem koordinatnem sistemu laboratorija. Na voziček je bil nameščen reflektor. Kinematične meritve so bile izvedene z elektronskim tahimetrom proizvajalca Leica Geosystems, ki omogoča samodejno prepoznavanje tarče in sledenje. Uporabljen je bil pripadajoči reflektor, proizvajalca Leica Geosystems. Meritve so bile zajete s programom Visual Basic, pripravljenim posebej za namen doktorske disertacije. V program je bila vključena knjižnica Leica Geosystems GeoCOM. Model vrednotenja in preizkušanje le-tega Glavni namen doktorske disertacije je bil razviti ustrezen model vrednotenja in ocena njegove učinkovitosti za ocenjevanje kinematičnega procesa, ki smo ga zajeli z elektronskim tahimetrom. Glede na dejstvo, da je bil kinematični proces zajet samo z enim merskim sistemom - elektronskim tahimetrom, ni bilo mogoče izvesti nadštevilnih meritev. Za pridobitev kinematičnega položaja premikajočega se reflektorja in predvsem parametrov, ki podajajo statistično oceno pridobljenega vektorja neznank, je potrebno uporabiti tehniko filtriranja. Na podlagi podrobne študije teoretičnih osnov filtriranja, še posebej v literaturi [Bar-Shalom in sod. (2001)] in [Simon (2006)], so bile preizkušene različne možnosti izvajanja KF. Čeprav sta algoritem in delovanje KF relativno enostavna za razumevanje, je vse prednosti in pomanjkljivosti mogoče razumeti le z dobrim znanjem matematičnih osnov in na podlagi praktičnih izkušenj. Na osnovi podrobnega poznavanja gibanja vozička vzdolž tirnice in zmogljivosti uporabljenega instrumenta za kinematične meritve bo razvit ustrezen model vrednotenja. V okviru zadnje in najpomembnejše naloge dela so bile izvedene različne kontrole modela, kot so ocene izhodnih parametrov modela, statistični testi in neodvisna kontrola. 1.4 Vsebina doktorske disertacije Doktorska disertacija vključuje 6 poglavij. Poglavje 1 opisuje razvrstitev deformacijskih modelov in kratek opis KF. Poglavje 2 podaja matematični model KF s poudarkom na linearnem KF in Kalmanovi matriki izboljšave. Podpoglavje 2.2 vključuje postopek diskretizacije - prehoda iz zveznega v diskretni časovni linearni model z definicijo kovariančne matrike procesnega šuma. V nadaljevanju, v podpoglavju 2.3, je na kratko opisan postopek v primeru nelinearnih sistemov. Podan je postopek linearizacije z razvojem v Taylorjevo vrsto. Poglavje 3 je namenjeno statistični analizi. Bralca seznani z osnovami inženirskega vrednotenja sistemov. V poglavju so predstavljeni statistični postopki, ki so lahko uporabljeni v geodetski analizi ocenjevanja stanja sistema v diskretnih linearnih časovno odvisnih modelih, kjer je vektor stanja sistema funkcija časa. Poglavje opisuje postopke vrednotenja KF, ki so uporabljeni za vrednotenje rezultatov v računskem delu doktorske disertacije. Poglavji 4 in 5 predstavljata osrednji del doktorske disertacije. Poglavje 4 - Kalmanov model meritev in sistema (računski primer) - opisuje zastavljen problem. V podpoglavju 4.2 je opisana praktična zasnova naloge: zagotovitev referenčnega okvirja, uporabljen in-strumentarij in podroben opis funkcionalnega in stohastičnega modela, t.j. diskretnega modela Wienerjevega procesa s pospeškom. Poglavje 5, z naslovom Vrednotenje diskretnega modela Wienerjevega procesa s pospeškom (DMWPP), razširi tehnike vrednotenja, teoretično predstavljene v poglavju 3, na primer časovno odvisnih količin v računskem primeru naloge. Obravnavan je problem ocenjevanja vektorja neznank v stohastičnem linearnem časovno odvisnem modelu. Vrednotenje diskretnega modela Wienerjevega procesa s pospeškom - DMWPP je izvedeno s pokazatelji notranje zanesljivosti - podpoglavje 5.2, s statističnimi testi - podpoglavje 5.3 in neodvisno referenčno tirnico - podpoglavje 5.4. Kot pokazatelja notranje zanesljivosti sta izračunani matrika nadzora in matrika opazovanosti. Opisane so lastnosti a posteriori kovariančne matrike vektorja neznank in matrike Kalmanove izboljšave. Problem skladnosti, ki je bistven pri vrednotenju optimalnosti modela vrednotenja v vsaki izvedbi, je obravnavan v podpoglavju 5.3. Izračun skladnosti je izveden na osnovi popravkov meritev - podpoglavje 5.3.2, in na osnovi popravkov vektorja neznank - podpoglavje 5.3.3. Kot neodvisen parameter za oceno modela je uporabljena referenčna tirnica. Rezultati so podani v poglavju 5.4. V poglavju 6 je podan zaključek dela kot povzetek dognanj in prispevkov - podpoglavje 6.1. Podani so predlogi za nadaljnje delo in raziskovanje - podpoglavje 6.2 ter na koncu -podpoglavje 6.3 - zaključne opombe. Doktorska disertacije vključuje dve prilogi. Prva - Priloga A - vključuje rezultate programa CAPLAN za določitev stojišča instrumenta. Druga priloga - Priloga B - grafično prikazuje razlike v meritvah. 2 KALMANOV FILTER Na splošno lahko filtriranje opišemo kot postopek za ločevanje ene komponente ali določene vrednosti od druge. Problem filtriranja v inženirstvu izhaja s področja elektrotehnike, t.j. določitve signala znotraj nekega območja frekvenc in izločitev frekvenc, ki so izven opredeljenega območja. K problemu je mogoče pristopiti na dva načina: - z izbiro ustreznega merskega instrumentarija ali - z modeliranjem primernega matematičnega algoritma za vrednotenje merjenih in iskanih količin. V 40. letih prejšnjega stoletja je Norbert Wiener pristopil k problemu z matematičnega vidika v smislu določitve filtrirane frekvence, ki bi ustrezala optimalni ločitvi signala in motenj - šuma. Rezultat dela je bil Wienerjev filter, ki ga označujejo naslednje značilnosti, [Brown, Hwang (1992)]: - predpostavka, da sta signal in šum slučajna procesa z znanimi spektralnimi značilnostmi oziroma poznano avtokorelacijsko in navzkrižno korelacijsko funkcijo, - kriterij pridobitve optimalne rešitve je najmanjši srednji kvadrat napake, - rešitev temelji na določitvi optimalne utežne funkcije filtra. Rezultat Wienerjevega filtra je utežna funkcija, ki določi uteži vhodnih količin na tak način, da je rezultat optimalna ocena izhodnih količin v danem trenutku. Vendar pa Wienerjev pristop ni podajal primerne rešitve za diskretne probleme, prav tako ne za zapletene časovno sprejemljive sisteme z večjim številom vhodnih in izhodnih količin. Leta 1960 je R. E. Kalman podal alternativno rešitev - filtriranje po kriteriju najmanjšega srednjega kvadrata napake. Glavni prednosti Kalmanovega zapisa in rešitve problema sta [Brown, Hwang (1992)]: - vektorsko modeliranje slučajnih procesov in - rekurzivna obdelava pogrešenih opazovanj kot vhodnih količin. Ideja KF je podana na sliki 2.1. Izračun Kalmanove matrike izboljšave Projekcija v naslednji korak VHODNE MERITVE Popravljena ocena neznank IZHODNA OCENA Popravljena kovariančna ocena Slika 2.1: Zanka KF, [Brown, Hwang (1992)] Figure 2.1: KF loop. [Brown. Hwang (1992)] KF je dobro poznano matematično orodje, ki daje odgovor na najpogosteje zastavljeno inženirsko vprašanje: kako pridobiti najboljšo oceno stanja sistema iz meritev, v katerih je prisoten šum. Odgovor na zastavljeno vprašanje ponuja algoritem za obdelavo podatkov, ki temelji na določitvi stanja sistema na osnovi pogrešenih meritev po metodi najmanjših kvadratov. Postopek zagotavlja pridobitev optimalne ocene stanja sistema in hkrati tudi mero. s kolikšno zanesljivostjo je ocenjeno stanje sistema tudi pravo stanje. KF je optimalen glede na katerokoli merilo v primeru linearnega procesa, kjer je v sistemu in meritvah prisoten nekoreliran. beli šum. z Gaussovo (normalno) porazdelitvijo in s srednjo vrednostjo enako nič [Xegenborn (2003)]: - nekoreliranost KF predpostavlja, da sta šum v sistemu in meritvah neodvisna in njuna jakost nima medsebojnega vpliva. Beli šum je šum. ki je nekoreliran v času. torej neodvisen od svojih preteklih stanj. Posledica je. da so napake nekorelirane v času. - srednja vrednost enaka nič Srednja vrednost šuma enaka nič pomeni, da je šum v sistemu in meritvah slučajen šum. kar pomeni, da lahko ima pozitiven ali negativen predznak, v povprečju pa njegova vrednost znaša nič. - Gaussova oz. normalna porazdelitev Predpostavka Gaussove porazdelitve pomeni, da je šum modeliran z Gaussovo krivuljo. Gaussov šum je v celoti opisan s srednjo vrednostjo in varianco ter tako podaja vse razpoložljive informacije o šumu in funkciji pogojne gostote verjetnosti. - beli šum KF predpostavlja, da je začetna vrednost stanja sistema prav tako modelirana kot slučajna spremenljivka, z Gaussovo porazdelitvijo s poznano srednjo vrednostjo in varianco. Sum v meritvah in sistemu ter začetna vrednost stanja sistema so prav tako medsebojno neodvisne količine. Najosnovnejša tehnika KF je linearni KF, kjer je zveza med meritvami in neznanimi količinami linearna. KF omogoča tudi obdelavo nelinearnih sistemov; bodisi se nelinearnost pojavi v enačbi sistema in/ali v merskih zvezah. V primeru nelinearnosti je potrebna linearizacija, kije izvedena z razvojem funkcije v Taylorjevo vrsto na dva načina [Brown, Hwang (1992)]: - linearizirani KF V primeru lineariziranega KF je linearizacija izvedena v poznanih približnih vrednostih, ki niso odvisne od meritev. - razširjeni KF V primeru razširjenega KF je linearizacija izvedena okoli ocenjenih stanj sistema, ki vključujejo informacije, pridobljene na osnovi meritev. Pri vseh tehnikah KF moramo poznati nekatere predhodne parametre modela, t.j. matrika prehoda stanja, merska matrika, ki podaja zvezo med vektorjem meritev in vektorjem neznank, kovariančna matrika procesnega šuma in kovariančna matrika meritev. 2.1 Linearni Kalmanov filter V tem poglavju je opisan model diskretnega linearnega KF (LKF). LKF obravnava problem ocenjevanja stanja sistema in pogrešenih opazovanj, ki so lahko v celoti opisana z linearnim modelom. LKF obravnava slučajni kinematični proces, ki ga lahko modeliramo z linearno enačbo sistema: Xk+i = Fk • xk + wfc (2.1) in linearno enačbo meritev. Zk = Hk • Xk + Vk, (2.2) kjer je: xk G Rnxl ... vektor neznank sis tema za tk zk G Rmxl ... vektor meritev za tk wk G Rnxl ... šum sistema vk G Rmxl ... merski šum Fk G Rnxn ... matrika prehoda stanja Hk G Rmxn ... merska matrika n m ... število meritev za tk Vrednosti n in m sta v celotnem procesu konstantni. Algoritem oceni vektor neznank Xk na osnovi diskretnih opazovanj zk, v časovnem trenutku tk. Matrika prehoda stanja Fk v enačbi sistema (2.1) povezuje dve zaporedni stanji sistema; Xk+1 v trenutku tk+1 in Xk v trenutku tk. Merska matrika Hk v enačbi meritev (2.2) povezuje trenutno stanje sistema xk z meritvami zk. Matriki Fk in Hk sta sicer lahko funkciji časa, običajno pa sta konstantni. Vektorja wk in vk predstavljata slučajni beli šum sistema in merski šum z normalno porazdelitvijo; wk ~ N(0, Qk) in vk ~ N(0, Rk). Pripadajoči kovariančni matriki Rk in Qk vektorjev wk m vk sta podani z: E W • wD = { Q ' = 1 (2.3) e ivk-n^ Rk ■ i=k (M E[wk • vT] = 0, za vse k in i. (2.5) V procesu sta matriki Rk in Qk konstantni ali pa sta funkciji časa. Kovariančna matrika meritev Rk G Rmxm je določena na osnovi poznane natančnosti merskega sistema ali pa je določena na osnovi predhodnih meritev. Določitev kovariančne matrike šuma sistema Qk G Rnxn je veliko težja, saj model, s katerim je opisan proces, ne moremo istočasno nadzorovati. Enačbe KF lahko zapišemo v dveh skupinah: enačbe napovedovanja ali enačbe časovne korekcije in enačbe meritev ali enačbe merske korekcije. V nadaljevanju tega poglavja bomo oceno stanja sistema označili z X. Kasneje bo znak zaradi poenostavljenega zapisovanja opuščen, saj imamo vedno opravka z dejanskimi numeričnimi vrednostmi, ki pa podajajo oceno stanja sistema. V različni literaturi lahko zasledimo oba zapisa. Proces KF zahteva poznavanje ocene začetnega stanja sistema oz. napovedano začetno vrednost, Xin pripadajočo a priori kovariančno matriko, P-, za časovni trenutek tk- Merske enačbe vključujejo nove meritve zk v a priori oceno X- za pridobitev izboljšane vrednosti oziroma za pridobitev a posteriori ocene stanja sistema, X+, in pripadajoče a posteriori kovarianče matrike, P+. Enačbe meritev sestavljajo: Kfc = P- • HT • (Hk • P- • HT + Rk)-1 (2.6) X+ = X- + Kk • (zk - Hk • X-) (2.7) P+ = (I - Kk • Hk) • P-. (2.8) A posteriori ocena v enačbi (2.7) je linearna kombinacija a priori ocene stanja sistema X-in utežene razlike med dejanskimi zk in napovedanimi Hk•X- meritvami. Uteži so določene s Kalmanovo matriko izboljšave, Kk G Rnxmy optimalno po izbranem kriteriju. Izbran je kriterij najmanjšega srednjega kvadrata napake. KF iterativno popravlja Kalmanovo Kk rešitvi. Enačbe napovedovanja projicirajo stanje sistema in njegovo oceno v naslednji korak ter omogočajo pridobitev a priori ocene za naslednji časovni trenutek tk+1} xfc+1 in P-+1. Enačbe napovedovanja sestavljajo: x-+i = Fk • x+ (2.9) P-+i = Fk • P+ • FT + Qk• (2.10) Celotna zanka KF je podana na sliki 2.2. Napovedovanje Popravljanje 4. Napovedovanje stanja sistema xk+1 5. Napovedovanje kovariančne matrike 1. Izračun Kalmanove matrike izboljšave K* 2. Popravljanje z meritvami zt 3. Izračun a-posteriori kovariančne matrike P J" Slika 2.2: Zanka KF Figure 2.2: KF loop Alternativne oblike določitve a posteriori kovariančne matrike P+ P k+ Josephova oblika a posteriori kovariančne matrike, podane z enačbo. [Bar-Shalom in sod. (2001)]: P+ = (I - Kk • Hk) • P- • (I - Kk • Hk)T + Kk • Rk • KT. (2.11) Izraz je sicer računsko zahtevnejši, vendar manj občutljiv na zaokrožitvene napake. Izraz ne daje negativnih lastnih vrednosti, kot se to lahko zgodi v izrazu (2.8) kot posledica odštevanja. Izraz (2.11) zagotavlja simetričnost in pozitivno definitnost matrike P+, dokler je P^T simetrična pozitivno definitna matrika. Izraz (2.8) je računsko enostavnejši kot izraz (2.11). Njegova oblika ne zagotavlja simetričnosti in pozitivne definitnosti matrike P+ Pk 2.1.1 Kalmanova matrika izboljšave Kk V tem podpoglavju je podana izpeljava Kalmanove matrike izboljšave, Kk, in njen pomen. Glavni cilj KF je pridobitev izraza za izračun a posteriori ocene vektorja stanja sistema X+ kot linearne kombinacije a priori ocene stanja sistema Xk in utežene razlike med pravimi vrednostmi meritev zk in napovedanimi meritvami Hk• Xk, podane v enačbi (2.7). Razlika dk = Zk - Hk • Xk (2.12) je izboljšava meritev ali popravek meritev. Vrednost izboljšave odraža neskladnost ali razliko med napovedanimi opazovanji Hk • Xk in pravimi opazovanji zk. Popravek ali izboljšava, ki ima vrednost enako 0 dk = 0, pomeni popolno skladnost med vektorjem napovedanih in dejanskih meritev [Welch, Bishop (2004)]. Kk bo zagotavljala najboljšo a posteriori oceno vektorja stanja sistema X+. Kot kriterij je izbran najmanjši srednji kvadrat napake. Kar pomeni, da želimo v po- Kk elemente vzdolž diagonale a posteriori kovariančne matrike vektorja neznank P+, ki predstavljajo variance (standardne deviacije) elementov vektorja stanja Xk A priori e'k in a posteriori e+ napako vektorja stanja sistema Xk zapišemo kot: ekk = Xk - Xk (2-13) e+ = Xk - X+. (2.14) Pripadajoči a priori in a posteriori kovariančni matriki sta: Pk = E[ek • (ek)T] = E[(x* - Xk) • (Xk - Xk)T] (2.15) in P+ = E[e+ • (e+)T] = E[(x* - X+) • (x* - X +)T], (2.16) kjer je E[•] pričakovana vrednost. Z vključitvijo enačbe meritev (2.2) v izraz za a posteriori oceno stanja sistema (2.7) in z nadomestitvijo pridobljenega izraza za X + v enačbi (2.16) dobimo razširjen izraz za P+: P + _ E[[(xk - x—) - Kk • (Hk • xk + Vk - Hk • x-)] • • [(xk - X-) - Kk • (Hk • xk + Vk - Hk • X-)] ]. (2.17) Z zapisom izraza za pričakovano vrednost in ob upoštevanju, daje a priori napaka (xk-xk ) Vk P + (I - Kk • Hk) • P- • (I - Kk • Hk) + Kk • Rk • K -T k. (2.18) Enačba (2.18) je splošen izraz za a posteriori kovariančno matriko in ustreza katerikoli matriki Kk. Cilj KF je minimizacija sled i matrike P+, saj le-ta predstavlja vsoto srednjih kvadratov popravkov v oceni vseh elementov vektorja neznank oz. vektorja stanja sistema Kk sled (tr) matrike P+ po matriki Kk: d(tr(P+)) dKk 0. (2.19) Za odvajanje gornjega izraza je potrebno poznati dve operaciji matričnega odvajanja, [Brown, Hwang (1992)]: d(tr(AB)) dA B T (2.20) d(tr(ACAT)) dA 2 AC, (2.21) AB C C = CT. Odvod skalarja po matriki je definiran z izrazom: ds dA Izraz (2.18) sedaj zapišemo kot: ds dan ds dai2 ds ds da2i da22 (2.22) P + _ p- P- - Kk • Hk • P- + Kk • (Hk • P- • HT + Rk) • K -T k. (2.23) Z uporabo enačb za matrično odvajanje pridobimo izraz za odvod sledi matrike P+ po Kk d(tr(P+)) dKk -(Hk • P-)T + 2 • Kk • (Hk • P- • HT + Rk) Kk d(tr (P+)) dKk (2.24) (2.25) in jo zapišemo kot: k k k 0 Kk = Pk • HT • (Hk • Pk • HT + Rk)k1. (2.26) Iz izraza Kalmanove matrike izboljšave v enačbi (2.26) lahko povzamemo [Welch, Bishop (2004)]: Rk Kk vergira proti H^J"1: lim Kk = H1. Rk k zk obravnava kot bolj natančna v primerjavi z napovedano vrednostjo meritve Hk • Xk- - Po drugi strani, v primeru, ko se a priori kovariančna matrika vektorja stanja sistema Pk- Kk lim Kk = 0. zk obravnava kot manj natančna v primerjavi z napovedano vrednostjo meritve Hk • Xk . Kalmanovo matriko izboljšave lahko imenujemo tudi utežna matrika. 2.2 Diskretizacija: Prehod iz zveznega problema na diskretni linearni problem Večina sistemov v realnem svetu je opisana z linearnim ali nelinearnim zveznim modelom. Ker imamo v praksi običajno na voljo direktna ali indirektna opazovanja stanja sistema izvedena v diskretnih časovnih trenutkih, je potrebno izvesti ustrezno transformacijo iz zveznega sistema v diskreten sistem. Sistem je opisan z zveznimi diferencialnimi enačbami na eni strani, medtem ko se vhodne količine spremenijo le v diskretnih časovnih trenutkih (npr., meritve pridobimo v časovnih trenutkih t0, t1y t2, ■ ■ ■, tN)- Cilj je pridobitev ocene stanja sistema in pripadajoče ocene natančnosti v teh časovnih trenutkih. Zvezni stohastični linearni sistem opišemo z diferencialno enačbo, [Simon (2006)]: X = A • X + B • u + w, (2.27) kjer je x vektor stanja sistema, A časovno zvezna, matrika prehoda stanja, u vektor Bw po Gaussovi porazdelitvi, s srednjo vrednostjo enako nič in kovariančno matriko Qk. V primeru, da sta matriki A in B konstantni, je rešitev enačbe (2.27) za X(t) v izbranem tk x(tk) = eA(tk-tk-1) • x(tk-i) + /tk eA(tk-T) • [B(t) • u(r) + w(t)] • dr. (2.28) Jtk-1 Pri modeliranju diskretnega sistema predpostavimo, da so vhodne količine odsekoma konstantne, t.j. u(t) = uk za t G [tk-i,tk]. Zapišimo: At = tk - tk-i (2.29) xk = x(tk) (2.30) uk = u(tk). (2.31) Ob upoštevanju gornjih definicij zapišemo enačbo (2.28) kot xk = eA'At • xk-i + / eA ničelno domnevo ne zavrnemo. Povzamemo lahko, da ni statistično značilnega odstopanja med modelom sistema in meritev ter ni potrebno poseganje v filter. V primeru, da je odkrita večja vrednost izboljšave, ničelno domnevo zavržemo. V naslednjem koraku je potrebno poiskati vzroke neskladja. Razlog za neskladnost je lahko v opazovanjih, v napovedanih vrednostih ali stohastičnem modelu filtra [Eichhorn (2005)]. 3.2.2 Testna statistika v domeni vektorja stanja sistema V domeni vektorja stanja sistema preizkušamo, ali je filtrirana vrednost stanja sistema primerljiva s predhodnim znanjem o sistemu. Da bi lahko odgovorili na zastavljeno vprašanje, moramo analizirati razliko med filtrirano x++ in napovedano vrednostjo X- stanja sistema. Z upoštevanjem enačb (2.7) in (2.12) je popravek stanja sistema vx,k zapisan kot: vx,k = X++ - X- = Kk • dk. (3.7) Pripadajoča kovariančna matrika je dana z [Lippitsch, Lasseur (2006)]: PVx,k = Kk • Dk • KT. (3.8) Popravek stanja sistema ima normalno porazdelitev vx,k ~ N(0, PVxk). Popravek stanja sistema predstavlja razliko med oceno stanja sistema pred meritvijo in po njej. V primeru, da ima popravek veliko vrednost, pomeni, da je napovedana vrednost slaba ali da je prisotno slabo opazovanje. Ničelna domneva je definirana kot H0: E(vx,k) = 0 in preizkušana z alternativno domnevo Ha: E(vx,k) = 0. Pripadajočo testno statistiko zapišemo z izrazom [Lippitsch, Lasseur (2006)]: vT P-1 v k = ^.^o^ , (3 9) kjer je a0 referenčna varianca. Testna statistika ^Vx k je porazdeljena po porazdelitvi x2 s številom prostostnih stopenj n (število neznank v časovnem koraku tk) in a tveganjem. Razmerje verjetnosti je podano z izrazom: p{^Vxk < xh-a\Ho) = 1 - a. (3.10) 3.3 Zmožnost nadzora in opazovanosti sistema Meri nadzor (angl. controllability) in opazovanost (angl. observability) sistema sta ključnega pomena v teoriji nadzora sistemov. Nadzor sistema nam pove, kako dobro se lahko poljubno stanje sistema približa želeni vrednosti. Opazovanost sistema definira začetne pogoje, t.j. določitev začetnih pogojev na osnovi meritev izhodnih količin [Simon (2006)]. 3.3.1 Zmožnost nadzora sistema Diskretni in deterministični sistem je popolnoma nadzorovan, če za poljubno stanje sistema obstaja vhodno zaporedje, ki s končnim številom korakov zagotavlja prehod sistema iz začetnega stanja do poljubno izbranega stanja sistema. Za diskretni linearni časovno invariantni sistem z dolžino vektorja stanja sistema n je pogoj nadzorovanosti izpolnjen, če ima matrika nadzora C = [ G F • G F2 • G ... F(n-1) • G ] (3.11) n s konstantnim časovnim intervalom, so matrike F, G in H neodvisne od časa [Grewal Andrews (2001)]. Pogoj nadzorovanosti določa, da je šum sistema prisoten v vsakem stanju sistema in preprečuje, da bi kovarianca stanja sistema konvergirala nič. Ta pogoj zagotavlja, da je kovariančna matrika pozitivno definitna, t.j. da so lastne vrednosti pozitivne [Bar-Shalom in sod. (2001)]. 3.3.2 Zmožnost opazovanosti sistema Deterministični sistem je popolnoma opazovan, če je začetna vrednost stanja sistema lahko ponovno pridobljena povsem enolično na podlagi končnega števila opazovanj izhodnih stanj in poznanih vhodnih stanjih. Za diskretni linearni časovno invariantni sistem z n vanosti O H H • F H • F2 H F(n-1) (3.12) n 4 KALMANOV MODEL MERITEV IN SISTEMA (RAČUNSKI PRIMER) 4.1 Opis problema Naloga doktorske disertacije je realizirati uporabo KF za obdelavo terestričnih geodetskih meritev. Glavni namen je razviti KF model za poseben kinematični proces, ki ga opazujemo z elektronskim tahimetrom, in ovrednotenje numeričnih rezultatov, kijih pridobimo s filtrom. Vrednotenje modela je izvedeno s statističnimi testi, parametri notranjega zaupanja in na osnovi poznanega neodvisnega referenčnega sistema. KF omogoča nov pristop pri obdelavi geodetskih meritev in postopkih deformacijske analize. Omogoča obdelavo geodetskih referenčnih točk kot kinematičnih, ne samo kot statičnih, kot je to običaj v klasični geodeziji. Na drugi strani je potrebno dati velik poudarek na pomembnost razvoja geodetskega terestričnega instrumentarija, tako v optičnem, mehanskem kot elektronskem smislu. V delu je poseben poudarek na modernih elektronskih tahimetrih (angl. Terrestrial Positioning Systems - TPS) proizvajalca Leica Geosystems. 5 teoretičnega vidika je glavni poudarek dela na statistični analizi časovno odvisnih procesov, t.j. statistični analizi razvitega modela, in potrditev rezultatov s poznanim referenčnim sistemom. Slednje je podrobneje predstavljeno v poglavju 5 - Vrednotenje diskretnega modela Wienerjevega procesa s pospeškom. 4.2 Vzpostavitev praktičnega primera Glavna ideja dela je, kot že omenjeno v poglavju 1.3: - simulirati kinematični proces, - zajeti proces z elektronskim tahimetrom, - razviti model vrednotenja in - preveriti ustreznost modela z neodvisnim referenčnim okvirjem. 4.2.1 Referenčni okvir Ustrezni neodvisni referenčni okvir za vrednotenje razvitega modela, t.j. vrednotenje izhodnih količin glede na različne vhodne količine, smo zagotovili z ustrezno referenčno tirnico. Geodetski laboratorij Katedre za geodezijo na Tehniški univerzi Miinchen (TUM), slika 4.1, ima v lasti tako tirnico. Geodetski laboratorij na TUM se ukvarja s številnimi geodetskimi nalogami. Poleg širokega spektra del na področju kalibracije geodetskih instrumentov se ukvarja tudi z inovativnimi raziskovalnimi nalogami, katerih rezultati so kot velik doprinos na področju poznavanja zmogljivosti sodobnih geodetskih instrumentov in postopkov. Trenutno je poudarek na kalibraciji in preizkušanju kvalitete natančnih geodetskih instrumentov [Ge-odeticLaboratory]: Slika 4.1: Geodetski laboratorij Katedre za geodezijo. TUM. z referenčno tirnico Figure 4.1: Geodetic laboratory at Chair for Geodesy. TUM. with reference trajectory - preizkušanje kinematičnih senzorjev, še posebej kalibracija tahimetrov. ki omogočajo sledenje in zajem kinematičnih procesov. - kalibracija preciznih invar nivelmanskih lat (s kodirano ali klasično razdelbo) za nivelman najvišje natančnosti. - kalibracija zelo dolgih merskih trakov za potrebe višinomerstva na jezovih. - preizkušanje terestričnih laserskih skenerjev. Kinematični proces je bil definiran kot problem sledenja. Sledenje je ocenjevanje položaja objekta, ki je v gibanju, na osnovi daljinskih meritev [Bar-Shalom in sod. (2001)]. Meritve so izvedene z uporabo enega ali več senzorjev na izbranih položajih. Kinematični proces je bil simuliran z vozičkom, sliki 4.1 in 4.2. ki se giblje vzdolž vodoravne kalibracijske tirnice dolžine 15 m. Tirnica je bila vzpostavljena posebej z namenom preizkušanja elektronskih tahimetrov s samodejnim prepoznavanjem reflektorja (avtomatsko prepoznavanje tarče - APT [Valh in sod. (2008)], angl. Automatic Target Recognition -ATR, oz. APT ) in sledenjem. Dosežena je bila zahtevana natančnost namestitve tirnice 1 mm v horizontalnem položaju . Višinska komponenta je konstantna znotraj zahtevane natančnosti vertikalne namestitve tirnice, ki je prav tako 1 mm (osebna komunikacija z Dr.-Ing. P.Wasmeier-jem, znanstveni asistent in vodja Geodetskega laboratorija TUM). Geometrija referenčne tirnice je določena na osnovi sistema teodolitov ECDS S (Electronic Coordinate Determination System), proizvajalca Kern/Leica [Scholderle (2006)]. Sistem ECDS S je bil uporabljen tudi za vzpostavitev referenčnega lokalnega koordinatnega sistema v laboratoriju, za definiranje večjega števila točk, ki so lahko uporabljene za določitev koordinat poljubnega stojišča instrumenta. Tirnica je nameščena tako, daje vzporedna osi y lokalnega koordinatnega sistema. Za kinematične meritve je najprimernejši 360° reflektor. V nalogi je bil uporabljen reflektor GRZJf. Leica Geosystems, slika 4.2, s konstanto 23.1 mm, kije bil nameščen na voziček. 360° 360° Izvedeno je bilo večje število kinematičnih opazovanj, pri različnih hitrostih vozička. Med meritvami je instrument postavljen na stabilnem stebru, katerega koordinate so bile določene s programom CAPLAN, slika 4.3, po metodi notranjega ureza z izravnavo. Uporabljene so bile meritve horizontalnih smeri, zenitnih razdalj in poševnih dolžin proti petim danim točkam. Rezultati so podani v preglednici 4.1 in Dodatku A. Y_StPt [m] X_StPt [m] Z_StPt [m] -14.640 1.116 -0.388 Oy [m] ax [m] a z [m] ±0.002 ±0.003 ±0.001 Preglednica 4.1: Koordinate stojišča instrumenta s standardnimi deviacijami. Priloga A Ml 8 M22 M 30 M60 X M50 M51 M63 ■M62 ■ M64 M61 L27 * L26 -L25 _> Y StandPunktS ... Stojišče instrumenta Slika 4.3: Položaj stojišča instrumenta (StandPunkt_S) in danih točk v referenčnem (y, x) sistemu laboratorija Figure 4.3: Instrument stand point (StandPunkt_S) and reference points, given in coordinate system (y, x) of laboratory 4.2.2 Instrumentarij Meritve so bile izvedene s sodobnim elektronskim tahimetrom ali sistemom TPS TCR.A1201 serije TPS1200 Professional Series proizvajalca Leica Geosystems (Heerbrugg. Švica). Tahimeter je geodetski instrument, v geodeziji najpogosteje uporabljen za terensko izmero za različne naloge, ki je sestavljen iz dveh glavnih komponent: - elektronskega teodolita za merjenje horizontalnih smeri in zenitnih razdalj ter - elektronskega razdaljemera (EDM) za merjenje poševnih dolžin. Moderni tahimetri merijo horizontalne smeri na osnovi elektro-optičnega skeniranja natančne digitalne kodne razdelbe. nameščene na steklenem limbu znotraj instrumenta. Osnovi princip merjenja dolžin z elektronskimi razdaljemeri je določitev velikosti dolžine na osnovi izmerjenega časa. v katerem elektromagnetno valovanje prepotuje razdaljo med začetno in končno točko. Instrument je izvor elektromagnetnega valovanja. Valovanje preko oddajne optike instrumenta usmerimo proti reflektorju. Valovanje pade na prizmo reflektorja, kjer se elektromagnetni valovi odbijejo nazaj proti instrumentu. Valovanje še enkrat prepotuje merjeno dolžino v obratni smeri in pade na sprejemno optiko instrumenta [Kogoj (2002)]. Signal, ki se vrne v instrument, je registriran in interpretiran v instrumentu. Podatki o oddanem in odbitem signalu so pretvorjeni v dolžino. Skoraj vsi sodobni elektronski instrumenti omogočajo merjenje dolžine brez reflektorja (angl. reflek-torless (E)). Podsistem za merjenje dolžine je eden najobčutljivejših sestavnih delov in pogojuje najvišjo mersko frekvenco pri kinematičnih meritvah. Avtomatski elektronski tahimetri (A) so opremljeni z motoriziranimi horizontalnimi in vertikalnimi pogoni ter senzorji za samodejno prepoznavanje reflektorja. V okviru tehnologije APT obstajata dva sistema: avtomatsko viziranje tarče (AVT) in avtomatsko sledenje tarče (AST) [Kogoj (2002)]. Sistem AST omogoča hitro in natančno sledenje reflektorja v prostoru. Natančnost določanja položaja je v največji meri odvisna od natančnosti merjenja razdalje, ki je funkcija merskega programa EDM, in od natančnosti merjenja kotov. Natančnost izračunanih koordinat je odvisna tudi od uporabljenega reflektorja, oddaljenosti, kalibriranega instrumenta, upoštevanja sistematičnih pogreškov in od atmosferskih pogojev. Med izvajanjem meritev v laboratoriju so bili merjeni atmosferski pogoji, ki so bili zelo konstantni (T = 22.1° C, p = 961.6 mbar, e = 60%). Vpliv spremenjenih pogojev je pri kratkih dolžinah majhen in zato popravki niso bili vključeni v izračun. Tehnični podatki uporabljenega elektronskega tahimetra TCRA1201 Leica Geosystems so podani v preglednici 4.2. Uporabljen instrument ima tako razdaljemer z infrardečo svetlobo kot z vidnim laserskim žarkom. Instrument TCRA1201 Leica Geosystems Lastnosti dvoosni kompenzator samodejno prepoznavanje reflektorja; angl. Automatic Target Recognition - ATR merjenje brez reflektorja; angl. reflectorless - R Delovno območje [—20°C, +50°C ] Masa 4.9 kg Kotna hitrost 45°/ s Povečava objektiva 30 krat Premer objektiva 40 Najkrajša merjena razdalja 1.5 m Občutljivost libele dozna libela: 6'/2 mm elektronska libela: ločljivost 22' Kotne meritve Natančnost hz, zr1, n 1 '' nhz, zr — 1 podana s standardno deviacijo ISO 17123-3 Natančnost prikaza 0.1'' Merjenje dolžin Merska frekvenca 100 MHz/1.5 m Način merjenja fazni način EDM način dvo-osni, infra-rdeč laserski žarek Nosilno valovanje 780 nm Natančnost za sledilni način, podana s standardno deviacijo Od = 5 mm + 2 ppm ISO 17123-4 čas meritev < 0.15 s 1 ^-horizontalna smer, zr-zenitna razdalja Preglednica 4.2: Tehnični podatki instrumenta TCRA1201 Leica Geosystems [Tech. Data] Odprta arhitektura serije Leica Geosystems TP S1200 Series Elektronski tahimetri serije TPS1200 so eni najsodobnejših instrumentov za izvajanje geodetskih meritev. Z njimi lahko opravimo večino geodetskih nalog, še posebej z njihovimi integriranimi aplikacijami. Za izvedbo širšega spektra nalog in aplikacij je na voljo dodaten vmesnik za dostopanje do posameznih funkcij senzorjev TPS1200. Programska oprema sistema TPS1200 organizira in nadzoruje medsebojno delovanje več senzorjev in tvori okvir za aplikacije, ki jih je mogoče izvajati. Z novim vmesnikom Leica Geosystems omogoča odprto arhitekturo tahimetrov serije TPS1200. Vmesnik lahko uporabimo za reševanje posebnih nalog, v kolikor obstoječe rešitve Leica Geosystems ne zagotavljajo potrebne funkcionalnosti ali zgolj za izboljšanje obstoječih rešitev [GeoCOM Reference Manual V 1.20]. Uporabniške programe za elektronske tahimetre lahko krmilimo neposredno v instrumentu (angl. on-board) ali iz zunanjega računalnika (angl. off-board). Komunikacijski protokoli Leica Geosystems so prikazani na sliki 4.4. Programi so razviti s komunikacijskimi standardi za sodobne instrumente Leica Geosystems TPS1200 Professional Series. Aplikacije On-board so razvite na zunanjem računalniku v programskem vmesniku Leica Geosystems GeoBASIC, nato pa so prenesene v instrument. V primeru aplikacij off-board so le-te razvite in krmiljene preko zunanjega računalnika v programskem vmesniku Leica Geosystems GeoCOM. Časovne zakasnitve prenosa podatkov so odvisne od uporabljenega vmesnika. V doktorski disertaciji je bila dvostranska komunikacija z instrumentom izvedena z vmesnikom Leica GeoCOM. V programsko okolje Visual Basic je bila vključena knjižnica Leica Geosystems GeoCOM. Knjižnica vključuje podatkovne datoteke z vsemi potrebnimi konstantami, tipi podatkov in definicijami funkcij. Za namen vzpostavitve koordinatnega sistema oziroma izračun koordinat stojišča instrumenta so bile meritve najprej izvedene proti znanim točkam. Dolžine so bile izmerjene v načinu merjenje dolžin brez reflektorja. V naslednjem koraku je bila zapisana koda za izvajanje kinematičnih meritev v načinu LOCK, ki instrumentu omogoča, da samodejno spremlja reflektor na gibajočem vozičku. Vse meritve so bile shranjene neposredno v računalnik. V nadaljevanju sta podana zgradba in delovanje vmesnika GeoCOM. Programska oprema sistema TPS1200 Aplikacije on-board -Razvoj samostojnih programov -Krmiljenje: iz instrumenta Osnovni način Code script Starejši instrumenti Višji nivo GeoBASIC Geo C + + Aplikacije off-board -Razvoj samostojnih programov za dvostransko komunikacijo -Krmiljenje: iz zunanjega računalnika Osnovni način GSI ONLINE Višji nivo GeoCOM ASCII VisualBasic GeoCOM C++ GeoCOM TPS1000 TPS1100 TPS2000 TPS1200 Slika 4.4: Komunikacijski protokoli Leica Geo systems Figure 4.4: Leica Geo systems communication protocols je izdelana v sklopu posameznih senzorjev, zato so funkcije organizirane kot podsistemi slika 4.5 [GeoCOM Reference Manual V 1.20]. Slika 4.5: Pregled aplikacije odjemalec - strežnik [GeoCOM Reference Manual V 1.20] Figure 4.5: Overview Client-Server application [GeoCOM Reference Manual V 1.20] AUS Podsistem AltUser vključuje funkcije, dostopne preko gumba SHIFT+USER: nastavitev statusa ATE in LOCK načina merjenja AUT Avtomatizacija (angl. Automatization); modul, ki zagotavlja funkcije, kot so nadzor samodejnega prepoznavanja reflektorja, spremembe krožne lege ali določitve položaja. BAP Osnovne aplikacije (angl. Basic Applications); funkcije, ki jih je mogoče zelo enostavno uporabljati za pridobitev podatkov meritev. BMM Funkcije, ki nadzirajo osnovne vhodne ali izhodne funkcije, npr. nastavitev zvočnega alarma, itd. (angl. Basic Man Machine) COMF Komunikacija (angl. Communication); modul, ki nadzira osnovne komunikacijske parametre. Večina funkcij se nanaša tako na del za odjemalca kot na del za strežnik. COM Komunikacija (angl. Communication); funkcije, ki prav tako nadzirajo proces komunikacije. Te funkcije se nanašajo bodisi na del za odjemalca bodisi na del za strežnik. CSV Osrednje funkcije (angl. Central Services); ta modul zagotavlja funkcije za pridobivanje ali nastavljanje osrednjih/osnovnih informacij o instrumentu TP S1200. CTL Nadzor (angl. Control Task); modul vsebuje funkcije za spremljanje delovanja sistema. EDM Merjenje dolžine (angl. Electronic Distance Meter); modul za merjenje dolžine. MOT Motorizacija (angl. Motorization); funkcije za nadzor premika, tudi nadzor hitrosti, instrumenta. SUP Nadzor (angl. Supervisor); funkcije za nadzor nekaterih splošnih vrednosti instrumenta TPS1200: npr. stanje baterije, čas sistema, stanje pomnilne kartice. TMC Meritve in izračuni (angl. Theodolite Measurement and Calculation); ključni modul za pridobivanje podatkov o meritvah. GeoCOM je izveden kot: - sinhrona komunikacija: Komunikacija se izvaja med dvema enotama - odjemalcem (zunanja naprava) in strežnikom (instrument TPS1200). Posamezna komunikacijska enota je sestavljena iz zahteve in ustreznega odziva, kar pomeni, da para zahteva/odziv ne more prekiniti drugi par zahteva/odziv. - protokol Remote Procedure Call (EPC): Vsakemu postopku, kije izveden na oddaljenem instrumentu, je dodeljena identifikacijska številka. - Komunikacija se vzpostavi preko serijskega vmesnika RS-232. Za učinkovito komunikacijo morajo biti nastavljeni ustrezni komunikacijski parametri. - Programiranje temelji na konceptu DLL (Dynamic Link Library). V doktorski disertaciji je bila za pridobivanje meritev v načinu sledenja uporabljena funkcija GeoCOM Leica Geosystems VB_TMC_QuickDist(OnlyAngle, dSlopeDi-stance). Funkcija najprej prične z merjenjem dolžine v sledilnem načinu in čaka. dokler razdalja ni izmerjena. Nato sporoči vrednost horizontalne smeri, zenitne razdalje in poševne dolžine, ne pa tudi koordinat. Če razdalje ni mogoče izmeriti, instrument sporoči vrednost horizontalne smeri in zenitne razdalje ter ustrezno kodo sporočila. Tahimetri s samodejnim prepoznavanjem reflektorja Sodobni tahimetri omogočajo poleg merjenja horizontalnih smeri, zenitnih razdalj in dolžin še veliko drugih funkcij. V zadnjem desetletju so instrumenti opremljeni s funkcijami za samodejno prepoznavanje reflektorja, natančno viziranje. sledenje in samodejno izvajanje meritev. Z vključitvijo elektronske kamere v tahimeter je omogočeno samodejno odkrivanje statičnih in kinematičnih objektov. Razvoj geodetskih instrumentov je opisan na sliki 4.6. Kotna opazovanja Meritve dolžin Slika 4.6: Razvoj geodetskih instrumentov Figure 4.6: Development of geodetic instruments Pri izvajanju kinematičnih meritev, kjer je reflektor nenehno v gibanju, ni mogoče izvajati nadštevilnih opazovanj, s katerimi bi zmanjšali ali odstranili vplive delovnega okolja, instrumentalnih pogreškov ali negativne vplive matematičnega modela, ki so posledica po- enostavitev funkcionalnega modela. Zato je potrebno za vsak kinematični proces določiti naslednje vrednosti in izpolniti zahteve: - velikost vektorja premika; t.j. hitrost in smer premika reflektorja/objekta, - frekvenca meritev: Frekvenca meritev mora zadoščati opisu gibanja v smislu smeri, hitrosti in pospeška. Predvsem pri zajemu nepravilnega gibanja mora biti frekvenca meritev dovolj visoka. - samodejno sledenje premikajočega objekta - sinhronizacija: Visoka kakovost sinhronizacije senzorjev - senzorja za merjenje horizontalnih smeri in zenitnih razdalj, senzorja za merjenje razdalje in senzorja za samodejno sledenje - je ključnega pomena za natančne meritve. Poleg tega je potrebna velika natančnost pri določanju časa vsake meritve. - kalibrirani instrument: Pri kinematičnem merjenju nimamo možnosti izvajanja meritev horizontalne smeri in zenitne razdalje v dveh krožnih legah in nadštevilnih meritev dolžine, ki bi omogočale odstranitev instrumentalnih pogreškov. Zato je še posebej pomembno, da je instrument redno kalibriran. - kratek čas prenosa: Pri shranjevanju meritev/koordinat na zunanji računalnik je potrebno določiti čas prenosa. Izračunani položaj, ki mu pripada določen časovni trenutek, zaradi časa prenosa merskih podatkov ne ustreza popolnoma izmerjenemu položaju. Določitev oz. upoštevanje časa prenosa je pomembna zlasti v sistemih, ki delujejo v dejanskem času. - natančnost reflektorja, - možnost samodejnega iskanja reflektorja, kar je zlasti pomembno v primeru neenakomernega gibanja ali ovir med instrumentom in reflektorjem, - brezhibna komunikacija med merskim in operacijskim sistemom ter sistemom za obdelavo opazovanj, - možnost krmiljenja merskega sistema z zunanjim računalnikom, - možnost transformacij med koordinatnimi sistemi. V splošnem so zahteve pri kinematičnih meritvah veliko kompleksnejše, izbira ustreznega instrumenta pa je zelo pomembna. Samodejno prepoznavanje tarč in sledenje omogočata visoko produktivnost, krajši delovni čas, sledenje reflektorja pa je omogočeno tudi pri slabši vidljivosti ali celo v temi. Meritve z avtomatskim sledenjem lahko izvajamo na tri različne načine: - statični način: meritve se izvajajo proti statičnemu objektu/reflektorju, - način stop-and-go: instrument sledi reflektorju, meritve pa se sprožijo, ko objekt miruje, - kinematični način: med sledenjem se izvajajo tudi meritve dolžine, zenitne razdalje in horizontalne smeri. Zadnji način je tehnično najzahtevnejši; potrebna je visoka stopnja sinhronizacije senzorjev ter kratek čas meritev, zapisa le-teh in izračuna. 4.3 Funkcionalni in stohastični model: Diskretni model Wienerjevega procesa s pospeškom Z razvitim programom je bila vzpostavljena komunikacija z izbranim elektronskim tahi-metrom, kjer so bile meritve poševne dolžine d zenitne razdalje zr in horizontalne smeri hz zajete v kinematičnem načinu. Za vsako mersko epoho je bila registrirana tudi časovna komponenta. Podatki meritev so bili ovrednoteni s programom Matlab. Glavni poudarek dela je pridobiti oceno horizontalnega položaja (y,x) in visi ne z reflektorja na premičnem vozičku za vsak registrirani časovni trenutek tk z uporabo KF. Meritve je mogoče obdelati s KF na dva načina: - z razširjenim KF (angl. extended KF), kjer so v zanki KF obdelane direktne meritve - hz, d zr, ali - z linearnim KF (angl. linear KF), kjer so v zanki KF obdelane indirektne meritve -xyz V doktorski disertaciji je bil razvit in preizkušan linearni model KF. Za opis je bil uporabljen kinematični model. V osnovi so kinematični modeli za oceno stanja sistema izpeljani iz diferencialnih enačb. Definiramo jih tako, da je vrednost določenega odvoda položaja enaka nič. Za te vrste modele predpostavimo [Bar-Shalom in sod. (2001)]: - gibanje vzdolž izbrane koordinate je neodvisno od druge koordinate in - šumi v koordinatah so medsebojno neodvisni. Ker so v vsakem modelu prisotne določene motnje, lahko le-te obravnavamo in modeliramo kot slučajne pogreške. Eden od načinov modeliranja je modeliranje s časovno zveznim belim šumom sistema. Ker pa se opazovanja skoraj vedno izvajajo v diskretnem času, potrebujemo ustrezne enačbe stanja sistema za diskreten čas. V splošnem sta predstavljena dva kinematična modela stanja sistema [Bar-Shalom in sod. (2001)]: - kinematični model drugega reda ali model pospeška, opisanega z belim šumom (angl. white noise acceleration model): hitrost je Wienerjev proces, - kinematični model tretjega reda ali model pospeška, opisanega kot Wienerjev proces (angl. Wiener process acceleration model): pospešek je Wienerjev proces. Ker se direktna ali indirektna opazovanja zveznega stanja sistema izvajajo v diskretnem času, potrebujemo ustrezne kinematične modele za diskretni čas. Glede na način pridobitve diskretnih kinematičnih modelov s šumom (angl. discrete time kinematic models) le-te ločimo v dve skupini [Bar-Shalom in sod. (2001)]: - modeli, pridobljeni z diskretizacijo zveznih modelov, z zveznim belim šumom, za določeni interval vzorčenja, - modeli, pridobljeni z direktnim definiranjem šuma sistema v diskretnem času kot deloma konstanten slučajni beli proces - za šum sistema predpostavimo, daje konstanten za določeni interval vzorčenja in neodvisen med posameznimi intervali. V doktorski disertaciji smo zaradi visoke frekvence meritev predpostavili, da gibanje reflek- a intervalu dolžine At = 0.125 s, ki je enak časovni registraciji med dvema zaporednima trenutkoma pridobitve meritev. Na osnovi opisane predpostavke je bil razvit diskretni kinematični model tretjega reda ali diskretni model Wienerjevega procesa s pospeškom (DMWPP) ali model z odsekoma konstantnim pospeškom, opisanim kot Wienerjev proces (angl. discrete third-order kinematic model ali discrete Wiener process acceleration model ali piecewise constant Wiener process acceleration model), s tremi parametri za vsako smer. Prednost te vrste modelov je, da s šumom sistema lahko opišemo način gibanja, za katerega se predpostavljajo majhne spremembe pospeška znotraj časovnega intevala. Enak model je bil uporabljen za vse tri koordinate. Linearna enačba sistema in kovariančna matrika šuma sistema Stanje sistema xk zapišemo v vektorski obliki: Xk =[ x Vx ax y Vy ay z vz az]Tk , (4.1) kjer so x, y, z koordinate lege premikajočega se reflektorja, vx,vy, vz in ax, ay, az pa hitrost in pospešek v vseh treh smereh. Za simulirani kinematični proces predpostavimo naslednje: - časovno spreminjanje stanja sistema je neodvisno od nadzornih signalov: diskretna matrika vhoda v enačbi (2.36) je enaka nič, Gk = 0, - šum sistema vstopi v sistem preko matrike Tk, enačba (2.40). Enačba sistema diskretnega stanja sistema tretjega reda (pospešek se obravnava kot približno konstanten, kar pomeni, daje tretji odvod položaja enak nič) v primeru simuliranega kinematičnega procesa je: Xk+i = Fk • Xk + Tk • Wk, (4.2) kjer je wk nekorelirani šum sistema in predstavlja spremembo pospeška v intervalu k. Za spremembo pospeška predpostavimo, da je bela, s srednjo vrednostjo nič. Pospešek je tako opisan kot diskretni Wienerjev proces: wk - (0,a2w). (4.3) rk je v primeru simuliranega kinematičnega procesa vektor izboljšave, ki pomnoži skalami šum sistema wk. Ce wk predstavlja konstantni pospešek v intervalu h, je sprememba hitrosti v tem intervalu podana z wk • At. Vpliv pospeška na položaj je podan z: At2 Wk ^. Vektor izboljšave rk je za posamezno smer podan kot: rk = [ ^ At 1 f (4.4) in za vse tri smeri kot: rk =[ ^ At 1 ^ At 1 ^ At 1 f. (4.5) Matriko prehoda stanja Fk je mogoče izpeljati iz linearne diferencialne enačbe prvega reda zveznega stohastičnega sistema, enačba (2.27), kjer šum sistema vstopi v sistem preko matrike B: A • x + B • u. (4.6) Za naš primer kinematičnega sistema s tremi parametri neznank (položaj, hitrost, pospešek), s pospeškom, opisanim kot belim šumom, in indirektnimi meritvami prostorskega položaja - ima vektorska oblika, enačba (4.6), tri-dimenzionalnega problema obliko: Vx 0 1 0 0 0 0 0 0 0 x 0 ax 0 0 1 0 0 0 0 0 0 Vx 0 a x 0 0 0 0 0 0 0 0 0 Čx 1 Vy 0 0 0 0 1 0 0 0 0 y 0 ay = 0 0 0 0 0 1 0 0 0 Vy + 0 Čl y 0 0 0 0 0 0 0 0 0 Čy 1 Vz 0 0 0 0 0 0 0 1 0 z 0 Čz 0 0 0 0 0 0 0 0 1 Vz 0 iz _ 0 0 0 0 0 0 0 0 0 Čz 1 u(t). (4.7) V primeru, daje u(t) = 0 v enačbi (4.7), enačba opisuje gibanje točke s konstantnim pospeškom. Linearni diskretni model za simulirani kinematični proces pridobimo na osnovi A Fk sistema, podano v enačbi (4.2): 1 At (At)2 2 0 0 0 0 0 0 0 1 At 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 At (At)2 2 0 0 0 Fk (At) = 0 0 0 0 1 At 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 At (At)2 2 0 0 0 0 0 0 0 1 At 0 0 0 0 0 0 0 0 1 V zapisu matrike Fk (At) predpostavimo konstantni časovni interval med dvema meritvama, At = tk+1 — tk. Matrika Fk je konstantna ves čas procesa, zato lahko indeksiranje opustimo. Ce ponovno zapišemo wk = rk • Wk, enačba (2.40), in ob upoštevanju, da Wk v primeru simuliranega kinematičnega procesa predstavlja spremembo pospeška, wk ~ (0,pW), ki vstopi v sistem preko vektorja izboljšave rk, lahko zapišemo kovariančno matriko šuma sistema za eno dimenzijo v primeru DMWPP kot: Qk(i) = rk • e[wk • wTk] • rTk = rk • a2w • rTk = p«. At4 At3 4 2 A3 At2 At At2 2 At 1 (4.9) Ob upoštevanju, da so šumi med posameznimi smermi medsebojno neodvisni, lahko zapišemo kovariančno matriko šuma sistema za DMWPP kot: Qk At4 4 At3 2 At2 2 0 0 0 0 0 0 At3 2 At2 At 0 0 0 0 0 0 At2 2 At 1 0 0 0 0 0 0 0 0 0 At4 4 At3 At2 2 0 0 0 0 0 0 At3 2 At2 At 0 0 0 0 0 0 At2 2 At 1 0 0 0 0 0 0 0 0 0 At4 4 At3 2 A2t2 0 0 0 0 0 0 At3 2 At2 At 0 0 0 0 0 0 At2 2 At 1 (4.10) Ponovno lahko opustimo indeks k ker je matrika Q konstantna ves čas procesa. Rezultat različnih praktičnih primerov je interval za najboljšo izbiro velikosti šuma pospeška -0.5 • AaM < pw < AaM, kjer j e AaM največja sprememba pospeška [Bar-Shalom in sod. (2001)]. V primeru doktorske disertacije je bilo izvedenih več meritev. Velikost šuma pospeška je bila določena za vsak posamezen primer na osnovi povprečnega pospeška, brez prisotnosti grobih pogreškov. 2 Linearna enačba meritev in zakon o prenosu varianc in kovarianc V modelu linearnega KF so opazovanja v vektorju meritev zk, enačba (2.2), indirektna opazovanja x, y in z - položajne komponente reflektorja za vsak časovni trenutek tk■ V nadaljevanju so ta opazovanja imenovana meritve KF. Merska matrika Hk ima obliko: H 100000000 000100000 000000100 (4.11) Matrika Hk je konstantna ves čas procesa, zato lahko opustimo indeks k. Vrednosti meritev KF, x, y in z, so izračunane izven zanke KF, na osnovi direktnih meritev horizontalne smeri hz, zenitne razdalje zr in poševne dolžine d, z upoštevanjem desno-sučnega koordinatnega sistema v vsakem časovnem trenutku tk, kot: x zk[3 ,1] = y = z X_StPt + dk • sin(zrk) • cos(hzk) Y_StPt + dk • sin(zrk) • sin(hzk) Z_StPt + dk • cos(zrk) (4.12) Vrednosti X_StPt, Y_StSt, Z_StSt so koordinate stojišča instrumenta in so podane v preglednici 4.1. Zakon o prenosu varianc in kovarianc je potrebno uporabiti za izračun kovariančne matrike meritev Rk vektorja meirtev KF zk, ki je časovno odvisna. Kovariančna matrika Rl direktnih opažovanj hz, d in zr je definirana s standardnimi deviacijami opazovanj - 6 kaže na slabo pogojeno matriko. Za DMWPP lahko rečemo, da je pogojno število še vedno sprejemljivo, saj se njegova vrednost stabilizira pri 9. Na sliki 5.6 lahko vidimo, kako že pogojno število lahko pokaže odstopanja v modelu. Vrednost pogojnega števila je močno odvisna od vrednosti skalarja intenzitete šuma sis- Slika 5.5: Pogojno število k(P+), brez grobih pogreškov Figure 5.5: Condition number k(P+), no gross errors Slika 5.6: Pogojno število k(P+), z grobimi pogreški Figure 5.6: Condition number k(P+), with gross errors tema aw. Pogojno število zato lahko podaja razmerje med šumom sistema in šumom meritev. Na sliki 5.7 je predstavljeno pogojno število za različne vrednosti šuma sistema aw, aw = 0.001 : 0.020 : 0.51. Spodnji graf prikazuje pogojno število za najmanjšo vrednost aw, t.j. aw = 0.001. 10 9 6 5 50 100 150 200 250 300 350 400 450 Slika 5.7: Pogojno število k(P+) za različne vrednosti av Figure 5.7: Condition number for different values of aw 8 7 4 0 m V delu je bila izračunana tudi Josephova oblika a posteriori kovariančne matrike stanja sistema po enačbi (2.11), ki predstavlja stabilizirano obliko enačbe a posteriori kovariančne matrike. Z Josephovo posodobljeno kovariančno matriko ne dobimo bistveno boljših rezultatov. Za vse izračune smo tako uporabili kovariančno matriko stanja sistema, izračunano po enačbi (2.8). Vse lastne vrednosti so pozitivne realne vrednosti, glej slike 5.2, 5.3 in 5.4. 5.2.4 Lastnosti matrike K Vključitev meritev v oceno stanja sistema določa matrika Kalmanove izboljšave K, na katero vpliva negotovost v predhodni oceni stanja sistema in negotovost v oceni meritev. V splošnem lahko rečemo, da meritve z visoko natančnostjo zagotavljajo natančno oceno stanja sistema. Meritve namreč merijo stanje sistema - direktno ali indirektno - in tako lahko ob dovolj visoki natančnosti KF računa nanje kot na dobre pokazatelje dejanskega stanja sistema. Slednje bi morala odražati Kalmanova matrika izboljšave. Interpretacija in pomen matrike Kalmanove izboljšave Iz enačbe Kalmanove matrike izboljšave, enačba (2.6), vidimo, daje izboljšava: - sorazmerna z variancami napovedi stanja sistema v matriki P", - obratno sorazmerna z variancami izboljšav meritev, D"1 = (Hfc • P" • HT + R^)"1. Izboljšava je torej - velika, če je napoved stanja sistema netočna (ima veliko varianco) in meritev natančna (ima relativno majhno varianco), - majhna, če je napoved stanja točna (ima majhno varianco) in meritev slaba (ima relativno veliko varianco). Velika vrednost izboljšave kaže hiter odziv na meritve pri posodobitvi stanja sistema, medtem ko majhna vrednost izboljšave daje počasnejši odziv na meritve. Na slikah 5.8, 5.9 in 5.10 so prikazani grafi elementov matrike uteži K za komponente položaja, hitrosti in pospeška v vseh treh dimenzijah. 1.5 1 0.5 0 -0.5 1.5 1 0.5 0 -0.5 I 1.5 1 0.5 0 -0.5 Elementi matrike K za polozaj x IV- i- 1 - - -.—----—TV--- I1 0 0 0 50 50 50 100 100 100 150 200 250 300 Elementi matrike K za polozaj y 350 400 150 200 250 300 Elementi matrike K za polozaj z 350 400 150 200 250 m 300 350 400 450 450 450 KK komponento x merskih popravkov, zeleno - eleme nt i matrike K za komponento y merskih popravkov, modro - elementi matrike K za kompone nto z merskih popravkov: Položaj KK x of measurement correction, green - elements of matrix K for component y of measurement correction, blue - elements of matrix K for compon ent z of measurement correction: Position Elementi matrike K za hitrost vx m KK komponento x merskih popravkov, zeleno - eleme nt i matrike K za komponento y merskih popravkov, modro - elementi matrike K za kompone nt o z merskih popravkov: Hitrost KK x of measurement correction, green - elements of matrix K for component y of measurement correction, blue - elements of matrix K for compon ent z of measurement correction: Velocity K 1. Za vsako smer - x, y ali z - je v izračunu popravka (corr) položaja, hitrosti ali pospeška v vektorju stanja sistema pripadajoča utež popravkov meritev KF (zk — Hk • X-) - v smislu izbrane smeri vedno največja. x corr^i = K11 • (zfc — Hfc • X-)n + K12 • (zfc — Hfc • X-)2i + K13 • (zfc — Hfc • X-)3i, kjer je element Kn, ki vključuje komponento x popravka meritev KF, največji. Na Kx Kx K 5.11, lahko vidimo, daje za vsako komponento položaja utež ustreznega popravka meritev KF približno enaka 1 in za druga dva popravka meritev KF približno enaka 0. Enak trend je viden pri hitrosti, slika 5.12, in pospešku, slika 5.13. Elementi matrike K za pospešek ax Elementi matrike K za pospešek ay 100 -2 50- 0 50 100 150 200 250 300 350 400 450 Elementi matrike K za pospešek az 100 -21 50- (D — 0 - 0 50 100 150 200 250 300 350 400 450 m KK komponento x merskih popravkov, zeleno - eleme nt i matrike K za komponento y merskih popravkov, modro - elementi matrike K za kompone nt o z merskih popravkov: Pospešek KK x of measurement correction, green - elements of matrix K for component y of measurement correction, blue - elements of matrix K for compon ent z of measurement correction: Acceleration 3. Ce primerjamo uteži za položaj - x, y, z, hitrost - vx, vy, vz in pospešek - ax ay, az, vidimo, da so uteži za položaj velikostnega reda 10°, za hitrost velikostnega reda 101 in za pospešek velikostnega reda 102. Razlog za to je verjetno v tem, ker opazujemo samo položaje, zato jih je tudi mogoče bolje napovedati. Za vrednost skalarja intenzitete šuma sistema aw, kjer je dosežena ustrezna skladnost med napako ocenjevanja in standardnimi deviacijami položajev reflektorja, je napoved hitrosti in predvsem pospeška slabša, zato daje DMWPP večjo utež indirektnim informacijah, ki jih dobi iz meritev KF. K ritev IvF, enačba (4.16), in jih ni mogoče izračunati vnaprej. K sistema aw, predvsem elementi, ki določajo popravke hitrosti in pospeška. Na spodnjih slikah, 5.11, 5.12 in 5.13, so podani grafi za elemente matrike K za aw = 0.1. Na slikah 5.14, 5.15 in 5.16 so podani elementi matrike K za ekstremno vrednost aw = 0.000001. Elementi matrike K za polozaj, r-x komp.,g-y komp.,b-z komp. Slika 5.11: Elementi Kalmanove matrike izboljšave K; rdeče - elementi matrike K za komponento x merskih popravkov, zeleno - elementi matrike K za kompone nto y merskih popravkov, modro - elementi matrike K za kompone nto z merskih popravkov: Položaj, aw =0.1 Figure 5.11: Elements of Kalman gain matrix K; red - elements of matrix K for component x of Ky - elements of matrix K for component z of measurement correction: Position, aw =0.1 Elementi matrike K za hitrost, r-x komp.,g-y komp.,b-z komp. 15 10 0 f----—........... *----"" [ - h 0 50 100 150 200 250 300 350 400 450 m Slika 5.12: Elementi Kalmanove matrike izboljšave K; rdeče - elementi matrike K za kompone nto x Ky matrike K za komponento z merskih popravkov: Hitrost, aw = 0.1 Figure 5.12: Elements of Kalman gain matrix K; red - elements of matrix K for component x of Ky Elementi matrike K za pospešek, r-x komp.,g-y komp.,b-z komp. Slika 5.13: Elementi Kalmanove matrike izboljšave K; rdeče - elementi matrike K za kompone nt o x merskih popravkov, zeleno - elementi matrike K za komponento y merskih popravkov, modro - elementi matrike K za kompone nto z merskih popravkov: Pospešek, aw =0.1 Figure 5.13: Elements of Kalman gain matrix K; red - elements of matrix K for component x of Ky - elements of matrix K for compon ent z of measurement correction: Acceleration, aw = 0.1 Elementi matrike K za polozaj, r-x komp.,g-y komp.,b-z komp. Slika 5.14: Elementi Kalmanove matrike izboljšave K; rdeče - elementi matrike K za kompone nto x Ky matrike K za komponento z merskih popravkov: Položaj, aw = 0.000001 Figure 5.14: Elements of Kalman gain matrix K; red - elements of matrix K for compon ent x of Ky Elementi matrike K za hitrost, r-x komp.,g-y komp.,b-z komp. 200 250 m Slika 5.15: Elementi Kalmanove matrike izboljšave K; rdeče - elementi matrike K za komponento x merskih popravkov, zeleno - elementi matrike K za kompone nto y merskih popravkov, modro - elementi matrike K za komponento 2 merskih popravkov: Hitrost, aw = 0.000001 Figure 5.15: Elements of Kalman gain matrix K; red - elements of matrix K for component x of Ky - elements of matrix K for compon ent 2 of measurement correction: Velocity, aw = 0.000001 Elementi matrike K za pospešek, r-x komp.,g-y komp.,b-z komp. 0 50 100 300 350 400 450 50 100 150 200 250 300 350 400 450 m Slika 5.16: Elementi Kalmanove matrike izboljšave K; rdeče - elementi matrike K za kompone nto x Ky matrike K za kompone nto 2 merskih popravkov: Pospešek, aw = 0.000001 Figure 5.16: Elements of Kalman gain matrix K; red - elements of matrix K for component x of Ky Iz slike 5.15 in slike 5.16 lahko vidimo, da je pri majhni vrednosti šuma sistema, t.j. aw = 0.000001, zaupanje v model veliko, napovedane vrednosti in s tem tudi popravki K K za položajno komponento so nesprejemljivo razpršeni. Za vrednost aw = 0.000001 je DMWPP neskladen. V nadaljevanju je izvedena statistična ocena v domeni meritev in v domeni stanja sistema, z namenom, ugotoviti razloge za odstopanja med oceno napake in standardno deviacijo. 5.3 Statistični testi Enačbe za kovariančno matriko izboljšave, enačba (3.2), matriko Kalmanove izboljšave, enačba (2.6), izboljšavo meritev, enačba (2.12), a posteriori oceno stanja sistema, enačba (2.7), in pripadajočo a posteriori kovariančno matriko, enačba (2.8), so rezultat modela filtriranja in so natančne le, če držijo vse predpostavke o modeliranju, uporabljene pri izpeljavi filtra. V praksi temu ni vedno tako, zato je potrebno dodatno preizkusiti izračunane količine, ki podajajo informacijo o natančnosti ocenjenega stanja sistema [Bar-Shalom in sod. (2001)]. Pridobljeni rezultati, ki vključujejo oceno položaja, hitrosti in pospeška, so odvisni od uravnavanja razmerja med šumom sistema in šumom meritev ter a priori informacij. Z uravnavanjem filtra je mogoče opazovati ocenjeno stanje sistema in ali je dosežena skladnost z a posteriori kovariančno matriko P++. Skladnost preverimo na osnovi intervala a in sicer je skladnost dosežena, če je napaka stanja sistema znotraj intervala 2a pri 95% aa korenom ustreznih diagonalnih elementov v izračunani a posteriori kovariančni matriki P++ (standardna deviacija ai = i^P^^). Če napaka ocene presega interval 2a, rezultati KF niso sprejemljivi ob izbranih vhodnih parametrih. 5.3.1 Položaj reflektorja Glavni namen dela je bil oceniti položaj, slike 5.17, 5.18, 5.19, 5.20, 5.21 in 5.22, hitrost in pospešek reflektorja v vseh treh smereh, x, y in z. Za vse tri komponente lahko opazimo večjo razpršenost napovedanih vrednosti (zelena) in boljše ujemanje filtriranih vrednosti z referenčnimi položaji (rdeča), kar sledi tudi iz numeričnih vrednosti v preglednici 5.1. I).\I\VPP daje več zanesljivosti meritvam kot modelu oz. napovedanim vrednostim. 1 1 1 } - izsek ? I1 I1 i1 - - - -------- 1 SI' ° i1 - 0 50 100 150 200 250 300 350 400 450 m xx xx xx xx 205 210 215 220 225 230 235 240 m xx xx xx xx yy yy yy yy yy yy yy yy zz zz zz zz 205 210 215 220 225 230 235 240 m yy yy yy yy -12 -10 Y [m] Slika 5.23: Horizontalna ravnina Figure 5.23: Ground plane x 10-3 q|-1-1-1-1-1-1-1-1- 0 50 100 150 200 250 300 350 400 450 m Slika 5.24: Standardne deviacije položajnih komponent Figure 5.24: Standard deviation of position components Variance položajnih komponent stanja sistema Variance ocen stanja sistema, kijih izračunamo s IvF. bi morale v idealnem primeru ves čas vključevati dejansko stanje sistema. Grafično lahko to prikažemo z izrisom standardnih 1a \napovedano — merjeno\ [m] \ filtrirano — merjeno\ [m] x 3.785 0.935 y 3.590 0.597 z 1.031 0.095 Preglednica 5.1: Vsota absolutnih razlik med napovedanimi oz. filtriranimi in merjenimi vrednostmi 2a, ki ga izračuna KF, nam pove, koliko je verjetno, da je pravo stanje znotraj določene razdalje od ocenjenega stanja [Xegenborn (2003)]. KF je: - 68% zanesljiv, da je pravo stanje sistema znotraj intervala 1a; t.j. znotraj oddaljenosti 1a od ocenjenega stanja sistema, ali - 95% zanesljiv, da je pravo stanje sistema znotraj intervala 2a; t.j. znotraj oddaljenosti 2a od ocenjenega stanja sistema. Ce narišemo standardne deviacije a skupaj z napakami, mora napaka ostati znotraj izbranega intervala a Za DMWPP in aw = 0.1 sta interval 1a in interval 2a prikazana za vsako komponento položaja posebej na slikah 5.25, 5.26 in 5.27, numerične vrednosti pa so podane v preglednici 5.2. Xa grafičnem prikazuje izpisana številka koraka, v katerem 2a 200 250 m Slika 5.25: Napaka položaja x (polna rdeča črta), interval 1a (zelena polna črta) in interval 2a aw = 0.1 x a 2a aw = 0. 1 0 50 100 150 300 350 400 450 - ■ A i fti ffllA- v /M „ ^ h . fVI If 0 50 100 150 200 250 300 350 400 450 m Slika 5.26: Napaka položaja y (polna rdeča črta), interval 1a (zelena polna črta) in interval 2a (prekinjena zelena črta), brez prisotnosti grobih pogreškov, aw = 0.1 Figure 5.26: Position error y (red solid), a-interval (green solid) and 2a-interval (green dashed solid), no gross errors, aw = 0.1 170 0 -56- 375 373 _17 200 250 m x 10 15 10 5 58 0 0 50 100 300 350 400 450 Slika 5.27: Napaka položaja z (polna rdeča črta), interval 1a (zelena polna črta) in interval 2a (prekinjena zelena črta), brez prisotnosti grobih pogreškov, aw = 0.1 Figure 5.27: Position error z (red solid), a-interval (green solid) and 2a-interval (green dashed solid), no gross errors, aw = 0.1 Interval 1a Interval 2a (za normalno porazdelitev (za normalno porazdelitev 68% znotraj) 95% znotraj) x y z x y z % % znotraj % znotraj % znotraj % znotraj % znotraj °w = 0.1 91.4 91.8 92.1 95.5 95.5 95.5 &w = 0.1 log10(cond. P) — 8.115 brez grobih pogreškov &w = 0.06 87.3 87.3 87.8 92.1 92.1 92.1 °w = 0.06 logw(condP) — 7.782 brez grobih pogreškov &w = 0.01 45.1 46.5 46.3 61.0 60.3 60.8 &w = 0.01 logio(condP) — 6.642 brez grobih pogreškov Preglednica 5.2: Odstotek filtriranih vrednosti znotraj intervala 1a in intervala 2a Iz preglednice 5.2 lahko opazimo zmanjšanje pogojnega števila (hkrati nastopi boljša konvergenca sledi matrike P^) z večjim zaupanjem v model (manjša vrednost aw). Hkrati pa je veliko manj napak položaja znotraj intervala 1a ali intervala 2a ocenjenih položajev, kar predstavlja neskladje z dejanskim stanjem. Za aw = 0.1 lahko sklepamo, da so ocene stanja sistema skladne z dejanskimi napakami. Rezultati za vrednost aw = 0.06 so še sprejemljivi za interval 2a in povsem sprejemljivi za interval 1a. Za aw = 0.01 lahko rečemo, da so ocene stanja položaja neskladne z izračunanimi standardnimi deviacijami. Samo približno 46% ocenjenih položajev je znotraj 1a-intervala in 61% ocenjenih položajev znotraj intervala 2a. V praktičnih primerih modelirani šum sistema ne odraža vedno pravega šuma sistema. Zato je težko pridobiti natančne vrednosti varianc komponent stanja sistema. Prav tako se lahko šum spreminja skozi proces. Zato moramo v vsakem modelu preučiti vpliv večjega oz. manjšega šuma sistema na rezultate KF. V idealnem primeru modelirani šum sistema popolnoma odraža šum sistema dejanskega procesa. KF uporabi modelirani šum sistema kot izraz nezanesljivosti v oceni stanja sistema zaradi nemodeliranih vplivov. Če je modelirani šum sistema večji kot dejanski šum sistema, potem KF predpostavi, da je v sistemu prisotno več šuma in zato po nepotrebnem preveč zmanjša zanesljivost ocene stanja sistema. V primeru določanja položaja pomeni modeliranje prevelikega šuma večjo nezanesljivost v ocenjenem položaju in odločanje kontrolorja v tem primeru je preveč pazljivo. Na voljo ima slabe informacije o oceni položaja. Večji šum sistema torej pomeni poslabšanje napovedanega stanja sistema, kar ima za posledico vključitev meritev [Ne-genborn (2003)]. Opazimo lahko tudi različne velikostne rede standardnih deviacij za različne položaj ne komponente, ki pa so znotraj izbranega intervala. Če bi v modelu DMWPP dali preveliko zaupanje napovedanim vrednostim, bi standardne deviacije za vse tri ocenjene položajne 0 ložaju boljša, ni pa skladna z dejansko napako. V primeru, da imamo opraviti z nižjo zanesljivostjo sistema, se lahko zgodi, da prava vrednost ni znotraj intervala zaupanja. Hitrost in pospešek reflektorja in pripadajoče standardne deviacije so podane v poglavju 5.4. 5.3.2 Skladnost DMWPP v domeni meritev Izboljšave meritev KF na sliki 5.29, s pripadajočo kovariančno matriko, izračunano v postopku filtriranja, so preizkušane na osnovi normirane kvadratne izboljšave, enačba (3.5), in razmerjem verjetnosti, podanim v enačbi (3.6). Za DMWPP je vrednost prostostne stopnje enaka r = 3. Statistični testi so izvedeni z a tveganje a = 0.05; t.j. P < x3, i_o.o5|Ho} > 1 — 0.05. Pri 95% stopnji zaupanja, za vektor meritev KF, ki vključuje tri komponente, je območje zaupanja za enostranski test definirano z x3,1-0.05 = 7.815. Normirana kvadratna izboljšava za DMWPP je grafično prikazana na sliki 5.28. Če filter ustrezno opisuje sistem, se mora ocenjena napaka nahajati znotraj intervala 2a 95% vseh merskih korakov [Gao in sod. (2005)]. 300 250 200 § 150 100 50 0 0 50 100 150 200 250 300 350 400 450 m Slika 5.28: Normirana kvadratna izboljšava Figure 5.28: Normalized innovation squared Za izboljšave meritev IvF, enačba (4.12). ki predstavljajo popravke meritev, se predpostavlja, da jih je mogoče opisati z normalno porazdelitvijo, s srednjo vrednostjo nič in pripadajočo kovariančno matriko Dk, dano v enačbi (3.2): dk -N (0, Dk). (5.4) Za napake se predpostavlja, da so normalno porazdeljene, saj normalna porazdelitev običajno dobro aproksimira porazdelitev merjenih količin. Za normalno porazdelitev je značilno, da so slučajne napake večjih vrednosti redke. Rezultati statističnih testov, kot so intervali zaupanja in intervali napovedovanja, zahtevajo za njihovo veljavnost normalno porazdelitev. Če je srednja vrednost napak enaka nič, potem so napake popolnoma slučajne. Če pa srednja vrednost napak ni enaka nič, se lahko zgodi, da model ni ustrezen za opis sistema ali pa da napake niso povsem slučajne in so prisotni sistematični pogre-ški, [Matlab Help]. Za naš model DMWPP je porazdelitev izboljšav meritev KF podana na slikah 5.30, 5.31 in 5.32. 0.1 0 -0.1 -0.2 0.2 0.1 0 -0.1 -0.2 0.01 0 -0.01 -0.02 -0.03 1 -1 cm|' 50 50 100 150 200 250 m 300 350 400 450 100 150 200 250 m 300 350 400 450 1 1 cm 1 i i i i i i '-pr/fJ IfMHHI' jw^tHM—^HHM \r-- -1 cm i i i i i i i P i 0 50 100 150 200 250 m 300 350 400 450 0 0 Slika 5.29: Izboljšave meritev KF Figure 5.29: KF-measurement innovation Razpršenost izboljšave x -0.6 -0.5 0.1 0 0.1 srednja vrednost Slika 5.30: Histogram izboljšave meritev KF za komponento x Figure 5.30: Innovations histogram for component x Razpršenost izboljšave y -0.1 0 0.1 srednja vrednost Slika 5.31: Histogram izboljšave meritev KF za komponento y Figure 5.31: Innovations histogram for component y Razpršenost izboljšave z i ilillmlliii i -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 srednja vrednost Slika 5.32: Histogram izboljšave meritev KF za komponento z Figure 5.32: Innovations histogram for component z 70 60 - 50 40 30 - 20 5.3.3 Skladnost DMWPP v domeni stanja sistema Drugi test skladnosti se izvede na osnovi normirane napake stanja sistema, enačba (3.9), kjer je razmerje verjetnosti podano v enačbi (3.10). Posamezni testi so bili izvedeni za različne vrednosti intenzitete skalarja šuma sistema aw, z namenom preveriti skladnost modela v domeni stanja sistema. V primeru vrednosti intenzitete šuma sistema aw = 0.1 je normirana napaka stanja sistema podana na sliki 5.33. Ob tveganju a = 5% in prosto-stni stopnji 9, ki ustreza dolžini vektorja stanja sistema, znaša kritična vrednost testne statistike, ki je porazdeljena po porazdelitvi X■> 16.919. V postopku filtriranja testna statistika preseže kritično vrednost za 26 točk o d N = 441 meritev, kar predstavlja 6% vseh točk. Situacijo lahko opišemo kot sprejemljivo. Testna statistika za preizkušanje skladnosti v domeni stanja sistema je bila izračunana po enačbi (3.9), kjer je a priori varianca izračunana na osnovi predhodnega filtriranja in znaša a0 = 0.01. Xa slikah 5.34, 5.35 in 5.36 so prikazani histogrami popravkov komponent stanja sistema. Slika 5.33: Normirana napaka stanja sistema Figure 5.33: Normalized system state error -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 srednja vrednost izboljšava z 0 0.01 srednja vrednost 80 60 40 20 Slika 5.34: Histogram popravkov komponent stanja sistema - položajne komponente Figure 5.34: System state incremental update histograms - position components izboljšava vx ......J L..... -2.5 -2 150 100 -4 -3 80 -1.5 -1 -0.5 0 0.5 1 srednja vrednost srednja vrednost izboljšava vz -0.3 -0.2 -0.1 0 0.1 srednja vrednost Slika 5.35: Histogram popravkov komponent stanja sistema - komponente hitrosti Figure 5.35: System state incremental update histograms - velocity components 80 r iz CO 60 C d ? 40 (D 20 0 -10 200 r iz CO 150 - C ? 100 - cu 50 - boljšava ax 0 srednja vrednost srednja vrednost 40 boljšava az -2.5 -2 -1.5 -1 -0.5 0 0.5 srednja vrednost 1 1.5 2 2.5 Slika 5.36: Histogram popravkov komponent stanja sistema - komponente pospeška Figure 5.36: System state incremental update histograms - acceleration components 100 50 2 5 10 15 20 5.4 Referenčna tirnica kot parameter ocenjevanja Proces uravnavanja razmerja med šumom sistema in šumom meritev, ki sta podana z matrikama Q in R, je nujen za pridobitev dobrih rezultatov modela ter skladnosti med ocenami stanja sistema in dejanskimi napakami. Xa osnovi simulacij različnih pogojev kinematičnega procesa (različne postavitve instrumenta, različne hitrosti vozička, različne vrednosti skalarja intenzitete šuma sistema) je bil razvit model, uporabljen za pridobitev ocene stanja sistema. Razen v primeru, ko so poznane prave vrednosti komponent stanja sistema, je zelo težko ovrednotiti, kako natančno in zanesljivo nam model KF izračuna vrednosti neznank, ki nas zanimajo. V ta namen je bila uporabljena referenčna tirnica, na osnovi katere so bili pridobljeni neodvisni referenčni parametri za oceno modela. Xa sliki 5.37 je grafično prikazana pravokotna oddaljenost merjenih in filtriranih/ocenjenih položajev od referenčne tirnice v horizontalni ravnini. Izrisan je tudi interval 1 cm, katerega vrednost je določena glede na pričakovano natančnost meritev, kije odvisna predvsem od uporabljene funkcije Leica, Geosystems GeoCOM, t.j. VB_TMC_QuickDist(OnlyAngle, dSlopeDistance). in hitrosti vozička. Vsi horizontalni položaji (y,x)k, določeni na osnovi direktnih meritev, ležijo znotraj itervala 1 cm Slika 5.37: Pravokotna oddaljenost filtriranih (rdeče) in merjenih (modra) vrednosti od 1 cm Figure 5.37: Perpendicular distance of filtered (red) and measured (blue) values from 1 cm Na sliki 5.38 je podana pravokotna oddaljenost z-komponente od referenčne tirnice. 0.015 S 0.01 C O C 0 C £ 1 0.005 "D O O R 0 E o 1, -0.005 CO ■O ■O 0 CO c 1 -0.01 > ci -0.015 0 50 100 150 200 250 300 350 400 m Slika 5.38: Pravokotna oddaljenost filtriranih (rdeče) in merjenih (modro) vrednosti z-komponente od referenčne tirnice in interval 1 cm Figure 5.38: Perpendicular distance of filtered (red) and measured (blue) values of z-component from the reference trajectory and 1 cm-bound + 1cm 374 modro: merjene vrednosti - rdece: filtrirane vrednosti —d_rt, 1 -1cm Na osnovi podatkov referenčne tirnice in časovnega intervala med meritvami, At = 0.125 s, sta bila izračunana referenčna hitrost in pospešek v vseh treh dimenzijah. Xa slikah 5.39 in 5.40 sta podana hitrost in pospešek za vse tri smeri, s pripadajočim intervalom 2a. Xa sliki 5.41 in sliki 5.42 so grafično podane napake hitrosti in pospeška v vseh treh smereh s pripadajočim intervalom a in intervalom 2a. Opazimo lahko veliko večje odstopanje napak hitrosti in napak pospeška kot v primeru položajev, podpoglavje 5.3.1, slike 5.25, 5.26 in 5.27. Razlog za to je v manjši opazovanosti komponent hitrosti in pospeška, ki jih ne merimo direktno, zato imamo na voljo manj informacij o teh komponentah stanja sistema. Hitrost in pospešek sta tako ocenjena le na osnovi modela. Prav tako je vzrok v nestabilnosti numeričnega odvajanja. Standardne deviacije hitrosti in pospeška v vseh treh dimenzijah so grafično prikazane na slikah 5.43 in 5.44. Iz slik 5.43, 5.44 in 5.24 lahko opazimo različne velikostne rede standardnih deviacij za hitrost, pospešek in položajne komponente, kakor tudi med komponentami za različne smeri. Posledica različnih velikostnih redov standardnih deviacij je slabša vrednost pogojnega števila. 0.2 „ 0.1 o 0)