Prostor kraj čas 15 GEOMETRIČNI POPRAVKI OPTIČNIH SATELITSKIH POSNETKOV Aleš Marsetič, Krištof Oštir in Mojca Kosmatin Fras PROSTOR, KRAJ, ČAS PROSTOR, KRAJ, ČAS 15 GEOMETRIČNI POPRAVKI OPTIČNIH SATELITSKIH POSNETKOV Aleš Marsetič, Krištof Oštir in Mojca Kosmatin Fras Uredil: Žiga Kokalj Izdajatelj: Inštitut za antropološke in prostorske študije, ZRC SAZU Za izdajatelja: Ivan Šprajc Založnik: Založba ZRC Za založnika: Oto Luthar Recenzenti: Mihaela Triglav Čekada, Andreja Gomboc, Dejan Grigillo in Domen Mongus Prva e-izdaja Ljubljana 2018. Knjiga je izšla s podporo Javne agencije za raziskovalno dejavnost Republike Slovenije, tudi prek razpisa za znanstvene monografije v letu 2017. Kataložni zapis o publikaciji (CIP) so pripravili v Narodni in univerzitetni knjižnici v Ljubljani. COBISS.SI-ID=292901120 ISBN 978-961-05-0039-1 (pdf) Način dostopa (URL): http://zalozba.zrc-sazu.si/p/P15 Vse pravice pridržane. Noben del te izdaje ne sme biti reproduciran, shranjen ali prepisan v kateri koli obliki oz. na kateri koli način, bodisi elektronsko, mehansko, s fotokopiranjem, snemanjem ali kako drugače, brez predhodnega pisnega dovoljenja lastnikov avtorskih pravic . © 2018 avtorji, Založba ZRC, ZRC SAZU. Naslovnica: © ESA/Belspo, mozaik izdelal VITO. PROSTOR, KRAJ, ČAS 15 GEOMETRIČNI POPRAVKI OPTIČNIH SATELITSKIH POSNETKOV Aleš Marsetič, Krištof Oštir in Mojca Kosmatin Fras Ljubljana 2018 Izvleček V knjigi je opisan popolnoma samodejen postopek geometričnih popravkov (ortorektifikacije) optičnih satelitskih posnetkov, ki deluje brez posredovanja operaterja med obdelavo. Dobljeni ortorektificiran posnetek je v državnem koordinatnem sistemu in predstavlja primerno podlago za prostorske analize, ki se uporabljajo pri številnih aplikacijah, na primer pri opazovanju naravnih nesreč, nadzorovanju aktivnosti na morju, spremljanju vegetacije in podobno. V knjigi je detajlno dokumentiran izviren postopek, ki povezuje več različnih metod v enoten in robusten sistem za samodejno izdelavo satelitskega ortofota. Opisan je celoten postopek samodejne ortorektifikacije, ki ga sestavljajo štirje osnovni moduli: modul za branje in pripravo metapodatkov, modul za samodejno določanje oslonilnih točk, modul za izračun parametrov geometričnega modela in modul za izvedbo ortorektifikacije. Osrednji del raziskave je modul za izračun parametrov geometričnega modela, ki je tudi najbolj zapleten del postopka. Modul deluje iterativno na več ravneh in z izravnavo po metodi najmanjših kvadratov ter dvema zaporednima metodama (RANSAC in robustna ocena), ki odstranita oziroma zmanjšata vpliv grobih napak pri izračunu parametrov notranje in zunanje orientacije posnetka. Opravljeno je bilo tudi testiranje modulov z visokoločljivimi posnetki RapidEye in zelo visokoločljivimi posnetki WorldView-2. Posnetki so bili posneti na različnih območjih Slovenije in obsegajo tako ravninska območja, kot gričevja, hribovja in gorovja. V knjigi predstavljeni poskusi so omogočili ovrednotenje delovanja samodejnega postopka določanja točk, geometričnega modela, izločanja grobih napak in oceno položajne točnosti ortofotov. Analiza rezultatov testov kaže, da lahko samodejni postopek izdela ortofote s položajno točnostjo okrog dveh pikslov ali manj, tudi če je med samodejno določenimi oslonilnimi točkami prisotnih več grobih napak. KLJUČNE BESEDE: samodejna ortorektifikacija, geometrični popravki, optični satelitski posnetki, geometrični model, RANSAC, robustna ocena 4 Abstract The monograph presents a fully automated procedure for geometric corrections (orthorectification) of optical satellite images, without any manual intervention during the processing. The resulting orthorectified image is in the national coordinate system and constitutes a suitable source for spatial analyses that are used routinely in many applications, for example, in monitoring of natural disasters and sea activities, observation of vegetation, etc. In the monograph the original procedure that connects several different methods into a single, robust system for automatic generation of orthoimages is documented in detail. The work describes the whole automatic orthorectification process, which comprises four basic modules: a module for extracting and preparing the metadata, a module for automatic extraction of ground control points, a module for calculation of parameters of the geometric model, and a module for orthorectification. The central part of the research is the module for the calculation of parameters of the geometric model, which is the most complicated part of the process. The module operates iteratively on several levels. It uses the least squares adjustment method and two consecutive methods (RANSAC, robust estimation) for removal and/or reduction of the impact of gross errors in the calculation of the interior and exterior orientation parameters of the image. Experiments with high-resolution RapidEye and very high-resolution WorldView-2 images were also carried out to test the modules. The images were acquired over different regions of Slovenia and include flat, hilly and mountainous areas. The experiments presented in the monograph evaluate the procedure for automatic extraction of points, the geometric model, the elimination of gross errors, and the positional accuracy of the orthoimages. The analisys of the results indicate that the automated procedure produces orthoimages with a positional accuracy of about two pixels or better, even if several gross errors are present among the automatically extracted ground control points. KEYWORDS: automatic orthorectification, geometric correction, optical satellite images, geometric model, RANSAC, robust estimation 5 Vsebina 1 UVOD ..................................................................................................................... 12 1.1 Pregled obstoječih raziskav na področju samodejne ortorektifikacije .................................................................................................. 15 1.1.1 Geometrični model ................................................................................. 15 1.1.2 Slikovna korelacija .................................................................................. 17 1.1.3 Samodejno georeferenciranje ............................................................ 18 1.2 Struktura knjige ................................................................................................. 19 2 IZHODIŠČA IN TEMELJNI POJMI .................................................................... 21 2.1 Geometrija zajema posnetkov glede na vrsto senzorja ....................... 21 2.2 Ravni obdelave posnetkov ............................................................................. 23 2.3 Ortorektifikacija ................................................................................................. 24 2.4 Samodejna procesna veriga .......................................................................... 25 3 UPORABLJENI PODATKI................................................................................... 29 3.1 RapidEye ............................................................................................................... 29 3.2 WorldView-2 ....................................................................................................... 32 4 SAMODEJNO DOLOČANJE OSLONILNIH TOČK ....................................... 35 4.1 Postopek samodejnega določanja oslonilnih točk ................................ 36 4.2 Implementacija .................................................................................................. 42 4.2.1 Modul za branje in pripravo metapodatkov .................................. 43 4.2.2 Modul za samodejno določanje oslonilnih točk ........................... 49 4.3 Testiranje .............................................................................................................. 50 4.3.1 Ocena točnosti samodejno določenih oslonilnih točk ............... 50 5 GEOMETRIČNI MODELI SENZORJEV ............................................................ 55 5.1 Senzorji na satelitih za opazovanje Zemlje .............................................. 56 5.1.1 Točkovni senzorji ..................................................................................... 56 5.1.2 Vrstični senzorji ........................................................................................ 57 5.1.3 Polnoslikovni senzorji ............................................................................ 59 5.2 Popravki koordinat oslonilnih točk ............................................................. 61 5.2.1 Atmosferska refrakcija ........................................................................... 61 5.2.2 Ukrivljenost Zemlje ................................................................................. 63 5.3 Geometrični model vrstičnega senzorja .................................................... 64 6 5.3.1 Izravnava po metodi najmanjših kvadratov ................................... 64 5.3.2 Funkcionalni model ................................................................................ 66 5.3.3 Stohastični model ................................................................................... 74 5.3.4 Reševanje sistema linearnih enačb ................................................... 75 5.3.5 Izločanje grobih napak .......................................................................... 76 5.4 Geometrični model polnoslikovnega senzorja ....................................... 81 5.4.1 Funkcionalni model ................................................................................ 81 5.4.2 Stohastični model ................................................................................... 85 5.4.3 Reševanje sistema linearnih enačb in izločanje grobih napak ........................................................................................................... 85 5.5 Implementacija .................................................................................................. 85 5.5.1 Modul za izračun parametrov geometričnega modela .............. 86 5.6 Testiranje .............................................................................................................. 89 5.6.1 Vpliv števila oslonilnih točk na izračun parametrov geometričnega modela ........................................................................ 89 5.6.2 Metode za robustno oceno ................................................................. 94 5.6.3 Detekcija in izločanje grobih napak .................................................. 97 6 ORTOREKTIFIKACIJA ....................................................................................... 105 6.1 Metode izdelave ortofotov ......................................................................... 105 6.1.1 Direktne metode ................................................................................... 106 6.1.2 Indirektne metode ................................................................................ 111 6.2 Postopek izdelave ortofotov ...................................................................... 113 6.3 Implementacija ............................................................................................... 115 6.3.1 Modul za izvedbo ortorektifikacije .................................................. 115 6.4 Testiranje ........................................................................................................... 117 6.4.1 Ocena položajne točnosti ortofotov ............................................... 117 7 ZAKLJUČKI ......................................................................................................... 122 8 VIRI IN LITERATURA ......................................................................................... 128 Seznam slik ..................................................................................................................... 136 Seznam preglednic ...................................................................................................... 139 Priloga ............................................................................................................................... 140 7 Slovar asinhrono snemanje (angl. asinchronous imaging) – Satelitsko snemanje z vzdolžnim senzorjem pri čemer je hitrost satelita večja od hitrosti snemanja posnetka. atmosferska višina homogene gostote (angl. homogeneous atmospheric density scale height) – Višina plasti ozračja s homogeno gostoto. bralnik podatkov (angl. parser) – Program za branje in izločanje podatkov iz datotek v standardnih zapisih (npr. XML). čas zaznave (angl. integration time) – Čas, v katerem slikovni elementi detektorja zbirajo naboj med prenosi do transportnih registrov. divergenca (angl. divergence) – Razhajanje vrednosti neznanke v sosledju iteracij. geometrična popačenja (angl. geometric distortions) – Popačenja posnetka, ki so posledica razgibanosti zemeljskega površja, ukrivljenosti Zemlje, napak senzorja in kota gledanja senzorja. iterativna bikonjugirana gradientna metoda (angl. iterative biconjugate gradient method, IBGM) – Iterativna metoda za reševanje sistemov linearnih enačb, ki temelji na uporabi gradienta. iterativna fotogrametrična metoda (angl. iterative photogrammetric method) – Iterativna direktna metoda za izračun prostorskih koordinat iz slikovnih koordinat, ki temelji na izračunu preseka slikovnega žarka z digitalnim modelom reliefa. iterativna metoda sledenja žarku (angl. iterative ray tracing method) – Iterativna direktna metoda za izračun prostorskih koordinat iz slikovnih koordinat, ki temelji na izračunu točke preseka z iterativnim podaljševanjem žarka gledanja, dokler ne doseže površja. kalibracija v tirnici (angl. in-flight calibration) – Kalibracija elementov detektorja z analizo posnetkov, ki so bili zajeti na kalibracijskih območjih. 8 kalibriran surov posnetek (angl. sensor-corrected image) – Posnetek, ki je bil obdelan do ravni, kjer so upoštevani le radiometrični popravki zaradi variacije elementov znotraj detektorja (ponavadi je to raven 1A). koren srednje kvadratne napake (angl. root mean square error, RMSE) – Kvadratni koren aritmetične sredine kvadratov odstopanj. Mera za odstopanje opazovanih vrednosti od pričakovanih vrednosti. kvaternion (angl. quaternion) – Matematičen zapis usmeritev in rotacij objekta v trirazsežnem prostoru v obliki kompleksnega števila. model afine transformacije (angl. affine transformation model) – Model afine transformacije je poseben primer modela racionalnih funkcij z osmimi parametri. model direktne linearne transformacije (angl. direct linear transformation, DLT) – Model DLT je poseben primer modela racionalnih funkcij z enajstimi parametri. model racionalnih funkcij (angl. rational function model, RFM) – Empirični model, ki nadomešča rigorozne modele in je v obliki racionalne funkcije z 80-imi parametri. način upočasnjevanja (angl. slow-down mode) – Zmanjšanje hitrosti vzdolžnega snemanja glede na hitrost leta satelita, s katerim lahko povečamo čas zaznave. nagib, naklon, zasuk (angl. roll, pitch, yaw) – Rotacije platforme glede na orbitalni koordinatni sistem (tirnico letenja). na Zemljo pritrjen geocentrični koordinatni sistem (angl. Earth Centered Earth Fixed, ECEF) – Ortogonalni koordinatni sistem z izhodiščem v centru Zemlje. objektni prostor (angl. object space) – Trirazsežni prostor, v katerem so podane prostorske koordinate oslonilnih točk. obrnjeno projiciranje posamičnega žarka (angl. single-ray backprojection) – Iterativne direktne metode za izračun prostorskih koordinat iz slikovnih koordinat in z uporabo digitalnega modela reliefa. 9 odsekovne polinomske funkcije (angl. piecewise polynomial functions) – Polinomske funkcije, sestavljene iz več medsebojno povezanih odsekov. podatkovna struktura (angl. data structure) – V jeziku IDL je to spremenljivka, ki lahko vsebuje različne podatkovne tipe. polnoslikovni snemalni sistem (angl. full-frame camera) – Snemalni sistem, ki naenkrat zajame celoten posnetek. popolni ortofoto (angl. true orthophoto) – Ortofoto, ki je popravljen tako za vplive geometričnih popačenj kot za vpliv visokih objektov, ki jih digitalni model reliefa ne vsebuje (npr. stavbe). prečni snemalni sistem (angl. whiskbroom scanner) – Snemalni sistem, ki površje opazuje s premikanjem vidnega polja senzorja pravokotno (prečno) na smer gibanja platforme. prečno na orbito (angl. across track) – Pravokotno na smer gibanja platforme. prostorski presek (angl. space intersection) – Postopek določevanja prostorskih koordinat točk, ki se pojavljajo na prekrivajočih se delih dveh ali več posnetkov z znano notranjo in zunanjo orientacijo. V matematičnem smislu je to izračun preseka dveh ali več prostorskih premic z izravnavo. racionalni polinomski koeficienti (angl. rational polynomial coefficients, RPC) – Koeficienti, ki se uporabijo v modelu racionalnih funkcij za pretvorbo slikovnih koordinat v prostorske koordinate ali obratno. ravni obdelave posnetkov (angl. image processing levels) – Ravni obdelave posnetkov, ki jih nudijo ponudniki posnetkov. Ravni si sledijo od surovih posnetkov do ortorektificiranih in mozaičenih posnetkov. razmerje signal-šum (angl. signal-to-noise ratio, SNR) – Vrednost, ki podaja stopnjo signala (jakost) glede na neželene pojave (šum). razpršena matrika (angl. sparse matrix) – Matrika, v kateri ima večina celic vrednost 0. reflektorji (angl. reflecting telescopes) – Teleskopi, ki za objektive uporabljajo zrcala. refraktorji (angl. refracting telescopes) – Teleskopi, ki za objektive uporabljajo leče. 10 robustna ocena (angl. robust estimation) – Metoda izravnave, ki temelji na zmanjševanju vpliva opazovanj z napakami z uporabo utežne funkcije. skoraj realni čas (angl. near-real-time) – Krajši časovni zamik, ki nastane pri samodejni obdelavi podatkov med zajemom izvornega podatka in dostavo izdelanega produkta. sledilec zvezd (angl. star tracker) – Naprava, ki določa rotacije platforme na podlagi sledenja položajev zvezd. slikovno ujemanje z metodo najmanjših kvadratov (angl. least squares matching, LSM) – Metoda slikovnega ujemanja za iskanje homolognih točk, ki temelji na minimiziranju razlik med radiometričnimi vrednostmi dveh slik z metodo najmanjših kvadratov. tehnologija TDI (angl. Time Delay and Integration) – Metoda vzdolžnega skeniranja, ki z večkratnim zajemom istih objektov zagotavlja dobro kakovost posnetka s hitrejšim zajemom vrstic. vzdolžni snemalni sistem (angl. pushbroom scanner) – Snemalni sistem, ki naenkrat zajame celotno vrstico posnetka v smeri gibanja platforme. vzdolž orbite (angl. along track) – V smeri gibanja platforme. zakrivanje (angl. occlusion) – Prikriti objekti ali deli objektov, ki na posnetku niso vidni zaradi kota snemanja in/ali geometričnih popačenj. značka (angl. tag) – Sestavina označevalnih jezikov, npr. HTML, SGML, XML. žarek gledanja (angl. view ray) – Žarek (linija), ki povezuje projekcijski center posnetka ter isto točko na posnetku in v prostoru. 11 1 UVOD V zadnjih desetletjih je bilo izstreljenih več satelitov za opazovanje Zemlje za potrebe različnih aplikacij (NSSDC 2017; CEOS 2017). Količina satelitskih posnetkov hitro narašča, njihova uporaba je prisotna na vedno več področjih, dostopni pa so tako znanstvenikom kot podjetnikom. V večini primerov uporabniki zahtevajo že predobdelane podatke oz. posnetke brez geometričnih popačenj (angl. geometric distortions), ki so posledica razgibanosti zemeljskega površja, ukrivljenosti Zemlje, napak senzorja in kota gledanja senzorja, in ki so v izbranem koordinatnem sistemu. Poleg tega v nekaterih aplikacijah potrebujemo geometrično popravljene posnetke v skoraj realnem času (angl. near-real-time). Na primer, pri opazovanju naravnih nesreč sta hiter zajem podatkov in njihova takojšnja analiza ključnega pomena (Veljanovski in dr. 2012; Joyce in dr. 2009). Tako rekoč vse aplikacije, ki temeljijo na opazovanju Zemlje, bi pridobile na uporabnosti, če bi bil postopek od zajema podatkov do dostave rezultatov končnemu uporabniku opravljen hitro in samodejno. Dosedanje metode obdelave satelitskih posnetkov temeljijo na postopkih, ki so dokaj samostojni, vendar nepovezani v celotni verigi od satelitskega snemanja do končnega izdelka (tj. rektificiran posnetek oz. karta). Želja tako raziskovalcev kot končnih uporabnikov satelitskih posnetkov je avtomatizirati in s tem pospešiti celoten postopek in se čim bolj približati izdelavi kart (predvsem rastrskih tematskih kart velikih meril) v realnem času (angl. real- time-processing). Določene faze v postopku obdelave posnetkov so se do nedavnega izvajale večinoma ročno in so zahtevale veliko časa, zato je bilo pomembno vse te faze med sabo čim bolj povezati in avtomatizirati. Predvsem je med členi v verigi avtomatizacije pomemben prav postopek geometričnih popravkov (ortorektifikacije) posnetkov. 12 Z razvojem digitalnih snemalnih sistemov, naprav za določanje položaja (globalni navigacijski satelitski sistemi, GNSS) in orientacije (inercialni navigacijski sistemi, INS), tehnologij prepoznavanja vzorcev, digitalnih baz podatkov in digitalnih modelov reliefa, lahko celoten postopek ortorektifikacije avtomatiziramo, kar zmanjša stroške in čas obdelave posnetkov (Slonecker in dr. 2009). V okviru raziskave smo razvili izviren postopek za samodejno ortorektifikacijo visokoločljivih optičnih satelitskih posnetkov, ki je del procesne verige STORM za samodejno obdelavo posnetkov v skoraj realnem času (Oštir in dr. 2014). Postopek (slika 1) brez posegov operaterja vhoden (t.i. kalibriran surov) satelitski posnetek, ki ga pripravi sprejemna postaja iz niza podatkov, prejetih s satelita, geometrično popravi tako, da dobimo ortofoto primerne kakovosti za zajem podatkov oz. za analize v geografskih informacijskih sistemih (GIS). Postopek je sestavljen iz več korakov: branje in priprava metapodatkov, samodejno določanje oslonilnih točk, geometrično modeliranje (izračun parametrov geometričnega modela) in ortorektifikacija. V sklopu raziskave je bil poudarek na korakih geometričnega modeliranja zajema posnetka in ortorektifikacije. izločanje grobih napak DIMAP PHR_SENSOR PRODUCT en RASTER_SENSOR DS_PH_0701_01272 image/jpeg image/jpeg Copyright@2013 CNES branje metapodatkov vhodni posnetek določanje oslonilnih točk geometrično modeliranje ortorektifikacija Slika 1: Samodejni postopek ortorektifikacije. Za izdelavo satelitskega ortofota je potrebnih več korakov. Čeprav je bilo v zadnjih letih več poskusov direktnega načina georeferenciranja posnetkov, ki za umestitev posnetka v prostor uporablja le podatke iz naprav za določanje položaja in orientacije, ta še vedno ne daje 13 dovolj natančnih rezultatov (Poli 2005). Zato smo v našem razvitem postopku ortorektifikacije uporabili še modul za samodejno določanje oslonilnih točk, ki je bil sicer izdelan za uporabo v postopku, ni pa bil v celoti predmet raziskav v knjigi. Oslonilne točke so pomemben vhodni podatek pri uporabi fizičnega oz. rigoroznega geometričnega modela, ki smo ga prilagodili za različne tipe visokoločljivih optičnih snemalnih sistemov. Geometrični model smo prilagodili tako za delo z vrstičnimi senzorji, ki so uporabljeni v večini večjih komercialnih satelitov, kot tudi polnoslikovnimi senzorji, ki so večinoma uporabljeni pri majhnih satelitih. Rezultat geometričnega modeliranja posnetka je določitev parametrov njegove notranje in zunanje orientacije. Ti se, skupaj z digitalnim modelom reliefa (DMR), uporabijo za samodejno ortorektifikacijo posnetkov. Osnovni namen raziskave je izdelati samodejen postopek ortorektifikacije, ki bo surovi posnetek ortorektificiral brez posredovanja operaterja najkasneje v dveh urah po prejetju, s čimer bo dostopen za uporabo v različnih aplikacijah, ki zahtevajo hitro odzivnost. Pri tem mora postopek zagotoviti položajno točnost izdelkov (ortofotov) v rangu dveh pikslov glede na referenčne podatke, s čimer bodo uporabni tudi za namene kartiranja v velikih merilih in združevanja podatkov (angl. data fusion). V raziskavi smo delovanje razvitega postopka prikazali s podatki satelitov serije RapidEye in WorldView-2, ki snemajo v vrstičnem načinu. Načrtovana je bila uporaba še več vrst senzorjev, tudi polnoslikovnih, vendar jih nismo uspeli pravočasno pridobiti. V raziskavi smo želeli pri ortorektifikaciji preizkusiti tudi nov državni lidarski DMR, ki pa je bil v času izvajanja raziskave še v obdelavi in za izbrana testna območja nedostopen. Kljub temu sta bila preizkušena senzorja, ki snemata v različnih prostorskih ločljivostih, primerna za vrednotenje delovanja vseh korakov razvitega postopka ortorektifikacije, kar je bil tudi naš glavni cilj. 14 1.1 Pregled obstoječih raziskav na področju samodejne ortorektifikacije V daljinskem zaznavanju in fotogrametriji se raziskovalci za potrebe ortorektifikacije posnetkov že vrsto let ukvarjajo z razvojem metod, ki bi bile čim bolj avtomatizirane oziroma samodejne. V zadnjih letih se je z uveljavitvijo zmogljivejših instrumentov, procesnih postaj in orodij povečalo število raziskav na tem področju. Čeprav narašča tudi število novih metod, je veliko poudarka tudi na popolnejši avtomatizaciji že preverjenih tehnik. V nadaljevanju poglavja so predstavljene dosedanje raziskave na področju samodejne ortorektifikacije optičnih satelitskih posnetkov. Postopek lahko v grobem delimo na tri dele, ki bodo tudi osnova za kratek pregled raziskav: geometrični model, slikovno korelacijo in samodejno georeferenciranje. 1.1.1 Geometrični model V zadnjih desetletjih je bilo izdelanih že veliko geometričnih modelov, ki so se razlikovali po zapletenosti in doseganju natančnosti. Splošno lahko vse geometrične modele delimo na dve vrsti (Toutin 2004): - fizični oz. rigorozni modeli in - empirični modeli. Rigorozni modeli poskušajo opisati fizične lastnosti zajema posnetkov. Za modeliranje zajema uporabljajo modificirane kolinearne enačbe, ki vključujejo še dodatne parametre. Empirični modeli so enostavnejši in teoretično tudi manj natančni. Odnos med slikovnimi in objektnimi koordinatami je tu vzpostavljen s polinomi, oziroma njihovimi kvocienti. Ob pomanjkanju podatkov o kalibraciji senzorja ali tirnici satelita lahko včasih empirični modeli nadomestijo rigorozni model. Enega prvih rigoroznih modelov sta predlagala Gugan in Dowman (1988) iz Londonske univerze UCL. Predlagata uvedbo dinamičnega modela orbitalnih parametrov (angl. dynamic orbital parameter model). V tem primeru se kolinearne enačbe razširi z dvema orbitalnima parametroma (prava anomalija in rektascenzija dvižnega vozla), s katerima se modelira premikanje satelita in rotacijo Zemlje. Geometrični model za orientacijo 15 posnetkov SPOT je razvil in v programu SPOTCHECK+ implementiral Kratky (1989). Model izračuna položaj satelita iz podatkov o tirnici, medtem ko se dogajanje med letom modelira s polinomskimi funkcijami (linearnimi ali kvadratnimi). Naslednji model uporablja princip orientiranja posnetkov, ki so ga razvili na Tehniški univerzi v Münchenu in implementirali v programu CLIC (Ebner in dr. 1992). Model temelji na razširjenih kolinearnih enačbah, parametri zunanje orientacije pa se določijo z blokovno izravnavo stereoposnetkov. Model so uporabljali za posnetke vrstičnega skenerja MOMS-2P. Na Inštitutu za fotogrametrijo in geoinformacije IPI v Hannovru so razvili program BLASPO za rektifikacijo satelitskih posnetkov (G. Konecny in dr. 1987). Za delovanje potrebuje le informacijo o satelitovi tirnici ter podatke o zasukih senzorja. Zelo fleksibilen rigorozen model za vrstične skenerje je izdelal tudi Toutin (2003). Model pri izvrednotenju upošteva distorzije platforme, senzorja, oblike Zemlje in kartografske projekcije. Westin (1990) je razvil geometrični model za rektifikacijo posnetkov SPOT. Uporabljen model tirnice satelita je še dodatno poenostavil iz elipse v krožnico, ki sicer med gibanjem spreminja polmer. Nekateri splošni geometrični modeli uporabljajo za modeliranje trajektorije in obnašanje senzorja kar polinomske funkcije. Poli (2005) je razvila model, ki uporablja polinomske funkcije druge stopnje, kar naredi model zelo splošen, vendar za določitev parametrov potrebuje opazovanje položaja in orientacije platforme. Splošen model senzorja za visokoločljive satelitske posnetke, kjer se tirnica in orientacija satelita modelirata s kubičnimi zlepki, so razvili tudi Weser in dr. (2008). Model lahko direktno vključi metapodatke v splošni model, kar omogoča neposredno georeferenciranje. Nekatere opisane geometrične modele so v zadnjih letih dopolnili in priredili za potrebe novih visokoločljivih satelitov (npr. Michalis in Dowman (2008) sta priredila dinamični model orbitalnih parametrov). Poleg rigoroznih modelov so za rektifikacijo satelitskih posnetkov v uporabi tudi empirični modeli, ki so pridobili na veljavi z uvedbo visokoločljivih satelitskih posnetkov, ki imajo v splošnem kompleksnejši način snemanja. Najbolj uporabljen je model racionalnih funkcij (angl. rational function model, RFM) z 80 parametri. Med prvimi, ki sta predstavila izčrpno študijo o modelu RFM, sta bila Tao in Hu (2001). Opisala sta metodi (odvisno 16 od reliefa ter neodvisno od njega) za izračun polinomskih koeficientov z iterativno neposredno izravnavo po metodi najmanjših kvadratov. Model sta preizkusila na aeroposnetkih in posnetkih SPOT ter dosegla primerljivo natančnost kot z rigoroznim modelom, vendar le z dodatnimi oslonilnimi točkami. Poseben primer empiričnih modelov so afine transformacije. Okamoto in dr. (1998) so afino transformacijo predlagali za probleme, ki nastanejo zaradi zelo ozkega kota gledanja snemalnega sistema satelita SPOT. Model uporablja lokalne koordinatne sisteme in elipsoidne višine, zato je pri izračunu potrebno upoštevati tudi ukrivljenost Zemlje. El-Manadili in Novak (1996) sta za geometrično modeliranje posnetkov SPOT predlagala model direktne linearne transformacije (angl. direct linear transformation, DLT). Metoda DLT uporablja 11 parametrov ter ne potrebuje parametrov notranje orientacije in podatkov o efemeridah satelita, saj rezultat dobi samo z oslonilnimi točkami. 1.1.2 Slikovna korelacija Pri postopku ortorektifikacije je zelo pomembno pridobiti kakovostne oslonilne točke, ki zanesljivo povežejo surov posnetek (preko geometričnega modela) s fizičnim Zemljinim površjem. Za samodejno določitev oslonilnih točk največkrat uporabljamo postopke slikovne korelacije oz. slikovnega ujemanja, ki slikovno korelacijo še nadgradi. Za avtomatizacijo postopka je samodejno iskanje oslonilnih točk s slikovno korelacijo ključnega pomena, zato je zelo pomembna izbira primerne metode slikovne korelacije, ki upošteva tip deformacij surovega posnetka v primerjavi z referenčnim. Metode slikovne korelacije za registracijo posnetkov daljinskega zaznavanja, kart, medicinskih posnetkov in drugih računalniško ustvarjenih slik obširno opišejo Brown (1992) ter Zitová in Flusser (2003). Metode, uporabne predvsem v daljinskem zaznavanju, opišeta in ovrednotita Fonseca in Manjunath (1996). V zadnjih letih so nekatere metode izboljšali in prilagodili različnim senzorjem. Izboljšan algoritem za določanje točk na visokoločljivih satelitskih posnetkih (IKONOS) sta predstavila npr. Xiong in Zhang (2009). Algoritem uspešno odkriva najbolj izrazite značilke (izrazite, bistvene lastnosti objektov 17 na posnetku), vendar je omejen na posnetke z ozkim kotom gledanja in upošteva samo premike in rotacije. Več študij je bilo izvedenih tudi na področju samodejne registracije posnetkov, ki uporabljajo razne tehnike slikovne korelacije za različne senzorje. Na srednje in nizkoločljivih posnetkih so delali Le Moigne in dr. (2002), ki so uspešno registrirali posnetke Landsat in AVHRR s križno korelacijo in valjčnimi transformacijami. Netanyahu in dr. (2004) so značilke na posnetkih Landsat pridobili z valjčnimi metodami, za registracijo pa so uporabili Hausdorffovo razdaljo. Algoritem za samodejno registracijo visokoločljivih posnetkov sta med drugimi razvila Zhang in Fraser (2007). V tem primeru sta značilke pridobila s Förstnerjevim operatorjem, korelacijo pa izvedla s hierarhičnimi piramidami. Posnetke IKONOS in QuickBird sta registrirala z afino transformacijo. Multimodalno metodo za samodejno registracijo optičnih (SPOT 4) in radarskih (ERS-2) posnetkov sta predstavila Inglada in Giros (2004). Za določanje homolognih točk sta uporabila razširjeno križno korelacijo. 1.1.3 Samodejno georeferenciranje Študij samodejnega georeferenciranja in ortorektifikacije različnih senzorjev ni zelo veliko. Gianinetto in Scaioni (2008) sta izdelala algoritem AGE za samodejno določitev oslonilnih točk. Adaptivno slikovno ujemanje najmanjših kvadratov (angl. adaptive least squares matching) se izvede med surovim posnetkom in geokodiranim posnetkom, iz katerega se dobi referenčne položajne koordinate točk. Višinske koordinate se samodejno dobi iz modela reliefa. Algoritem je bil, kljub malo slabšim doseženim točnostim ortofotov, uspešno preizkušen na posnetkih QuickBird, SPOT5 in IKONOS. Uspešen postopek samodejnega ortorektificiranja so razvili za radarske satelitske posnetke (Sheng in Alsdorf 2005). Najprej sta avtorja s filtriranjem dobila značilke, katerim sta določila homologne točke z normaliziranim korelacijskim koeficientom. V naslednjem koraku je sledila rektifikacija posnetka. Ta je potekala po posameznih kosih (angl. piecewise) posnetka, saj je vsak kos uporabljal svojo transformacijo. Rektifikacija po kosih je zelo uspešna pri posnetkih z zelo kompleksnimi deformacijami. Končne 18 ortorektificirane posnetke sta dobila z modelom reliefa SRTM (angl. Shuttle Radar Topography Mission). Skoraj sočasno ortorektifikacijo posnetkov iz letalske videokamere daljinsko vodenega letala je razvil Zhou (2009). Metoda deluje samodejno, vendar pri tem uporablja z GPS-om izmerjene oslonilne točke na območju, ki je določeno za snemanje. Značilke se izračunajo samo na prvem posnetku, na vseh naslednjih pa se jim skuša slediti (so vezne točke). Orientacijski parametri se izračunajo z blokovno izravnavo. Ortorektificirani posnetki se na koncu še samodejno mozaičijo. Na srednje in nizkoločljivih posnetkih sta delala tudi Eugenio in Marques (2003). Njuna metoda registracije temelji na določitvi robov, ki služijo kot značilke, ter nadaljnji afini transformaciji. Skoraj povsem samodejen postopek za ortorektifikacijo optičnih satelitskih posnetkov so izdelali Leprince in dr. (2007) za uporabo v raziskavah premikov (deformacij) zemeljskega površja. Za korelacijo med posnetki so izbrali metodo fazne korelacije, ki se zanaša na Fourierjev teorem premikov. To pomeni, da se relativne premike površja med podobnimi posnetki lahko dobi iz faznih razlik njihovih Fourierjevih transformacij. Pred ortorektifikacijo je treba ročno določiti oslonilne točke na vseh posnetkih. Celoten postopek je bil prilagojen le za obdelavo posnetkov SPOT 5. Model je bil kasneje implementiran v modul, ki deluje v okolju programa ENVI. Uspešen postopek za samodejno ortorektifikacijo visokoločljivih posnetkov Formosat-2 sta razvila Liu in Chen (2009). Njuna metoda deluje na korelaciji ploskovnih elementov iz letalskega ortofota in satelitskega posnetka. Ortorektifikacija poteka z rigoroznim transformacijskim modelom, dosežena natančnost pa je pod 1,5 piksla. Samodejno geometrično korekcijo sta razvila tudi Devaraj in Shah (2014). Ortorektifikacija srednje ločljivih posnetkov CBERS-2 in CBERS-2B poteka z referenčnimi posnetki Landsat. Pred korelacijo se ortorektificirane posnetke Landsat transformira v surovo obliko (posnetek ima znova geometrična popačenja). Nato se registracija izvede z afino transformacijo. 1.2 Struktura knjige Knjiga je razdeljena na sedem poglavij. V prvem so predstavljeni motivacija za raziskavo in pretekle raziskave na področju samodejne ortorektifikacije. 19 Drugo poglavje je namenjeno razlagi izhodišč in temeljnih pojmov. Opisane so vrste senzorjev in ravni obdelave satelitskih posnetkov. V grobem sta predstavljena ortorektifikacija in postopek razvitega ortorektifikacijskega postopka. V tretjem poglavju sledi opis uporabljenih podatkov ter satelitov, s katerimi so bili podatki zajeti. Četrto poglavje predstavlja algoritem za samodejno določanje oslonilnih točk, implementirana modula verige za branje in pripravo metapodatkov ter določitve oslonilnih točk. Prikazano je tudi testiranje razvitega postopka z visokoločljivimi posnetki RapidEye in zelo visokoločljivimi posnetki WorldView-2. Poglavje vsebuje pregled rezultatov in analizo dosežene točnosti določevanja točk. Peto poglavje je namenjeno obširni predstavitvi geometričnih modelov vrstičnega in polnoslikovnega senzorja, ki sta uporabljena v postopku. Predstavljena je sestava modelov in njihovo delovanje od priprave koordinat oslonilnih točk do izravnave in končnih orientacijskih parametrov. Tudi v tem poglavju je opisan implementiran modul za izračun parametrov geometričnega modela in testiranje postopka z razpoložljivimi podatki ter pregled rezultatov in njihove točnosti, ki jih dosežemo z geometričnim modelom. Šesto poglavje je namenjeno obravnavi metod ortorektifikacije in predstavitvi indirektnega modela, ki je bil uporabljen v raziskavi. Poleg teoretičnega dela je predstavljen tudi implementiran modul za izvedbo ortorektifikacije ter izvedeni testi na podatkih. Analiza doseženih rezultatov o točnosti ortorektifikacije je prikazana na koncu poglavja. Zadnje, sedmo poglavje je namenjeno končnim ugotovitvam raziskave ter razpravi o možnostih za nadaljnje delo in nadgradnjo razvitega postopka. 20 2 IZHODIŠČA IN TEMELJNI POJMI V poglavju so predstavljena osnovna izhodišča raziskave in temeljni pojmi, ki so potrebni za razumevanje samega postopka. Najprej so predstavljene vrste senzorjev in geometrije zajema posnetkov ter ravni obdelave podatkov. Nato je razložen osnovni postopek ortorektifikacije in opisana samodejna procesna veriga za satelitske posnetke, katere del je tudi postopek, opisan v knjigi. 2.1 Geometrija zajema posnetkov glede na vrsto senzorja Optični satelitski posnetek nastane v snemalnem sistemu satelita, ki sprejema svetlobo (radiacije različnih valovnih dolžin), ki prihaja od Sonca in se odbije od opazovanega predmeta (površja Zemlje). Sam snemalni sistem je sestavljen iz več delov, med njimi pa je za zapis slike ključen slikovni detektor, ki zajame odbito svetlobo in jo pretvori v električni naboj. Magnituda električnega naboja je sorazmerna z intenziteto prejete svetlobe. Za satelitsko snemanje poznamo več različnih izvedb detektorjev: najbolj sta v uporabi tehnologiji CCD (Charge-Coupled Device) in CMOS (Complementary Metal- Oxide Semiconductor). Detektor CCD je sestavljen iz več slikovnih elementov, ki zbirajo, shranjujejo in prenašajo električni naboj iz enega elementa na drugega. Dobljeni naboji se prenesejo v element za odčitavanje in se jih pred vhodom v analogno-digitalni pretvornik še dodatno ojača. Vsak slikovni element bo predstavljal en piksel na posnetku. Tudi detektor CMOS je sestavljen iz več slikovnih elementov, vendar ima vsak element svoj ojačevalec in se ga lahko obravnava in bere ločeno od ostalih. Vsaka od tehnologij ima svoje prednosti in slabosti, vendar se je v preteklosti zaradi zgodnje iznajdbe, nižje cene in boljše izdelave slik večinoma uporabljalo detektorje CCD (Schenk 1999). Konkurenčni detektorji 21 CMOS (Bigas in dr. 2006) so se pojavili kasneje in se do pred kratkim v vesoljske namene niso uporabljali. Ker detektorji CMOS porabijo manj energije in so bolj naklonjeni miniaturizaciji, so trenutno zelo priljubljeni za vgradnjo v majhne satelite, kjer sta količina energije in prostor zelo omejena. Izbira tehnologije izdelave detektorja je odvisna tudi od oblike oziroma izvedbe detektorja, ki je lahko vzdolžen (vrstičen) ali polnoslikovni (matričen). Zaradi boljše preizkušenosti in zanesljivega delovanja se v večje satelite (vsi delujejo v vrstičnem načinu snemanja) za opazovanje Zemlje vgrajuje le detektorje CCD. Kako bo v prihodnje, pa bo odvisno od možnosti prilagoditve obeh tehnologij za prihajajoče potrebe daljinskega zaznavanja. Pri modeliranju geometrije zajema posnetka sta zelo pomembni dve lastnosti sistema za zajem posnetkov: velikost piksla in goriščna razdalja. Prva se navezuje na sam detektor oziroma na velikost slikovnega elementa (piksla) iz katerih je sestavljen detektor. Manjši kot je element, večjo prostorsko ločljivost ima posnetek pri enakih ostalih parametrih. Pri tem je potrebno paziti na količino sprejete svetlobe, saj so lahko ob neprimerni optiki in premajhnih elementih detektorja vrednosti na posnetku polne šuma in zato neuporabne. Z detektorjem, ki ima slikovne elemente velikosti nekaj μm in z optiko oz. teleskopom, ki ima primerno goriščno razdaljo in premer vstopne odprtine, lahko sateliti dosežejo prostorske ločljivosti tudi pod 0,5 metra. Druga lastnost je povezana s teleskopom sistema za zajem posnetkov. Od zgradbe teleskopa je odvisna velikost goriščne razdalje, ki tudi določa prostorsko ločljivost posnetka. Večja kot je goriščna razdalja, boljša je ločljivost. Teleskopi večjih satelitov imajo goriščno razdaljo velikosti več metrov, kar dosežejo predvsem z zrcali (t.i. reflektorji) in lečami (t.i. refraktorji). Večji sateliti, s težo nekje od 300 kg do 3000 kg, ponavadi uporabljajo teleskope tipa Cassegrain (npr. Landsat), Korsch (npr. Pleiades) ali TMA (Three-Mirror Anastigmat; npr. RapidEye), manjši, s težo pod 300 kg, pa večinoma refraktorje. Geometrija snemanja je odvisna predvsem od načina zajemanja posnetka. Zaradi potrebe po visoki prostorski ločljivosti in večji velikosti snemanega območja večina sodobnih satelitov uporablja tehniko vzdolžnega snemanja, kjer senzor zajema vrstico za vrstico v določenih časovnih intervalih. Skupek vseh zajetih vrstic nato tvori posnetek. V tem primeru je lahko slikovna 22 ravnina sestavljena iz več detektorjev, ki tvorijo vrstico slikovnih elementov. Geometrija je zelo kompleksna, saj ima vsaka vrstica med snemanjem svoj unikaten čas in perspektivni center. Drugi, manj pogosti način snemanja, je polnoslikovni. V tem primeru celoten posnetek nastane v centralni projekciji z enim perspektivnim centrom, saj ima detektor obliko matrike. Ker celoten posnetek nastane v zelo kratkem času, je geometrijo zajema lažje modelirati. Zaradi večje deformiranosti posnetka (kot gledanja) in občutljivosti na šum se ta način ni uveljavil za snemanje iz vesolja. Več o obeh načinih zajemanja posnetka je zapisanega v poglavju 5. 2.2 Ravni obdelave posnetkov Satelit posneto sliko pošlje preko oddajnika v sprejemno postajo, kjer se posnetek dodatno obdela. Komercialna podjetja, ki imajo v lasti satelite ali razpolagajo s posnetki, obdelajo posnetke do različnih stopenj (angl. processing levels) glede na potrebe uporabnikov. Vsak ponudnik posnetkov uporablja svojo notacijo in definicije obdelave (npr. Gutman in Ignatov 1997), vendar so si definicije različnih proizvajalcev podobne in se ponavadi razlikujejo le v detajlih. Če povzamemo skupne značilnosti obstoječih definicij, lahko ravni razdelimo v naslednje razrede: - Raven 0: surov, neobdelan posnetek. - Raven 1A: upoštevani radiometrični popravki zaradi variacije elementov znotraj detektorja. Popravke se računa glede na kalibracijo elementov v laboratoriju in po izstrelitvi. - Raven 1B: radiometrično (kot raven 1A) in geometrično popravljen posnetek. Popravljene so sistematične napake, ki nastanejo zaradi hitrega premikanja satelita in rotacije Zemlje. Projekcija na referenčni elipsoid ali horizontalno ploskev. - Raven 2A: prevzorčen posnetek v standardno kartografsko projekcijo glede na položaj satelita v času zajema. - Raven 2B: georeferenciran oz. ortorektificiran posnetek. Za rektifikacijo so uporabljene oslonilne točke in DMR. - Raven 3: mozaik posnetkov in specifični izdelki, za katere je bil določen instrument namenjen. 23 Nekateri ponudniki satelitskih posnetkov označijo ortorektificiran posnetek kot raven 3A. Tudi cene posnetkov so odvisne od ravni obdelave. Na primer, arhivski pankromatski posnetek WorldView-2 ločljivosti 0,5 m pri podjetju e- GEOS (e-GEOS 2017) stane 14,5 $/km2 za ravni 1B ali 2A in 27,5 $/km2 za ortorektificiran posnetek (raven 2B). V raziskavi smo kot vhodne posnetke v postopku uporabljali le posnetke ravni 1A oziroma le radiometrično popravljene posnetke. Le geometrično nepopravljeni posnetki se lahko uporabijo v izdelanem rigoroznem geometričnem modelu, ki opiše geometrijo zajema posnetka. Za ortorektifikacijo naslednjih stopenj so potrebni drugi modeli, kot je recimo RFM, ki uporabljajo racionalne polinomske koeficiente (angl. rational polynomial coefficients, RPC). 2.3 Ortorektifikacija V knjigi se bomo ukvarjali z ortorektifikacijo optičnih satelitskih posnetkov. Izdelava ortorektificiranih posnetkov je nujna, če želimo posnetke uporabljati za prostorske analize z uporabo različnih podatkov (npr. v geografskih informacijskih sistemih). Potreba po ortorektifikaciji satelitskih posnetkov je nastala takoj po zmožnosti snemanja satelitov v višji prostorski ločljivosti v sedemdesetih letih. Ortorektifikacija je pretvorba posnetka v ortogonalno projekcijo na izbrano ravnino (običajno je to ravnina izbrane kartografske projekcije), kar dosežemo s postopkom zmanjšanja različnih geometričnih popačenj, ki jih vsebuje izvorni (surovi) posnetek. Najpogostejši dejavniki za nastanek geometričnih popačenj so: - orientacija senzorja, - napake senzorja, - kot gledanja senzorja, - relief posnetega območja in - ukrivljenost Zemlje. S postopkom ortorektifikacije se iz surovega posnetka izdela ortorektificiran posnetek oz. ortofoto. Pri tem se uporabijo parametri notranje in zunanje 24 orientacije, ki jih dobimo z izravnavo, in digitalni model reliefa za obravnavano območje. Po ortorektifikaciji je vsaka meritev, ki se jo opravi na ortorektificiranem posnetku, enaka meritvi, ki bi se jo opravilo na zemeljski površini (ob upoštevanju lokalnega merila, ki je odvisen od lokalnega popačenja kartografske projekcije) (Marsetič in Oštir 2007). Pri ortorektifikaciji se posnetek z geometrično transformacijo (georeferenciranjem ali ortorektifikacijo) umesti v izbran koordinatni sistem. Radiometrične vrednosti v ortofotu dobimo s prevzorčenjem vrednosti pikslov surovega posnetka. Najpogostejše tehnike prevzorčenja so metoda najbližjega soseda, bilinearna interpolacija in kubična konvolucija. Ponavadi se za izhodni posnetek uporabi enaka velikost pikslov glede na vhodni posnetek (tolerira se 10 % –20 % razlika). S posebnimi tehnikami pa se lahko doseže tudi manjšo velikost celic brez znatnih izgub informacij in napak (Planet Labs 2017). 2.4 Samodejna procesna veriga Postopek samodejne ortorektifikacije optičnih satelitskih posnetkov je sestavni del samodejne procesne verige (Oštir in dr. 2014), ki je bila razvita v okviru Centra odličnosti Vesolje, znanost in tehnologije (Vesolje-SI). Veriga, imenovana STORM, obsega vse korake obdelave posnetkov od t.i. kalibriranega surovega posnetka (angl. sensor-corrected image), ki ga ponavadi označujemo z ravnijo 1A, do obdelanih posnetkov in izdelkov, dostopnih preko spleta. Veriga obsega: - modul za samodejno ortorektifiacijo, ki je predmet knjige, - modul za samodejne radiometrične popravke, ki je sestavljen iz atmosferskih (uporablja se programsko opremo ATCOR švicarskega podjetja ReSe Applications Schläpfer) in topografskih popravkov (delujejo na osnovi modela anizotropne osvetlitve), ter - module za samodejni izračun enostavnih izdelkov (vegetacijskega indeksa NDVI, njegovih sprememb in karto sprememb). Iz vhodnega posnetka torej nastane več končnih tematskih izdelkov, ki se skupaj s parametri obdelave shranijo v podatkovno bazo. Na koncu procesne 25 verige več povezanih servisov dostavlja podatke do končnih uporabnikov v obliki rastrske spletne karte. Razvoj samodejne procesne verige se je začel leta 2011 s posnetki satelitov RapidEye, trenutno pa lahko veriga obdeluje več različnih posnetkov (RapidEye, WorldView-2, Pleiades in SPOT 6). Veriga torej podpira tako visokoločljive posnetke (prostorska ločljivost med 2 m in 20 m) kot zelo visokoločljive posnetke (prostorska ločljivost 2 m in manj). Njeno delovanje je sicer omejeno na območje Slovenije z okolico, kjer imamo na voljo vse zahtevane podatke (DMR, digitalne podatke o cestah in državni ortofoto) s primerno točnostjo. Kljub temu, da je bila prvotno veriga namenjena uporabi v Sloveniji, bodo nadgradnje podpirale tudi podatke drugih držav in različne koordinatne sisteme. Samodejni postopek ortorektifikacije je zelo pomemben člen procesne verige STORM. Namenjen je odpravljanju geometričnih popačenj in skrbi za točno pozicioniranje posnetka v izbran koordinatni sistem. Postopek je razdeljen v več korakov (slika 1). Nekateri so medsebojno zelo povezani in delujejo v iteracijah (geometrični model in izločanje grobih napak), drugi pa so samostojni in lahko tvorijo tudi samostojen modul (npr. določanje točk). Cel postopek je popolnoma samodejen in se začne izvajati, ko procesni strežnik prejme vhodni posnetek. Najprej postopek prebere metapodatke posnetka in izloči ter shrani tiste, ki se bodo uporabili v nadaljnjem poteku postopka. Sledi samodejno določanje (detekcija in merjenje koordinat) oslonilnih točk, geometrično modeliranje ter na koncu ortorektifikacija. Vsi omenjeni postopki so implementirani v štiri osnovne module (slika 2): - modul za branje in pripravo metapodatkov (glej podpoglavje 0), - modul za samodejno določanje oslonilnih točk (glej podpoglavje 4.2.2), - modul za izračun parametrov geometričnega modela (glej podpoglavje 0) in - modul za izvedbo ortorektifikacije (glej podpoglavje 6.3.1). Vsi moduli so bili napisani (kodirani) v programskem jeziku IDL (Interactive Data Language). Izjema je modul za iskanje oslonilnih točk, ki je delno napisan v programskem jeziku C++, delno pa v IDL. IDL omogoča pisanje 26 hitrih programov, ki lahko delujejo z velikimi matrikami, kar je v drugih razvojnih okoljih pogosto težava. Izdelani programi uspešno obdelajo datoteke velikosti več gigabajtov. Pozorni moramo biti le na razpoložljivi prostor na bralno-pisalnem pomnilniku, kjer poteka shranjevanje vseh začetnih, vmesnih in končnih rezultatov. Slika 2: Samodejni postopek za ortorektifikacijo je sestavljen iz štirih osnovnih modulov. Razviti algoritmi potrebujejo za delovanje več pomožnih podatkov, ki jih mora zagotavljati namenski strežnik z zanesljivo in hitro povezavo do procesnega računalnika. Teh podatkov (DMR, državni ortofoto, podatki o cestah ipd.) je samo za območje Slovenije približno 450 GB. Z načrtovano 27 uporabo podatkov boljših ločljivosti in podatkov drugih držav se bo v prihodnje količina precej povečala. Za potrebe Slovenije je trenutna rešitev zadovoljiva, v primeru širitve pa bo potrebno celostno načrtovanje shranjevanja in rokovanja tako s pomožnimi podatki, kot s posnetki in rezultati. Rezultati postopka so ortorektificiran posnetek ter podatki o obdelavi, vrednosti parametrov zunanje in notranje orientacije in dosežene točnosti teh parametrov, ki se zapišejo v tekstovno datoteko. 28 3 UPORABLJENI PODATKI Delovanje izdelanega ortorektifikacijskega postopka smo ocenili s poskusi na podatkih, ki smo jih imeli na voljo in so bili primerni za ta namen. Vse poskuse smo opravili na treh visokoločljivih posnetkih satelitov RapidEye in treh zelo visokoločljivih posnetkih satelita WorldView-2. Sateliti kot tudi njihovi posnetki se razlikujejo v več lastnostih, za ortorektifikacijo pa sta pomembna predvsem dva – prostorska ločljivost in širina pasu snemanja. Preizkušeni posnetki omenjenih satelitov so reprezentativni predstavniki v svojih ločljivostnih razredih za večino posnetkov, ki se trenutno uporabljajo v raznih aplikacijah (npr. kartiranje, analize v GIS). Oba uporabljena satelita snemata v vrstičnem načinu. Trenutno v polnoslikovnem načinu snemajo le nekateri majhni sateliti, njihovi posnetki pa so v veliko primerih neuporabni za ortorektifikacijo (pomanjkanje metapodatkov, slaba kakovost ipd.) ali jih je težko pridobiti. Posnetki slovenskega majhnega satelita, ki deluje v polnoslikovnem načinu in ki bo v prihodnje uporabljal model, razvit v sklopu raziskave, še niso na voljo. Zaradi tega delovanje in rezultati obdelave polnoslikovnega geometričnega modela v knjigi ne bodo prikazani, smo pa v knjigi opisali teorijo, ki je potrebna za njegovo implementacijo. 3.1 RapidEye Sateliti RapidEye so del komercialnega sistema za daljinsko zaznavanje Zemlje. Konstelacija je sestavljena iz petih minisatelitov, ki lahko v enem dnevu posnamejo do 5 milijonov km2. Sateliti krožijo v sončno sinhroni tirnici na višini okoli 630 km in snemajo z vrstičnimi senzorji v petih spektralnih kanalih. Prostorska ločljivost posnetkov na tleh je okoli 6,5 m, kar RapidEye 29 uvršča med visokoločljive satelite. Sistem satelitov, kontrolni segment in zemeljske postaje je izdelalo podjetje MDA (MacDonald, Dettwiler and Associates Ltd) s podizvajalci. Satelite so izstrelili 29. avgusta 2008 iz izstrelišča Baikonur v Kazahstanu. Sistem je trenutno last podjetja Planet Labs (Planet Labs 2017), ki z njim tudi upravlja. Glavne lastnosti satelitov so prikazane v preglednici 1. Razvoj samodejnega ortorektifikacijskega postopka se je začel z 19 posnetki RapidEye, ki so bili posneti nad Slovenijo med letoma 2011 in 2013. Postopek trenutno dosega dobre rezultate z vsemi posnetki, ki pokrivajo vse mesece v letu in imajo zaradi tega tudi zelo različno osvetljenost in stanje rastja. V knjigi so predstavljeni rezultati za tri izbrane posnetke, ki so bili posneti na območjih z različno pokrovnostjo in razgibanostjo terena ter v različnih letih, žal vsi v poletnem času. Uporabljeni posnetki predstavljajo dober nabor za preizkušanje zmogljivosti in robustnosti samodejnega postopka. Preglednica 1: Lastnosti satelitov RapidEye (eoPortal 2017a). Tirnica sončno sinhrona, višina leta ~630 km Prostorska ločljivost nadir ~6,5 m Širina pasu snemanja 77 km Detektor vrstični CCD, 12000 pikslov (elementov) Optika TMA, goriščna razdalja 633 mm Čas ponovnega obiska 1 dan z zasuki (vsi sateliti) Radiometrična ločljivost 12 bitov Masa satelita ob izstrelitvi 156 kg modra: 440-510 nm zelena: 520-590 nm Kanali rdeča: 630-685 nm rob rdeče: 690-730 nm bližnja IR: 760-850 nm Prostor za shranjevanje podatkov 48 Gbit Hitrost prenosa podatkov 80 Mbit/s Preglednica 2: Lastnosti uporabljenih posnetkov RapidEye. Posnetek Datum Kot gledanja [°] Velikost [piksel] Oblačnost [%] Radgona 11. 6. 2011 3,65 11.802  7604 0 Koper 18. 8. 2012 6,74 11.796  10.753 1 Bohinj 1. 8. 2013 0,12 11.802  11.223 3 30 Slika 3: Uporabljeni surovi posnetki RapidEye (spodaj) in njihov položaj (zgoraj). 31 Nekatere lastnosti testnih posnetkov prikazuje preglednica 2. Posnetke, ki smo jih poimenovali po krajih oz. območjih, ki se pojavijo na posnetkih, vsebujejo zelo malo oblakov. Posnetek Bohinj je bil posnet skoraj v nadirju, medtem ko imata preostala dva le manjši kot gledanja. Zaradi bližine državne meje je posnetek Radgona manjši, s površino 3800 km2, Bohinj in Koper pa obsegata 5500 km2. Posnetki prikazujejo tri morfološko različne regije v Sloveniji. Posnetek Koper obsega večinoma gričevnato obalno območje, Radgona ravninski teren z območji gričevja, posnetek Bohinj pa je najbolj razgiban z visokimi gorami in ozkimi dolinami. Vsi posnetki poleg ozemlja Slovenije vsebujejo tudi dele sosednjih držav (slika 3). 3.2 WorldView-2 Satelit WorldView-2 so izstrelili 8. oktobra 2009 iz izstrelišča Vanderberg v Kaliforniji, ZDA. Izdelalo ga je podjetje Ball Aerospace and Technologies Corporation, lastnik in upravljalec pa je podjetje DigitalGlobe Inc. (DigitalGlobe 2017a). Snema v zelo visoki ločljivosti in do izstrelitve njegovega naslednika (WorldView-3 leta 2014) je veljal za komercialni satelit z najboljšo prostorsko ločljivostjo (0,46 m). WorldView-2 kroži v sončno sinhroni tirnici na višini okrog 770 km. Snemalni sistem lahko pridobi posnetek v pankromatskem in/ali večspektralnem načinu, pri čemer ima skupno devet kanalov. V primerjavi s sateliti RapidEye je WorldView-2 veliko večji, ima več načinov snemanja in tudi boljšo gibljivost. Ena večjih razlik med satelitoma je prostorska ločljivost. Na zelo visokoločljivih posnetkih je prisotnih veliko več detajlov, kar vpliva tudi na vrsto uporabljenih metod za obdelavo in način obdelave posnetkov. Glavne lastnosti satelita so prikazane v preglednici 3. Za raziskave smo imeli na razpolago tri posnetke WorldView-2, ki so bili pridobljeni v okviru Centra odličnosti Vesolje-SI. Prvi je bil posnet poleti leta 2010, druga dva pa v letu 2012 in sicer v istem preletu. Zaradi tega imata posnetka tudi podobne značilnosti (preglednica 4). Vsi uporabljeni posnetki so večspektralni s prostorsko ločljivostjo okoli 2 m. Manjša ločljivost od nominalne je posledica kota gledanja, ki v vseh primerih presega deset stopinj. Največji kot gledanja ima posnetek Ljubljana, satelit pa je bil med 32 snemanjem usmerjen pravokotno na tirnico. Pri ostalih posnetkih je satelit snemal pod kotom v smeri leta (vzdolž tirnice). Vsi posnetki (slika 4) imajo približno enako velikost in sicer okrog 255 km2, kar je precej manj od posnetkov RapidEye. Pokritost z oblaki je majhna in nima vpliva na obdelavo. Preglednica 3: Lastnosti satelita WorldView-2 (eoPortal 2017b). Tirnica sončno sinhrona, višina leta ~770 km 0,46 m PAN Prostorska ločljivost nadir 1,84 m MS Širina pasu snemanja 16,4 km Detektor vrstični CCD, 35000 pikslov PAN, 9300 pikslov MS Optika TMA, goriščna razdalja 13,3 m Čas ponovnega obiska 1,1 dan z zasuki Radiometrična ločljivost 11 bitov Masa satelita ob izstrelitvi 2800 kg pan: 450-800 nm svetlo modra: 400-450 nm modra: 450-510 nm zelena: 510-580 nm Kanali rumena: 585-625 nm rdeča: 630-690 nm rob rdeče: 705-745 nm bližnja IR: 770-895 nm bližnja IR: 860-1040 nm Prostor za shranjevanje podatkov 2199 Gbit Hitrost prenosa podatkov 800 Mbit/s Posnetki se razlikujejo tudi glede pokrovnosti tal. Posnetek Ljubljana prikazuje del mesta Ljubljane, ravninsko in poljedelsko Ljubljansko barje ter gričevnat svet. Posnetka Jesenice in Bled obsegata gorata območja z dolinami in planotami. Gorat je predvsem posnetek Jesenice, katerega severni del zaobjema območje Karavank. Preglednica 4: Lastnosti uporabljenih posnetkov WorldView-2. Posnetek Datum Kot gledanja [°] Velikost [piksel] Oblačnost [%] Ljubljana 1. 8. 2010 18,8 7603  8766 0 Jesenice 11. 9. 2012 12,8 7415  8472 3,5 Bled 11. 9. 2012 12,5 7409  8464 5 33 Slika 4: Uporabljeni surovi posnetki WorldView-2 (spodaj) in njihov položaj (zgoraj). 34 4 SAMODEJNO DOLOČANJE OSLONILNIH TOČK Oslonilne točke (OT) so v večini metod ortorektifikacije nujne za doseganje visoke položajne točnosti ortofotov. V tem poglavju je predstavljena metoda za samodejno določevanje oslonilnih točk, ki sicer večinoma ni bila razvita v okviru raziskave, je pa nujna za samodejen potek celega postopka. Zato je v knjigi njen potek in osnovno delovanje predstavljeno bolj na kratko. Vsi sodobni sateliti za opazovanje Zemlje so opremljeni s sistemi, ki jim določajo položaj (npr. globalni navigacijski satelitski sistemi (GNSS)) in orientacijo (npr. sledilec zvezd oz. angl. star tracker) satelita. Omenjeni sistemi so že zelo zmogljivi in se lahko uporabijo za prostorsko pozicioniranje posnetkov. Čeprav je ta direktni način zelo preprost in hiter, pa ne daje visoko natančnih rezultatov (Poli 2005), ki jih v določenih aplikacijah potrebujemo. Na primer, položajna točnost (CE901) posnetkov RapidEye ravni 1B je okoli 50 m (Planet Labs 2017), posnetkov WorldView-2 ravni obdelave Basic (enakovredno ravni 1A) pa 5 m (DigitalGlobe 2017b). Pri tem je treba poudariti, da so omenjeni posnetki že delno obdelani (uporabljene so tudi oslonilne točke, ki pa niso določene z visoko natančnostjo) in da navedene natančnosti veljajo le za majhne kote gledanja ter ne upoštevajo vpliva reliefa. Ker podatki o položaju in orientaciji satelita ne dosegajo želene natančnosti, ki je uveljavljena za tradicionalne geodetske topografske izdelke, jih je treba popraviti z natančnimi referenčnimi podatki. Za najbolj primerne 1 Circular Error 90 (CE90) se uporablja za opis položajne napake, povezane z 90 % vsebine posnetka glede na uporabljene referenčne točke. Če je CE90 na primer 50 m, potem ima 90 % vsebine posnetka napako manjšo od 50 m. 35 so se izkazale trirazsežne (3R) točke v izbranem koordinatnem sistemu, ki ležijo na Zemljinem površju (točke z ravninskimi koordinatami x in y ter višinsko koordinato z). Takim točkam pravimo oslonilne točke in so v fotogrametriji v uporabi že vrsto let, tudi za izdelavo državnega ortofota. Njihova naloga je, da povežejo vhodni surovi posnetek z referenčnimi podatki, ki so lahko v obliki rastrskih (npr. ortofoto), vektorskih (npr. podatki o cestah) ali tekstovnih podatkov (npr. niz točk z GNSS-om izmerjenih koordinat). Oslonilne točke se ponavadi dobijo ročno iz referenčnih virov ali s terenskimi meritvami, kar je lahko zelo zamudno. V zadnjem obdobju so se razvile tehnike, ki točke pridobivajo polsamodejno ali celo povsem samodejno. Zaradi spremenljivosti tako vhodnih posnetkov kot referenčnih podatkov, tudi razvite samodejne tehnike niso bile povsem zanesljive za določanje oslonilnih točk. Boljšo točnost in zanesljivost postopka pa smo dosegli z novim postopkom, ki je opisan v naslednjem poglavju. Uporabnost določenih oslonilnih točk je sicer odvisna tudi od nadaljnjega postopka obdelave oz. geometričnega modela. 4.1 Postopek samodejnega določanja oslonilnih točk V raziskavi uporabljen postopek za samodejno določanje točk je razdeljen v dve fazi: določanje z vektorskimi podatki o cestah in določanje z državnim ortofotom (GURS 2017a) s slikovnim ujemanjem po metodi najmanjših kvadratov. Visokoločljivi posnetki potrebujejo le prvo fazo, posnetki zelo visoke prostorske ločljivosti pa obe. Omenjeni postopek se je v obsežnih preizkušanjih izkazal kot najbolj robusten in zanesljiv. Prvotna ideja o uporabi državnega ortofota za referenco in korelacijo s satelitskimi posnetki z algoritmoma SIFT (Scale-Invariant Feature Transform; Lowe 1999) ali SURF (Speeded Up Robust Features; Bay in dr. 2008) se je zaradi prevelike časovne variabilnosti med viri in različnih prostorskih ločljivosti izkazala za neuspešno. V okviru Centra odličnosti Vesolje-SI smo za potrebe procesne verige razvili nov postopek samodejnega določanja oslonilnih točk (Zaletelj in dr. 2013), ki ga prikazuje slika 5. Postopek je plod sodelovanja med znanstveniki s Fakultete za elektrotehniko, kjer je bil algoritem tudi narejen, in znanstveniki (večinoma iz ZRC SAZU), ki smo delali na drugih delih procesne verige. Veliko 36 rešitev je bilo posledica zelo tesne interakcije med sodelujočimi skupinami. Razvita metoda temelji na digitalnih podatkih o cestah, ki so stabilne strukture na zemeljskem površju, se redko spreminjajo, njihov položaj pa je dobro določen. Poleg tega so dobro vidne na satelitskih posnetkih in jih je možno zaznati z razmeroma preprostimi algoritmi za obdelavo slik. Metoda določanja oslonilnih točk uporablja za referenco vektorske podatke o cestah, do katerih smo lahko dostopali v okviru projekta. Za Slovenijo smo uporabili dva državna vira za digitalne podatke o cestah: podatke o cestah iz DTK5 (topografski podatki merila 1:5000; GURS 2017b) in podatke o cestah iz ZKGJI (Zbirni kataster gospodarske javne infrastrukture; GURS 2017c). Osnovni vir so bili podatki o cestah iz DTK5, ki so bili večinoma zajeti med letoma 2002 in 2008 iz stereoparov cikličnega aerosnemanja. Kljub temu, da podatki že nekaj let niso bili posodobljeni in pokrivajo le približno 60 % ozemlja Slovenije, imajo večjo gostoto in položajno točnost cest (ocenjeno na 1 m) kot drugi uporabljeni viri. Podatke smo rastrirali na podlagi atributa širine cest. Vektorske podatke o cestah iz ZKGJI smo uporabili le na območjih, kjer ni bilo podatka iz DTK5. ZKGJI je redno vzdrževan, podatki za avtoceste in pomembnejše ceste pa so pridobljeni z meritvami GNSS, ki imajo deklarirano točnost pod 1 m, večinoma pa so iz državnega ortofota in imajo deklarirano točnost od 1 m do 5 m (GURS 2017c). V ZKGJI ceste nimajo podatka o širini, zato smo jo ocenili na podlagi podatka o kategoriji cest in šele nato rastrirali. Za okoliške države smo uporabili prosto dostopne podatke organizacije OpenStreetMap (OSM; OpenStreetMap 2017). OSM izdelujejo in dopolnjujejo prostovoljci, zato ima neznano in spremenljivo točnost (razlikuje se glede na območje). Tudi v tem primeru smo podatek o širini ocenili iz podatka o kategoriji cest. Referenčne vektorske podatke o cestah smo torej rastrirali glede na njihovo širino in jih pripravili za uporabo v postopku. Pred začetkom iskanja korelacije med podatki o cestah in posnetkom je treba vhodne podatke dodatno obdelati. Najprej se iz posnetka odstrani oblake. Masko oblakov se lahko dobi iz histograma, v nekaterih primerih pa je dosegljiva tudi v metapodatkih. Nato se iz baze prebere le ceste, ki se prekrivajo s posnetkom in se jih razdeli na več manjših pravokotnih območij (izsekov), katerih število je odvisno od velikosti posnetka. Rastrirane podatke o cestah se za potrebe 37 korelacije transformira v rastrske podatke oddaljenosti od cest, kjer vsaka vrednost piksla pomeni razdaljo do najbližjega cestnega piksla. Tako pripravljeni cestni podatki so vhodni podatek za korelacijo. Na vhodnem posnetku se najprej z morfološkim operatorjem top-hat, ki zaznava svetle linije, izloči ceste. Čeprav je ta operator pri izločanju cest zelo uspešen, je med rezultati običajno tudi nekaj napačnih detekcij in manjkajočih cestnih odsekov, ki otežujejo natančno korelacijo z referenčnimi cestami. Kljub temu je količina napačnih detekcij oz. nezaznanih cest majhna in v večini primerov ne vpliva na nadaljnji postopek. Nato se detektirane ceste zasuka za kot, ki je izračunan iz metapodatkov posnetka in predstavlja azimut tirnice satelita. Po zasuku so referenčne in detektirane ceste približno poravnane. Pripravi cest sledi poravnava po izsekih v več korakih. V prvem koraku se izvede groba poravnava izsekov, s katero dobimo niz perspektivnih transformacijskih parametrov, ki povežejo koordinate pikslov na izsekih z zaznanimi cestami in izsekih s podatki oddaljenosti od cest. Devet parametrov transformacije izračunamo z uporabo parov točk. Enakomerno razporejene točke na podatkih zaznanih cest imajo koordinate t  [ ( t ) x , ( t ) y ] . k k k Njihov predviden položaj na podatkih oddaljenosti od cest označimo z r  [ ( r ) x , ( r ) y ] . Položaj točk dobimo z izračunom korelacije med podatki k k k zaznanih cest in podatki oddaljenosti od cest na izbranem območju ter izbiro najboljših korelacij, s katerimi dobimo kandidate točk. Korelacijo med izsekom podatkov zaznanih cest in podatkov oddaljenosti od cest se izračuna za vsak možen položaj znotraj iskalnega okna s središčem v rk. Ta kriterij korelacije cest predstavlja povprečno razdaljo piksla zaznanih cest do najbližjega piksla oddaljenosti od cest na položaju, ki je za vrednost Δ oddaljen od rk. Kriterij korelacije računamo z enačbo (Zaletelj in dr. 2013): S 1 C S x y r   P x( ) x, y( ) y D x( ) x, y( ) y , (1) k    t  t  k k   r   r   k x k y  NR x0 y0 kjer je NR število pikslov cest na izseku in sta sx, sy velikosti izseka v smereh x in y. 38 Slika 5: Potek modula za samodejno določanje oslonilnih točk. 39 V iskalnem oknu se pojavi več lokalnih minimumov kriterijske funkcije. Za vsako točko se izbere le nekaj položajev Δ z najboljšimi korelacijami, ki predstavljajo kandidate za končne položaje oslonilnih točk. Te dobimo z minimiziranjem razdalje, katero se izračuna s korelacijo. Korak se konča s preverjanjem rezultatov z vnaprej znanimi podatki o položaju, rotacijah in merilu. Drugi korak je namenjen le izsekom, ki se v prvem koraku (zaradi oblakov, majhnega števila cestnih izsekov, slabši detekciji ipd.) niso pravilno umestili. Te izseke se popravlja z rezultati sosednjih izsekov, ki imajo podobne transformacijske parametre. Parametre nepravilnih izsekov se nato ponovno preračuna in preveri njihovo natančnost. Če tudi novi parametri ne zadoščajo določeni natančnosti, se jih opusti in se jih za določanje točk ne uporablja. + Slika 6: Položaj samodejno določene oslonilne točke. Točka se nahaja na križišču izseka zaznanih cest (siva), v ozadju pa je prikazan izsek podatkov oddaljenosti od cest. Naslednji in zadnji korak je namenjen le izsekom, ki izpolnjujejo izbrano mejo natančnosti. Najprej se glede na transformacijske parametre na vsakem izseku določi izbrano število kandidatov za oslonilne točke. Točke se izbere le na križiščih cest, katerih položaj se lahko zanesljivo in nedvoumno določi in katere se ponavadi uporablja tudi pri ročnem merjenju točk. Nato se končno lokacijo oslonilnih točk določi z minimizacijo oddaljenosti zaznanih cestnih 40 pikslov od dejanske mreže cest na manjšem iskalnem območju. Slika 6 prikazuje primer končnega položaja oslonilne točke na križišču zaznanih cest. Končni izbor primernih oslonilnih točk se izvede s preverjanjem natančnosti s kriterijem povprečne oddaljenosti cest (oddaljenost do najbližje ceste glede na povprečje vseh zaznanih cestnih pikslov). Z določitvijo primernega kriterijskega praga se običajno izloči večina slabih točk z grobimi napakami; ostanejo le zelo kakovostne točke ali točke z majhnimi napakami. Večja kot je prostorska ločljivost posnetkov, več detajlov je prisotnih na cestah, ki so pogosto široke tudi več pikslov (širina cest je lahko tudi več deset pikslov). Zato je samo uporaba metode z digitalnimi podatki o cestah premalo za natančno določanje oslonilnih točk in smo ji dodali še drugo metodo korelacije (ujemanja), ki uporablja državni ortofoto (DOF) in še izboljša rezultate. Ta dodatna metoda se izvaja le pri zelo visokoločljivih posnetkih in deluje na osnovi ujemanja med satelitskim posnetkom in državnim ortofotom. Uporabljen državni ortofoto ima ločljivost 0,5 m in ocenjeno položajno točnost 1 m (GURS 2017a). Na nerazgibanem terenu je točnost boljša, medtem ko je na drugih območjih lahko odstopanje tudi več metrov. Poznamo več metod slikovne korelacije oz. slikovnega ujemanja. Uporabili smo metodo slikovnega ujemanja najmanjših kvadratov (angl. least squares matching, LSM), ki modelira ujemanje med dvema izsekoma posnetkov s šest parametrično (afino) transformacijo. Radiometrične vrednosti izseka prvega posnetka privzamemo kot opazovanja, medtem ko radiometrične vrednosti izseka drugega posnetka kot konstante. Postopek slikovnega ujemanja zapišemo kot iskanje rešitve enačbe: g(x, y) = k1  h(x', y') + k2 , (2) kjer sta koordinatna sistema surove h(x', y') in referenčna slike g(x, y) povezana z afino transformacijo: x' = a1 x + a2 y + a3 , (3) y' = b1 x + b2 y + b3 , (4) k1 in k2 pa sta radiometrična parametra za vrednosti kontrasta in svetlosti. 41 Metoda temelji na izravnavi po metodi najmanjših kvadratov, ki je prikazana v podpoglavju 5.3.1. Raziskave so pokazale, da slikovno ujemanje najmanjših kvadratov konvergira le, če je že začetno ujemanje med izsekoma slik znotraj nekaj pikslov (Schenk 1999). Zaradi tega se metodo LSM uporablja le za izboljšanje rezultatov, ne pa za korelacijo med surovimi podatki. Če uporabljamo surove podatke, je slikovno ujemanje uspešno le v primeru uporabe slikovnih piramid. V našem primeru imamo približno ujemanje zagotovljeno že s korelacijo z digitalnimi podatki o cestah, zaradi česar je metoda LSM zelo primerna za naše potrebe. Po končanem zadnjem koraku se izbrane oslonilne točke zapišejo v tekstovno datoteko. Vsaka točka v seznamu je sestavljena iz oznake, pikselskih koordinat, projekcijskih (prostorskih) koordinat in ocene točnosti, ki temelji na kriteriju povprečne oddaljenosti cest. Približno število želenih oslonilnih točk se sicer poda že pred začetkom obdelave, vendar je končno število odvisno predvsem od količine oblakov in prisotnosti ter razporeditve cestnega omrežja. 4.2 Implementacija Metoda za samodejno določanje oslonilnih točk je bila implementirana v samostojen procesni modul znotraj verige. Čeprav je bila samodejna določitev točk namensko izdelana za ortorektifikacijski postopek, večina tega dela ni bila narejena v okviru raziskave, zato bo njena implementacija prikazana le na kratko. Podrobneje bo v naslednjem podpoglavju predstavljen modul za branje in pripravo podatkov, ki v postopku ortorektifikacije nastopi prvi in je osnova za vse sledeče module. Ker se izvede pred modulom za samodejno določanje oslonilnih točk, je njegov opis umeščen v tem poglavju. 42 4.2.1 Modul za branje in pripravo metapodatkov Modul v dveh zaporednih korakih prebere metapodatke in jih pripravi za vstop v geometrični model. 4.2.1.1 Branje metapodatkov Pred vstopom posnetka v geometrični model je potrebno pripraviti metapodatke, ki so nujni za delovanje modela. V ta namen je bil izdelan modul, ki prebere metapodatke in izloči tiste parametre, ki se kasneje uporabijo pri določevanju oslonilnih točk in v postopku za ortorektifikacijo posnetka: podatke o geometriji senzorja, lastnostih posnetka ter začetne približke parametrov zunanje in notranje orientacije. Trenutno ponudniki posnetkov metapodatke večinoma že dostavljajo v standardiziranem zapisu XML (angl. EXtensible Markup Language), vendar je sama struktura datoteke XML različna (nestandardna) za skoraj vsak satelit. Izjeme so sateliti, ki so proizvod istega podjetja, vendar tudi v tem primeru lahko pride do razhajanj zaradi različnih verzij. Zaradi tega je bil bralnik za izločanje podatkov (angl. parser) izdelan za vsako vrsto posnetkov (metapodatkov) posebej. V nekaterih primerih vseh potrebnih podatkov ni v priloženi metapodatkovni datoteki (ponavadi so to podatki o snemalnem sistemu), zaradi česar je potrebno manjkajoče metapodatke pridobili iz drugih virov (medmrežje, lastnik satelita ipd.). Metapodatki, ki smo jih uporabili pri ortorektifikaciji preizkušenih posnetkov, so: - število vrstic posnetka, - število stolpcev posnetka, - število kanalov posnetka, - goriščna razdalja snemalnega sistema, - velikost piksla (elementa) v detektorju, - kot gledanja satelita v času zajema, - datum in ura zajema, - čas zajema prve vrstice, - rotacije satelita, - položaj satelita (v sistemu ECEF; glej str. 44) in - hitrost satelita (v sistemu ECEF). 43 Edina večja razlika med metapodatki posnetkov RapidEye in WorldView-2 je vrsta rotacij. Posnetki RapidEye imajo rotacije podane v obliki trojice kotov nagib, naklon in zasuk (angl. roll, pitch, yaw), medtem ko so pri WorldView-2 rotacije podane v obliki kvaternionov (glej str. 46). Bralnik metapodatkov je prirejen za zapis XML, kar pomeni, da zaporedno bere podatke iz metapodatkovne datoteke glede na najdeno značko XML (angl. tag). Ko podatek najde, ga shrani v zanj določeno spremenljivko. Vse dobljene spremenljivke se shranijo v podatkovno strukturo (angl. data structure), ki služi za enostavno izmenjavo podatkov med moduli. Podatkovna struktura je vstopni podatek za naslednji korak, ki prebrane podatke pripravi za geometrični model. 4.2.1.2 Priprava metapodatkov za vstop v geometrični model Drugi korak modula služi za preračun v prvem koraku prebranih podatkov v koordinatni sistem, ki se uporablja v modelu in v katerem bo tudi končni ortofoto. Pridobljeni podatki bodo vstopili v model kot približne vrednosti iskanih parametrov (neznanke), ki se dokončno določijo v modelu. Ker je model generičen oz. splošen, so približne vrednosti vedno enake (istega tipa) za vse senzorje. V pripravo metapodatkov vstopijo trirazsežne položajne koordinate satelita (tirnica), ki so v večini primerov že podane v koordinatnem sistemu ECEF ( angl. Earth Centered Earth Fixed) oziroma v njegovi realizaciji – koordinatnem sistemu WGS84. Koordinate so podane v časovnem intervalu, ki pokrije čas zajema posnetka ter čas pred začetkom in po koncu snemanja. V istem časovnem intervalu so ponavadi podane tudi rotacije satelita. V primeru rotacij je med senzorji večja razlika, saj so včasih podane s koti nagib, naklon in zasuk, drugič z rotacijami v sistem ECEF itn. V zadnjem času je tendenca podajanja kotov v obliki kvaternionov, ki zahtevajo še dodatno preračunavanje za vstop v model. V primeru vrstičnega skenerja se vhodni metapodatki uporabijo za izračun 24 parametrov zunanje orientacije. Ker je tirnica razdeljena na dva segmenta, se parametri dvakrat ponovijo. Vsak segment ima devet parametrov za opis položaja in tri za rotacije satelita. Pri polnoslikovnem senzorju pa imamo le tri parametre za položaj in tri za rotacije satelita (poleg parametrov notranje 44 orientacije). Našteti parametri so torej izhodne vrednosti modula. Ker ima vsak senzor svojo obliko zapisa metapodatkov, je bilo potrebno izdelati prav toliko modulov za branje in pripravo metapodatkov. Pri pripravi se metapodatkom najprej uskladi enote (pretvorba v metre, decimalne stopinje, sekunde). Pri izdelavi modula je bilo potrebno biti posebno pozoren pri enotah (tudi predznakih) za goriščno razdaljo teleskopa in velikosti elementov (pikslov) na detektorju. V nadaljevanju začne modul s preračunavanjem podatkov položaja satelita v slovenski koordinatni sistem in obliko, ki jo zahteva geometrični model. Najprej pretvori čas prve posnete vrstice in čase za podane položaje satelita v sekunde. Nato se izračuna čas začetka tzač, sredine tsre in konca tkon za oba segmenta tirnice, ki jo je satelit prepotoval v času snemanja. Šele potem je na vrsti preračun položajnih koordinat iz ECEF v slovenski koordinatni sistem D48 oz. Gauss-Kruegerjevo projekcijo. Koordinate ECEF, ki se nanašajo na elipsoid WGS84, se najprej spremenijo v kartezične koordinate na elipsoidu Bessel2. Sledi transformacija v elipsoidne koordinate in v končno projekcijo (Borčić 1976). Tako izračunane koordinate so vstopni podatek za izračun približnih vrednosti. V kolikor se ortofoti računajo za druge države oziroma v drugi projekciji, je potrebno transformacijo prilagoditi za to projekcijo. V bolj splošnem primeru se lahko izravnavo računa tudi v sistemu ECEF (npr. Poli 2005), vendar je potrebno pred uporabo vse ortofote še preprojicirati v želeni koordinatni sistem. V vseh primerih je nadaljnja obdelava enaka. Vrednosti se v modulu najprej interpolirajo s kubičnimi zlepki, kar daje najboljše rezultate, posebno, če so podatki o tirnici bolj redki (z večjimi časovnimi presledki). Iz zlepkov modul poišče vrednosti za čas začetka, sredine in konca za oba segmenta tirnice. Na koncu se iskane vrednosti dobijo z enačbami (Poli 2005): X  t  X t kon   zač  d  1 t  t kon zač 2 Obdelava lahko deluje tako v starem koordinatnem sistemu D48, kot v novem D96. Trenutno je operativna različica še v sistemu D48. 45 2 X  t  X t  X t sre   zač   kon d  2  t  t kon zač 2 X  X (5) 0 zač d 2 X  2 2 X  d  2  X . 1 1 2 V enačbah (5) je podan primer za koordinato X. Začetne vrednosti za koordinati Y in Z se dobi na enak način. Pri rotacijah iščemo le tri t.i. fotogrametrične kote ( ω, φ, κ) za vsak segment. Koti so rotacije slikovne ravnine snemalnega sistema oz. senzorja satelita v objektni koordinatni sistem. Ker v modulu iščemo samo približne vrednosti rotacij, smo izračun poenostavili tako, da smo računali rotacije v lokalni sistem, ki je tangenten na elipsoid in ima izhodišče v središču posnetka. Pri izračunu rotacij so se pokazale večje razlike pri podanih metapodatkih za posamezen satelit. Na primer, posnetki RapidEye imajo podane kote oz. rotacije glede na tirnico satelita, medtem ko imajo posnetki satelitov WorldView-2 ter Pleiades 1A in 1B rotacije glede na sistem ECEF. Poleg tega imajo nekateri novejši sateliti, tudi zadnji trije, kote podane v obliki kvaternionov. Kvaternioni so alternativen matematični zapis (v obliki kompleksnega števila) orientacij ali rotacij nekega objekta v trirazsežnem prostoru. V primerjavi z Eulerjevimi koti so bolj enostavni in jih je lažje izračunati. Trenutno se vedno bolj uporabljajo v robotiki, računalniški grafiki, navigaciji, orbitalni mehaniki satelitov itn. (Kuipers 2002). Čeprav imajo kvaternioni nekaj več prednosti, smo jih morali pred vključitvijo v model pretvoriti v Eulerjeve kote, ki so še vedno večkrat prisotni v metapodatkih in so tudi vhodni podatki v izravnavo. Transformacija poteka po enačbah (6), (7) in (8) (Airbus Defence and Space 2017). Kvaternione Qn se najprej normalizira z: 2 2 2 2 norm  Q  Q  Q  Q , (6) 0 1 2 3 46  Q 0    w  norm    Q 1   x    norm . (7)  y   Q 2       z   norm  Q 3   norm Nato se izračuna rotacijsko matriko med koordinatnim sistemom slikovne ravnine in sistemom ECEF z:  2 w  2 x  2 y  2 z 2 ( x  y  w  z) 2 ( x  z  w  y)   2 2 2 2  Rot  2 ( x y w z) w x y z 2 ( y z w x) . (8) ECEF               2( x  z  w  y) 2 ( y  z  w  2 x) w  2 x  2 y  2   z  Za izračun končne približne rotacije smo rabili še rotacijo iz sistema ECEF v lokalni tangentni koordinatni sistem z izhodiščem na površini elipsoida. Omenjena rotacija je sestavljena iz zasukov okoli osi X in Y z: cos(   90 ) 0  sin(  90  )  cos(   180 ) sin(   180 ) 0     Rot  0 1 0    . (9) LG    sin(   180 ) cos(   180 ) 0 sin(   90 ) 0 cos(   90 )   0 0      1 Končno rotacijo pa pridobimo s produktom obeh rotacij: Rot  Rot Rot (10) LG ECEF Če končno rotacijsko matriko Rot predstavimo kot:  r r r 11 12 13    Rot   r r r , (11) 21 22 23   r r r   31 32 33  dobimo Eulerjeve kote z enačbami (Kraus 1993): 47  r    atan 23    r 33    asin  r (12) 13   r    atan 12  .    r 11  V primeru, ko imamo podane kote za nagib, naklon in zasuk (npr. RapidEye), se postopek prične z interpolacijo kotov s kubičnimi zlepki. Po interpolaciji se poišče kote za sredino obeh segmentov. Rotacijo v sistem ECEF dobimo z množenjem rotacije iz sistema senzorja v orbitalni koordinatni sistem in iz orbitalnega v sistem ECEF. Rotacijo iz koordinatnega sistema satelita v orbitalni sistem dobimo z enačbo, kjer so vstopni podatki koti nagib, naklon in zasuk:  r r r 11 12 13  R    r r r , (13) SO  21 22 23   r r r 31 32 33  kjer so elementi matrike: r  cos( pitch)cos( yaw)  sin( pitch)sin( roll)sin( yaw) 11 r  cos( roll)sin( yaw) 12 r  sin( pitch)cos( yaw)  cos( pitch)sin( roll)sin( yaw) 13 r  cos( pitch)sin( yaw)  sin( pitch)sin( roll)cos( yaw) 21 r  cos( roll)cos( yaw) (14) 22 r  sin( pitch)sin( yaw)  cos( pitch)sin( roll)cos( yaw) 23 r  sin( pitch)cos( roll) 31 r  sin( roll) 32 r  cos( roll)cos( pitch). 33 Rotacijo v ECEF pa dobimo z vektorji položaja p in hitrostmi h v sistemu ECEF, ki so del metapodatkov. Vektorja najprej interpoliramo s kubičnimi zlepki in določimo vrednosti na sredini segmentov. Z enačbami (Jovanovic in dr. 1999): 48 p r   1 p r  h 1 r   (15) 2 r  h 1 r  r  r 3 2 1 izračunamo vektorje s katerimi dobimo rotacijsko matriko: R  r r r . (16) O E   3 2 1 Dobljena skupna rotacijska matrika se uporabi za rotacijo v lokalni koordinatni sistem. Na koncu se izvede še izračun končnih kotov po enačbah (12). Pri delu s kvaternioni se interpolacija izvede šele po preračunu v kote ω, φ in κ. Modul tudi v tem primeru uporablja kubične zlepke. Vstopni podatki v model so koti za sredine segmentov. Vse pripravljene podatke in druge uporabne podatke, ki smo jih prebrali iz metapodatkovne datoteke, se na koncu modula zapiše v podatkovno strukturo. Podatkovna struktura se v postopku uporablja za prenos podatkov med moduli, ki so napisani v jeziku IDL. Podatki iz podatkovne strukture se uporabijo v modulu za izračun parametrov geometričnega modela. 4.2.2 Modul za samodejno določanje oslonilnih točk Tudi modul za samodejno določanje oslonilnih točk je razdeljen na dva zaporedna koraka. Prvi korak poleg surovega posnetka uporablja še metapodatke in digitalne podatke o cestah. Digitalni podatki o cestah so na podatkovnem strežniku v rastrirani obliki in so sestavljeni iz treh virov vektorskih podatkov: cestnih podatkov iz DTK5, cestnih podatkov iz ZKGJI in podatkov OpenStreetMap za območja izven meja Slovenije. Glede na podatka o položaju in velikosti posnetka, ki ju dobimo iz metapodatkov, se izreže ustrezno območje digitalnih podatkov o cestah, iz katerih se takoj izračuna podatke oddaljenosti od cest. Za tem ali vzporedno lahko že poteka zaznava in izločanje cest iz surovega posnetka z morfološkim operatorjem. Pridobljene podatke se uporabi za določanje oslonilnih točk, kot je opisano v 49 poglavju 4. V primeru dela z visokoločljivimi posnetki se postopek po tem koraku ustavi in točke se zapišejo v tekstovno datoteko, ki je rezultat modula. Če pa ima vhodni posnetek zelo visoko ločljivost, se izvede še drugi korak. Drugi korak uporablja državni ortofoto, s katerim dodatno izboljšamo točnost oslonilnih točk. Modul za vsako točko, ki jo dobimo po prvem koraku, izreže izsek iz surovega posnetka in ortofota. Središče izsekov je iskana točka. Oba izseka se najprej prevzorčita na enako prostorsko ločljivost in obrežeta tako, da je izsek surovega posnetka velikosti 100  100 pikslov in za 20 % večji od ortofota. Najboljši položaj (največja korelacija) korelacijskega okna (izsek ortofota) znotraj iskalnega okna (izsek iz posnetka) predstavlja končni položaj točke. Pridobljene koordinate in natančnosti korelacije se zapišejo v tekstovno datoteko, ki je vhodni podatek za geometrični model. 4.3 Testiranje Uspešnost in robustnost delovanja razvitega modula za samodejno določanje oslonilnih točk smo preverili z uporabo empiričnega poskusa na podatkih, ki so bili predstavljeni v poglavju 3. 4.3.1 Ocena točnosti samodejno določenih oslonilnih točk Pred preverjanjem delovanja geometričnega modela smo preverili kakovost samodejno določenih oslonilnih točk. Ker so lahko v točkah prisotne grobe napake, je nujno določiti približno število in velikost napak. Le tako lahko kasneje pravilno ocenimo kakovost modela in ortorektifikacije. Zaradi velikega števila pridobljenih točk je vizualna kontrola zelo zamudna in nepraktična. Zato smo natančnost točk preverjali z racionalnimi polinomskimi koeficienti RPC, ki jih ponudniki posnetkov priložijo ostalim metapodatkom. Koeficienti se ponavadi uporabljajo za ortorektifikacijo, saj predstavljajo empirični model RFM za preslikavo točk iz slikovnega sistema v prostorski ali obratno. Pred uporabo je potrebno koeficiente popraviti z nekaj dobrimi oslonilnimi točkami. V praksi so popravki v obliki translacij dvorazsežnih (2R) koordinat (x in y). 50 Pred uporabo racionalnih polinomskih koeficientov RPC smo le-te preverili z ročno izmerjenimi in dobro razporejenimi kontrolnimi točkami (KT). Na vsakem posnetku smo izmerili pikselske koordinate približno 30 enakomerno porazdeljenih točk, ki so ležale na dobro vidnih križiščih cest. Iz državnega (letalskega) ortofota z ločljivostjo 0,5 m smo pridobili položajne koordinate istih lokacij. S pridobljenimi koordinatami točk smo preverili natančnost racionalnih polinomskih koeficientov RPC in dobljene rezultate uporabili za ovrednotenje poskusov s samodejno določenimi točkami. Čeprav so koeficienti običajno zelo natančni (Grodecki in Dial 2003), je njihov vpliv na določitev natančnosti točk vseeno vprašljiv na posnetkih z zelo razgibanim terenom. Natančnost točk smo preverili s primerjavo koordinat, pridobljenih prek koeficientov, in izmerjenih (določenih) koordinat. Pri tem smo za mero uporabili koren srednje kvadratne napake (angl. root mean square error, RMSE). 4.3.1.1 Poskusi s posnetki RapidEye Pred oceno točnosti samodejno določenih točk smo preverili točnost racionalnih polinomskih koeficientov RPC z ročno določenimi točkami, ki smo jih v izravnavi uporabili kot kontrolne točke. Rezultati so prikazani v preglednici 5. Če upoštevamo napako merjenja točk, ki smo jo ocenili na 0,5 piksla, lahko iz rezultatov vidimo, da so racionalni polinomski koeficienti RPC posnetkov Koper in Radgona zelo točni. Z racionalnimi polinomskimi koeficienti RPC posnetka Bohinj dobimo RMSE nad enim pikslom, kar smo pričakovali glede na topografijo območja. Na slabšo točnost racionalnih polinomskih koeficientov RPC tega posnetka kažeta tudi vrednosti največje napake in odstotek napak nad 2 σ. Preglednica 5: Točnost ročno določenih oslonilnih točk glede na koordinate izračunane iz racionalnih polinomskih koeficientov RPC za posnetke RapidEye. Število Največja Skupni Napake RMSE X RMSE Y Posnetek oslonilnih napaka RMSE nad 2 σ [%] [piksel] [piksel] točk [piksel] [piksel] Radgona 22 1,15 0 0,45 0,49 0,66 Koper 26 1,61 0 0,37 0,86 0,94 Bohinj 27 2,26 7,4 0,75 0,90 1,17 51 Samodejno določene točke smo verificirali z istimi racionalnimi polinomskimi koeficienti RPC. Za vsak posnetek smo preverili dva niza oslonilnih točk, pri čemer je prvi vseboval vse točke (večji niz), drugi pa okoli 50 točk (manjši niz). Točke v manjšem nizu smo izbrali glede na njihovo ocenjeno točnost in razporejenost po posnetku (poglavje 4). Preglednica 6 prikazuje RMSE samodejno detektiranih oslonilnih točk. Tudi rezultati teh točk kažejo na dobro ujemanje z racionalnimi polinomskimi koeficienti RPC. Za posnetka Radgona in Koper je RMSE pod enim pikslom, za posnetek Bohinj pa napaka presega en piksel. Glede na razgibano območje, ki ga pokriva posnetek Bohinj, so rezultati pričakovani in v vseh primerih kažejo na kakovostne oslonilne točke. Zanimiva je primerjava med nizi točk, saj v dveh primerih manjši niz vsebuje slabše točke, kar kaže na nezanesljivost podatka o točnosti točk, ki ga dobimo v postopku določanja. Tudi pri ohranjanju optimalne razporeditve točk se v nizu obdržijo točke z večjimi napakami. Preglednica 6: Točnost samodejno določenih oslonilnih točk glede na koordinate izračunane iz racionalnih polinomskih koeficientov RPC za posnetke RapidEye. Število Največja Napake Skupni RMSE X RMSE Y Posnetek Niz oslonilnih napaka nad 2 σ RMSE [piksel] [piksel] točk [piksel] [%] [piksel] Večji 501 5,86 1,6 0,51 0,41 0,65 Radgona Manjši 50 5,68 4,0 0,79 0,55 0,96 Večji 467 3,46 3,9 0,48 0,73 0,88 Koper Manjši 56 1,70 5,3 0,47 0,77 0,90 Večji 325 3,63 4,3 0,75 1,24 1,45 Bohinj Manjši 53 3,02 5,7 0,69 1,08 1,28 Za pravilno validacijo celotnega ortorektifikacijskega postopka je pomembno oceniti velikost odstopanj na določenih oslonilnih točkah. Dobra kazalnika za validacijo in prisotnost grobih napak sta največja napaka in odstotek napak, ki so za 2 σ večji ali manjši od srednje napake. Čeprav lahko iz preglednice 6 sklepamo, da so v vseh nizih prisotne večje grobe napake, je njihov odstotek majhen in nikoli ne presega 6 %. Izračunane vrednosti potrjujejo tezo, da v geometričnem modelu potrebujemo dobre metode za odkrivanje in izločanje grobih napak. 52 4.3.1.2 Poskusi s posnetki WorldView-2 Zaradi zelo visoke ločljivosti posnetkov WorldView-2 so bile oslonilne točke pridobljene v dveh korakih (glej poglavje 4.1). Njihovo točnost pa smo ocenili z enakim postopkom, kot pri posnetkih RapidEye. Ocene točnosti racionalnih polinomskih koeficientov RPC z ročno določenimi točkami, ki smo jih v izravnavi uporabili kot kontrolne točke, so prikazane v preglednici 7. Če upoštevamo napako merjenja točk, je iz rezultatov razvidno, da so racionalni polinomski koeficienti RPC zelo točni, najbolj za posnetka Jesenice in Ljubljana. Koeficienti posnetka Jesenice so bili preverjeni samo na spodnji polovici posnetka, ki leži v Sloveniji. Racionalni polinomski koeficienti RPC posnetka Bled so nekoliko slabši, se pa glede na RMSE bistveno ne razlikujejo od ostalih dveh posnetkov. Nekoliko slabša točnost racionalnih polinomskih koeficientov RPC tega posnetka se kaže tudi v vrednosti največje napake in odstotku napak nad 2 σ. Preglednica 7: Točnost ročno določenih oslonilnih točk glede na koordinate izračunane iz racionalnih polinomskih koeficientov RPC za posnetke WorldView-2. Število Največja Skupni Napake RMSE X RMSE Y Posnetek oslonilnih napaka RMSE nad 2 σ [%] [piksel] [piksel] točk [piksel] [piksel] Jesenice 24 1,41 4,2 0,56 0,29 0,63 Bled 30 1,93 6,7 0,70 0,63 0,94 Ljubljana 30 1,64 3,3 0,52 0,55 0,76 Zaradi manjšega števila cest oz. križišč na posnetkih, lege (nekateri deli posnetkov so izven Slovenije) in razgibanega terena, je bilo na posnetkih WorldView-2 samodejno določenih manj točk. Posledično je bilo na voljo tudi manj nizov za nadaljnje poskuse. Samodejno določene točke smo najprej preverili z racionalnimi polinomskimi koeficienti RPC. Za vsak posnetek smo izbrali dva niza oslonilnih točk, pri čemer je prvi vseboval vse točke (večji niz), drugi pa okoli 50 oz. 30 točk (manjši niz). Točke v manjšem nizu smo izbrali glede na njihovo natančnost in razporejenost po posnetku. Preglednica 8 prikazuje dobljene RMSE samodejno določenih oslonilnih točk. Točnosti glede na ročno določene točke so nekoliko slabše, posebno za posnetka Jesenice in Ljubljana, v vseh primerih pa nekoliko presegajo velikost enega 53 piksla. Najslabši rezultati so v primeru posnetka Jesenice, kjer smo (predvsem zaradi lege posnetka) dobili najmanj točk. Preglednica 8: Točnost samodejno določenih oslonilnih točk glede na koordinate izračunane iz racionalnih polinomskih koeficientov RPC za posnetke WorldView-2. Število Največja Napake Skupni RMSE X RMSE Y Posnetek Niz oslonilnih napaka nad 2 σ RMSE [piksel] [piksel] točk [piksel] [%] [piksel] Večji 75 5,05 2,67 0,87 0,90 1,25 Jesenice Manjši 31 4,82 3,23 0,86 1,19 1,47 Večji 158 3,73 4,43 0,83 0,74 1,11 Bled Manjši 50 3,52 6,0 0,88 0,91 1,26 Večji 206 5,25 3,88 0,76 0,86 1,15 Ljubljana Manjši 48 2,17 6,25 0,84 0,63 1,05 Če primerjamo nize točk, ugotovimo, da v dveh primerih manjši niz vsebuje slabše točke, kar tudi pri tem senzorju kaže na nezanesljivost podatka o točnosti točk, ki ga dobimo v postopku določanja. Tudi pri ohranjanju optimalne razporeditve točk se v nizu obdržijo točke z večjimi napakami. Čeprav lahko iz preglednice 8 sklepamo, da so v vseh nizih prisotne večje grobe napake, je njihov odstotek majhen in le redko preseže 6 %. Kljub temu je točnost določanja nekoliko slabša kot v primeru posnetkov RapidEye in zaradi tega lahko pričakujemo tudi slabše rezultate geometričnega modela. 54 5 GEOMETRIČNI MODELI SENZORJEV Za izdelavo ortofotov z visoko položajno točnostjo je potrebno uporabiti geometrični model, ki čimbolj natančno opiše geometrijo zajema posnetka. Geometrija zajema je odvisna od izvedbe snemalnega sistema; različne izvedbe so opisane v prvem podpoglavju. Sledi opis sistematičnih napak oslonilnih točk, ki so posledica atmosferske refrakcije in ukrivljenosti Zemlje. V tretjem in četrtem podpoglavju sta podrobno opisani izvedbi geometričnega modela za vrstične in polnoslikovne snemalne sisteme, ki sta bili izdelani v sklopu raziskave. Prikazani so sestavni deli geometričnega modela, ki vključuje tudi samodejno izločanje grobih napak v dveh korakih. Razvitih je bilo že več različnih geometričnih modelov za različne satelite (podrobnejši opis v podpoglavju 1.1.1). Nekateri geometrični modeli so bili narejeni in prilagojeni le za posamezne satelite, veliko pa je bilo poskusov izdelave splošnega geometričnega modela za vse primere, ki uporabljajo določen tip senzorja (npr. vrstični, polnoslikovni, točkovni). Geometrični modeli se delijo glede na uporabljen tip senzorja in glede na vrsto uporabljenih podatkov. T.i. direktni (neposredni) modeli izračunajo parametre zunanje orientacije samo iz podatkov predvidene tirnice satelita in vrednosti meritev položaja in orientacije satelita z napravami GNSS in sledilci zvezd. Vrednosti meritev se pred izračunom parametrov interpolirajo na čas snemanja vrstic pri vrstičnem senzorju ali celotnega posnetka pri polnoslikovnem senzorju. Sam izračun je hiter, vendar zaradi netočnih meritev rezultati še ne dosegajo želenih natančnosti. Drugi tip geometričnih modelov so t.i. indirektni (posredni) modeli, ki delujejo z oslonilnimi točkami in so še vedno največkrat uporabljen način izračuna parametrov zunanje orientacije posnetkov. Čeprav je iskanje oslonilnih točk lahko dolgotrajen postopek in izračun modela zahtevnejši, so rezultati zelo natančni in primerni 55 tudi za najzahtevnejše namene uporabe. Z indirektnimi modeli je mogoče tudi zelo natančno določanje in popravljanje parametrov notranje orientacije (npr. goriščna razdalja, položaj in oblika detektorja ipd.), ki pa so ponavadi že zelo natančno kalibrirani (v laboratoriju in tudi že v tirnici – angl. in-flight) in jih običajno dobimo v metapodatkih posnetka. Knjiga se osredotoča na dve izvedbi geometričnega modela, ki sta bili izdelani za delo z vrstičnimi ali polnoslikovnimi senzorji. Samo jedro delovanja modela je bilo izbrano iz množice že razvitih modelov. Končna izbira je bila narejena na podlagi stopnje splošnosti, enostavnosti in doseženih točnosti modela. Izbrano jedro modela smo nato dopolnili in mu dodali funkcionalnosti za obdelavo podatkov v samodejnem postopku ortorektifikacije. Obe izvedbi modela sta indirektni in splošni, kar pomeni, da delujeta le z oslonilnimi točkami in sta prilagojeni za modeliranje geometrije vseh snemalnih sistemov, ki snemajo v vrstičnem ali polnoslikovnem načinu. 5.1 Senzorji na satelitih za opazovanje Zemlje Večina satelitov za opazovanje Zemlje uporablja eno od treh vrst snemalnih sistemov: prečni (angl. whiskbroom), vzdolžni (angl. pushbroom) in polnoslikovni (angl. fullframe). Vsaka vrsta snemalnega sistema pa za svoje delovanje potrebuje določeno obliko detektorja (CCD ali CMOS). 5.1.1 Točkovni senzorji Prečni snemalni sistemi oz. t.i. prečni skenerji Zemljo snemajo z zaporedji vrstic, ki so pravokotne na smer gibanja platforme. Vsako vrstico zabeležijo s premikanjem senzorja z ene strani na drugo, pri čemer uporabljajo nihajoče ali vrteče se zrcalo. Ker se nosilna platforma (satelit) premika naprej, lahko z zaporednimi vrsticami sestavimo dvorazsežni posnetek površine Zemlje (Oštir 2006). Prečni snemalni sistemi uporabljajo točkovni (angl. point-based) detektor, ki ima malo elementov in jih je zaradi tega tudi lažje kalibrirati in vzdrževati. Glavna omejitev prečnih sistemov je omejeni čas za branje detektorja, saj se ogledalo zelo hitro premika. Za dosego zadostnega razmerja signal-šum (angl. signal-to-noise ratio) imajo taki sistemi precej široke 56 spektralne kanale, zaradi kompleksnosti sistema in veliko premikajočih se delov pa so sateliti, opremljeni s prečnimi skenerji, večji in težji. Znani snemalni sistemi, ki uporabljajo (ali so uporabljali) točkovne senzorje, so MSS, TM in ETM+ na satelitih Landsat, AVHRR na satelitih POES in SeaWiFS na satelitu SeaStar. Zaradi omenjenih omejitev se za opazovanje Zemlje v visoki prostorski ločljivosti prečnih snemalnih sistemov ne uporablja. 5.1.2 Vrstični senzorji Vzdolžni sistemi oz. skenerji uporabljajo vrstične detektorje, ki so sestavljeni iz vrste elementov, ki se nahajajo v goriščni razdalji snemalnega sistema. Ker je majhne detektorje (nekaj cm) za vesoljske potrebe težko izdelati, je število elementov v vrsti omejeno. Zaradi tega je lahko na goriščni ravnini tudi več detektorjev z več tisoč elementi, ki so lahko razporejeni tudi v več nivojih. Dobljene posnetke posameznih detektorjev se kasneje združi v "umetni posnetek", ki je osnova za posnetek ravni obdelave 1A. Vzdolžni snemalni sistemi naenkrat zajamejo celotno vrstico na Zemlji in z izkoriščanjem gibanja snemalnega sistema ustvarijo posnetek v smeri gibanja (slika 7). Ker njihovo delovanje spominja na pometanje z metlo, vzdolžne skenerje včasih imenujemo tudi pometajoči. Vsak kanal potrebuje poseben sistem detektorjev. Ti zaznavajo svetlobo (energijo) po vrsticah, jo nato spremenijo v električne impulze in pretvorijo v digitalno obliko (Oštir 2006). Registriranje količine naboja in "praznjenje" detektorja mora biti dovolj hitro, da dovoli neodvisno sprejemanje svetlobe za naslednjo vrstico. Čas, ki ga porabi sistem za snemanje ene vrstice, imenujemo čas integracije (angl. integration time). Ker je velikost elementov na novejših detektorjih že zelo majhna, čas integracije pa kratek, detektorji zaznajo zelo malo svetlobe. Če občutljivost detektorjev ni dovolj velika za doseganje potrebnega razmerja signal-šum, obstajate dve možnosti za premostitve te ovire (Sandau 2004): - uporaba tehnologije TDI (angl. Time Delay and Integration), ki z večkratnim zajemom istih objektov zagotavlja dobro kakovost posnetka s hitrejšim zajemom vrstic, 57 - ali uporaba t.i. načina upočasnjevanja (angl. slow-down mode), ki zmanjša hitrost premikanja projekcije vrstice na površju glede na hitrost satelita. Zaradi enostavnosti in zanesljivosti se v sodobnih satelitih večinoma uporablja tehnologija TDI, ki temelji na večkratnem zajemu iste vrstice oz. z zajemom v več stopnjah (angl. TDI stages). Vsaka stopnja podaljša čas zaznavanja svetlobe na detektorju. Od števila stopenj je odvisna tudi izboljšava razmerja signal-šum. Slika 7: Zajem posnetka z vrstičnim sistemom. Vzdolžni sistemi imajo več prednosti v primerjavi s prečnimi sistemi. Množica elementov detektorja omogoča vsakemu od njih, da dlje časa opazuje površino pod seboj. Tako lahko detektor zazna več svetlobe, kar omogoča manjše vidno polje (angl. field of view) in ožje spektralne pasove, zato lahko s takimi detektorji dosežemo boljšo prostorsko in spektralno ločljivost. 58 Detektorji so preproste mikroelektronske naprave, ki so majhne, lahke in zahtevajo malo energije. Poleg tega nimajo gibljivih delov, zaradi česar so bolj zanesljivi, manj občutljivi na mehanske okvare in imajo daljšo življenjsko dobo. Po drugi strani pa moramo umeriti več tisoč elementov detektorja, da lahko dosežemo enakomerno občutljivost po celotnem posnetku (Oštir 2006). Trenutno vzdolžne snemalne sisteme uporabljajo skoraj vsi boljši sateliti za opazovanje Zemlje, kot so na primer WorldView-1, WorldView-2, WorldView-3, GeoEye, Pleiades, QuickBird ipd. 5.1.3 Polnoslikovni senzorji Polnoslikovni detektorji imajo obliko pravokotne matrike, kjer je vsak element matrike tudi element detektorja. Pri snemanju vsi elementi detektorja sprejmejo svetlobo v istem trenutku. Na ta način celoten posnetek nastane v enem času integracije, ki je zelo kratek (običajno pa daljši kot pri vrstičnih senzorjih). Posnetki se zajamejo v centralni projekciji z enim samim perspektivnim centrom (slika 8). Zaradi tega je posnetek nepopačen le v nadirju, v vseh smereh od te točke pa se popačenja povečujejo z naraščanjem oddaljenosti. Popačenja so manjša, če je posnetek pridobljen iz visoke višine leta in ima ozko vidno polje. Prav te lastnosti lahko izkoriščajo sateliti. Polnoslikovni senzorji so se do pred kratkim pretežno uporabljali le pri letalskih snemanjih. Z uveljavitvijo novih in izpopolnjenih tehnologij ter miniaturizacijo komponent pa so se začeli uporabljati tudi na majhnih satelitih. Najprej se je na majhnih satelitih uporabljalo kar detektorje in optiko sodobnih digitalnih fotoaparatov, s katerimi pa se je dobilo slabšo prostorsko ločljivost. V zadnjih letih se specializirane polnoslikovne detektorje uporablja tudi v visokoločljivih snemalnih sistemih (Grocott in dr. 2013), ki dokazujejo, da se bo njihova uporaba v bodočih misijah še povečala. Glavna prednost pred drugimi vrstami snemalnih sistemov je njihova enostavnost in zmožnost vgradnje tudi v majhne satelite. Z njimi je iz vesolja možno tudi zajemanje videa. Kljub temu je za izdelavo profesionalnega sistema z (zelo) visoko prostorsko ločljivostjo in dobrim razmerjem signal- šum potreben kakovosten teleskop in natančno umerjanje detektorja. 59 Primer satelita s polnoslikovnim sistemom je tudi prvi slovenski satelit, ki bo predvidoma izstreljen leta 2018. Ta majhen satelit bo letel na višini 600 km, imel bo maso okoli 70 kg in bo nosil dva snemalna sistema: visokoločljiv sistem z ločljivostjo 2,8 m in štirimi kanali ter srednjeločljiv sistem z ločljivostjo 40 m. Oba sistema bosta lahko snemala tudi video z ločljivostjo 1920  1080 pikslov (Grocott in dr. 2013). Slika 8: Zajem posnetka s polnoslikovnim sistemom. 60 5.2 Popravki koordinat oslonilnih točk V raziskavi izdelana geometrična modela sta indirektna, kar pomeni, da za izračun parametrov zunanje orientacije uporabljata oslonilne točke. Oslonilne točke se določijo samodejno in v model vstopijo v obliki slikovnih (pikselskih) koordinat surovega posnetka (slikovni prostor) in objektnih (prostorskih) koordinat referenčnih podatkov (objektni prostor). Za vsako točko sta torej podana dva para položajnih koordinat. Pri modeliranju geometrični model poveže položaj točke v slikovnem prostoru s homologno točko v objektnem prostoru z linearnimi povezavami kolinearnih enačb. Ker sateliti snemajo veliko območje visoko nad tlemi, zaradi ukrivljenosti Zemlje in atmosferske refrakcije uporaba linearnih povezav med točkami ni možna. Omenjena dejavnika povzročata odstopanja od matematičnega modela, ki niso naključna. Taka odstopanja spadajo med sistematične napake in se jih lahko popravi s premikom položaja točk za vrednosti obeh omenjenih vplivov. Ponavadi se popravlja le slikovni koordinati, objektni se pusti nespremenjeni. Popravki so običajno majhni, naraščajo pa z večanjem kota gledanja. Posebno pomembni so pri visokoločljivih posnetkih, kjer lahko majhno odstopanje pomeni veliko izgubo točnosti. Sistematične popravke se odpravlja v obratni smeri nastanka (Mikhail in dr. 2001). 5.2.1 Atmosferska refrakcija Žarek, ki prehaja skozi atmosfero, se giblje po Snellovem zakonu. Zaradi prehoda med plastmi zraka z različno gostoto se žarek lomi do vstopa v snemalni sistem. Zaradi tega sistematičnega vpliva je opazovana točka premaknjena glede na njen pravi položaj. Popravek zaradi atmosferske refrakcije dobimo po postopku, ki ga je opisal Noerdlinger (1999). Če točko P na površju Zemlje opazujemo iz točke D na satelitu, v resnici vidimo točko P' , ki je za razdaljo d premaknjena od njenega pravega položaja (slika 9). Razdalja d pomeni tudi popravek koordinat, ki jih dobimo iz surovega posnetka. Pri izračunu popravkov so pomembni trije zenitni koti: - zenitni kot Z0, ki je kot med navpičnico na točki P in idealnim žarkom, na katerega ne vpliva refrakcija, 61 - zenitni kot Z, ki je kot med navpičnico na točki P' in idealnim žarkom, na katerega ne vpliva refrakcija ter - zenitni kot Z' , ki je kot med normalo na točki P' in žarkom pod vplivom refrakcije. Kot Z0 dobimo iz metapodatkov o kotu gledanja satelita, za izračun popravka pa rabimo ostala dva kota. Premik zaradi refrakcije izračunamo z enačbo: d  R ( Z  Z ) , (17) 0 kjer je R radij Zemlje, razlika med zenitnimi koti pa je podana v radianih. Če predpostavimo, da je atmosferski model sferično simetričen, lahko kot Z' izračunamo kot funkcijo kota Z0:  sin( Z )  Z' = asin  0  , (18)  0  kjer je μ0 lomni količnik v točki P' . Slika 9: Geometrični prikaz atmosferske refrakcije. 62 Ker sta pri satelitskih posnetkih kota Z in Z0 zelo majhna, lahko predpostavimo, da je atmosfera v območju posnetka vzporedna z ravnim površjem in z atmosfersko višino homogene gostote (angl. homogeneous atmospheric density scale height) W=8592 m (glede na vrednosti v Allen 1973). Kot Z dobimo s kotom Z' in refrakcijo kot jo je predlagal Allen (1973):  1 Z  Z' tan Z' 0  ,00117tan Z 3 ' , (19) W 1 R kjer je W atmosferska višina homogene gostote in μ lomni količnik, ki je definiran kot:   1  1  , (20) 0  0 kjer je ρ/ρ0 razmerje med gostoto atmosfere in gostoto na morski gladini ter se ga lahko dobi iz atmosferskega modela. Vrednosti atmosferske refrakcije na morski gladini in izhajajoči linearni premiki so prikazani v preglednici 9. Preglednica 9: Vrednosti atmosferske refrakcije na morski gladini (Noerdlinger 1999). Zenitni kot v Zenitni kot na Atmosferska Linearni premik d vesolju Z0 [°] površju Z' [°] refrakcija ( Z0 - Z' ) [°] [m] 10 9,9971 0,0029 0,55 20 19,9939 0,0061 1,22 30 29,9904 0,0096 2,22 40 39,9860 0,0140 3,98 45 44,9834 0,0166 5,46 50 49,9802 0,0198 7,73 55 54,9762 0,0238 11,40 60 59,9712 0,0288 17,85 65 64,9643 0,0357 30,37 5.2.2 Ukrivljenost Zemlje Zaradi ukrivljenosti Zemlje je opazovana točka prikazana bližje nadirju kot je v resnici. Popravek ukrivljenosti Zemlje računamo po enačbi (Mikhail in dr. 2001): 63 r 3  H d  (21) E 2  f 2  R kjer je dE popravek koordinate, H je višina orbite satelita, f je goriščna razdalja, R je radij Zemlje in r je razdalja od centra posnetka do točke na detektorju. 5.3 Geometrični model vrstičnega senzorja V okviru raziskave smo najprej izdelali geometrični model za vrstične senzorje oz. snemalne sisteme. Implementiran geometrični model je vrsta matematičnega modela, ki poskuša oceniti fizično realnost pri snemanju satelitskih posnetkov. Kot vsak matematični model je razdeljen na funkcionalni in stohastični model. Funkcionalni model v splošnem določa fizične lastnosti, stohastični pa verjetnostne lastnosti. Oba modela povežemo in izvrednotimo v postopku izravnave. 5.3.1 Izravnava po metodi najmanjših kvadratov Matematični model povezuje opazovanja s količinami, ki nas zanimajo (neznanke). Za izračun zahtevanih neznank moramo določiti minimalno število spremenljivk, ki jih potrebujemo za enolično določitev modela. V primeru, da rešujemo model z minimalnim številom opazovanj, lahko napaka kateregakoli opazovanja povzroči napačne vrednosti iskanih neznank. Zaradi tega vedno uporabimo več opazovanj (nadštevilna opazovanja), kot jih je za enolično določitev modela nujno potrebnih. V primeru, ko imamo nadštevilna opazovanja, moramo za pridobitev enolične rešitve uporabiti izravnavo. Med vsemi metodami izravnav se je zaradi splošnosti in sistematičnosti najbolj uveljavila metoda najmanjših kvadratov (MNK). Metoda temelji na iskanju množice popravkov opazovanj, ki ima najmanjšo vsoto kvadratov popravkov. Izravnana opazovanja morajo poleg prileganja modelu zagotavljati tudi pogoj minimalne vsote kvadratov odgovarjajočih popravkov. Izravnavo po metodi najmanjših kvadratov lahko računamo z različnimi tehnikami. Najbolj pogosti sta tehniki pogojne in posredne izravnave. V 64 splošnem so vse tehnike izravnave po MNK enakovredne in zagotavljajo identične rezultate, kljub temu se v fotogrametriji najbolj pogosto uporablja tehnika posredne izravnave. Pri tej tehniki je število enačb popravkov enako skupnemu številu opazovanj (vsaka enačba vsebuje samo eno opazovanje). Ker je število pogojev enako številu nadštevilnih opazovanj, tehnika posredne izravnave vključuje dodatne spremenljivke oz. t.i. neznanke. Neznanke so količine, ki nimajo vnaprej znane vrednosti in se določijo šele v postopku izravnave. Splošno obliko enačb popravkov v izravnavi posrednih opazovanj po MNK lahko zapišemo v matrični obliki kot: v  B  f , (22) kjer je v vektor popravkov, B matrika koeficientov, Δ vektor neznank in f vektor opazovanj. Če označimo število opazovanj z n in neznanke z u, ima matrika B dimenzije n × u, vektor Δ dimenzije u × 1, vektorja v in f pa n ×. Zaradi različne natančnosti opazovanj je potrebno v izravnavo uvesti pojem uteži. Utež opazovanja se določi glede na njegovo natančnost, ki je ponavadi podana z vrednostjo variance 2 ali standardnega odklona . Visoka natančnost pomeni majhno varianco in obratno. Ker je vrednost variance nasprotna natančnosti, je utež opazovanja definirana kot količina, ki je obratno sorazmerna varianci: k p  , (23) 2  kjer je p utež, k pa poljubna konstanta. Za vrednost konstante se večkrat uporabi kar referenčno oz. a-priori varianco 2  . 0 Ko imamo opravka z več opazovanji, so lahko le-ta medsebojno odvisna oziroma korelirana. To odvisnost med pari opazovanj označimo s kovariancami. V primeru neodvisnih opazovanj se lahko uteži opazovanj združi v diagonalno matriko uteži opazovanj P, ki se jo nato vključi v izravnavo za končno določitev iskanih količin (neznank). Z matrično obliko enačb popravkov se vektor neznank določi z minimumom kvadratov popravkov: v T Pv(min)   f  B T P f  B (24) 65 z rešitvijo:   BT 1 -  PB  BT Pf  . (25) Enačbo (25) lahko zapišemo tudi v obliki linearnega sistema:  N t 1 -  , (26) kjer je N matrika koeficientov normalnih enačb, t pa vektor konstantnih členov z enačbama: N BT  PB in (27) t BT  Pf . (28) Pred izračuni se neznankam določi približne vrednosti, med izravnavo pa se izračuna le njihove popravke. Po prvi izravnavi se dobi le prve ocene neznank. Za pridobitev točnih vrednosti je treba postopek večkrat iterativno ponoviti. Ocenjene neznanke prejšnje iteracije se uporabi kot približke vrednosti neznank za novo iteracijo. Postopek se ponavlja, dokler popravki približnih vrednosti neznank ne izpolnijo vnaprej postavljenih mejnih vrednosti, ki se lahko nanašajo na popravke približnih vrednosti neznank, popravke opazovanj ali varianco a-posteriori. 5.3.2 Funkcionalni model V večini primerov, ki se uporabljajo v fotogrametriji, je funkcionalni model izražen s kolinearnimi enačbami. S kolinearnimi enačbami lahko povežemo slikovne koordinate točke digitalnega posnetka s homologno točko s prostorskimi koordinatami (Mikhail in dr. 2001): m  X  X  m  Y  Y  m  Z  Z 11  S  12  S  13  S  x  f m  X  X  m  Y  Y  m  Z  Z 31  S  32  S  33  S  (29) m  X  X  m  Y  Y  m  Z  Z 21  S  22  S  23  S  y  f m  X  X  m  Y  Y  m  Z  Z 31  S  32  S   S  , 33 kjer so X, Y in Z koordinate točke v prostorskem (objektnem) koordinatnem sistemu; XS, YS in ZS so koordinate perspektivnega centra v prostorskem koordinatnem sistemu; x in y sta koordinati točke v slikovnem koordinatnem 66 sistemu; mxy so elementi rotacijske matrike za rotacijo iz slikovnega v prostorski koordinatni sistem in f je goriščna razdalja. Enačbi (29) lahko zapišemo tudi kot: Nx x   f D (30) Ny y   f , D kjer so: Nx  m  X  X  m  Y  Y  m  Z  Z 11  S  12  S  13  S  Ny  m  X  X  m  Y  Y  m  Z  Z (31) 21  S  22  S  23  S  D  m  X  X  m  Y  Y  m  Z  Z 31  S  32  S   S . 33 Rotacijska matrika je oblike:  m m m 11 12 13    M   m m m , (32) 21 22 23   m m m   31 22 33  kjer imajo elementi vrednosti: m  cos( ) cos( )  sin( ) sin() sin( ) 11 m  cos() sin( ) 12 m  sin( ) cos( )  cos( ) sin() sin( ) 13 m  cos( ) sin( )  sin( ) sin() cos( ) 21 m  cos()cos( ) (33) 22 m  sin( ) sin( )  cos( ) sin() cos( ) 23 m  sin( ) cos() 31 m  sin() 32 m  cos() cos(). 33 Ker snemalni sistemi uporabljajo različne vrste senzorjev in načine snemanja, se geometrični modeli, kljub temu, da temeljijo na enakih enačbah, med seboj razlikujejo glede števila parametrov, kompleksnosti izračunov ipd. Zaradi snemanja posnetka po vrsticah, ki lahko traja več sekund, je geometrični model vrstičnega senzorja zelo zapleten. 67 V okviru raziskave je bil izdelan splošen in rigorozen geometrični model za vrstične senzorje. S pravilno predpripravo metapodatkov model deluje praktično z vsemi snemalnimi sistemi sodobnih satelitov, ki uporabljajo vrstične senzorje. Poleg tega model rigorozno (fizično) opiše razmerje med posnetkom in objektnim prostorom s kolinearnimi enačbami. Uporabljen geometrični model vrstičnega senzorja ima 24 parametrov zunanje orientacije, ki v modelu nastopajo kot neznanke. Ker imajo sodobni sateliti z vrstičnim snemanjem zelo stabilno notranjo orientacijo (Grodecki in Dial 2003), lahko parametre notranje orientacije obravnavamo kot konstante in se med izravnavo ne spreminjajo. Njihove vrednosti so rezultat natančnih kalibracij in jih dobimo v metapodatkih. Geometrični model vrstičnega senzorja je podrobno opisan v naslednjih podpoglavjih. 5.3.2.1 Modeliranje zunanje orientacije Zunanjo orientacijo senzorja, ki pove položaj in orientacijo senzorja v času snemanja, smo modelirali z odsekovnimi polinomskimi funkcijami (angl. piecewise polynomial functions), ki so odvisne od časa. Pri tem se tirnica, ki jo satelit preleti v času snemanja posnetka, razdeli na več medsebojno povezanih odsekov. V geometričnem modelu smo uporabili dva odseka (slika 10), za katera Poli (2005) navaja, da sta optimalna za dosego natančnih rezultatov za večino satelitskih sistemov z vrstičnimi senzorji. V vsakem odseku oz. segmentu se parametre zunanje orientacije modelira s šestimi enačbami, ki vsebujejo 12 parametrov. Tri enačbe, ki definirajo položaj satelita, so v obliki polinomskih funkcij drugega reda in so odvisne od časa t, ostale so konstantne. Enačbe imajo konstantne ( X0, Y0, Z 0, ω0, φ0, κ0), linearne ( X1, Y1, Z1) in kvadratne ( X2, Y2, Z2) člene: 68 X( t ) 2  X  X t  X t 0 1 2 Y( t ) 2  Y  Y t  Y t 0 1 2 Z( t) 2  Z  Z t  Z t 0 1 2 (34)    0    0    , 0 pri čemer se čas t izračuna z enačbo: t  t v zač t  , (35) t  t kon zač kjer je tkon čas na koncu segmenta, tzač čas na začetku segmenta in tv čas vrstice. Parametri orientacije se v času snemanja praktično ne spremenijo in jih lahko obravnavamo kot konstantne. V točki, kjer se segmenta spojita, je treba zagotoviti zveznost tirnice. Zvezno nadaljevanje se zagotovi z enačenjem vrednosti funkcij ter njihovih prvih in drugih odvodov. Da zagotovimo zveznost tirnice, morajo imeti funkcije na koncu prvega in na začetku drugega segmenta enako vrednost. Če enačimo funkciji segmentov 1 in 2 za X, dobimo (Poli 2005): 2 2 X 0  X 1 t  X 2 t  X 0  X 1 t  X 2 t (36) 1 1 1 2 2 2 Slika 10: Modeliranje tirnice z dvema segmentoma. V praksi je prehod med segmentoma veliko bolj gladek, kot je prikazano na sliki, kar zagotovimo z enačenjem vrednosti funkcij segmentov ter njihovih prvih in drugih odvodov. 69 Če pri tem upoštevamo, da je ob koncu segmenta 1 t=1 in na začetku segmenta 2 t=0, dobimo: X 0  X 1  X 2  X 0 (37) 1 1 1 2 Po enakem principu dobimo enačbe za vse parametre zunanje orientacije: X 0  X 1  X 2  X 0 1 1 1 2 Y 0  Y 1  Y 2  Y 0 1 1 1 2 Z 0  Z 1  Z 2  Z 0 1 1 1 2 (38) 0  0 1 2 0 0 1 2  0  0 1 2 Za gladek in zvezen prehod enega segmenta na drugega je treba zagotoviti tudi enakost prvega in drugega odvoda funkcij zunanjih parametrov. Funkcije odvajamo po času t, zaradi česar odpadejo funkcije rotacij. Če enačimo prve odvode funkcij segmentov 1 in 2 za X, dobimo: X 1  2 X 2 t  X 1  2 X 2 t (39) 1 1 2 2 Če pri tem upoštevamo, da je ob koncu segmenta 1 t=1 in na začetku segmenta 2 t=0, dobimo: X 1  2 X 2  X 1 (40) 1 1 2 Po enakem principu dobimo enačbe za vse položajne parametre zunanje orientacije: X 1  2  X 2  X 1 1 1 2 Y 1  2  Y 2  Y 1 (41) 1 1 2 Z 1  2  Z 2  Z 1 . 1 1 2 Enak postopek je tudi za drugi odvod. Če enačimo druge odvode funkcij segmentov 1 in 2 za X, dobimo: 2 X 2  2 X 2 (42) 1 2 oziroma: 70 X 2  X 2 (43) 1 2 Po enakem principu dobimo enačbe za vse položajne parametre zunanje orientacije: X 2  X 2 1 2 Y 2  Y 2 (44) 1 2 Z 2  Z 2 . 1 2 5.3.2.2 Linearizacija Kolinearne enačbe niso linearne glede na neznanke, zato jih moramo pred uporabo v modelu linearizirati. Linearizacija enačb poteka z razvojem funkcij v Taylorjevo vrsto, kjer ohranimo le prva dva člena (do prvega odvoda). Če je v funkciji n neznank, je Taylorjeva vrsta oblike (Kraus 1993):  f    f   f ( x , x ,, x )  f ( x , x ,, x )   x    x   1 2 n 1 2 n 0 1 2     x  x   1   2  x 10 x 20 (45)  f     x  . n   x   n  xn 0 V primeru kolinearnih enačb in funkcij parametrov zunanje orientacije se odvodi računajo za vse neznanke ( X0, Y0, Z0, X1, Y1, Z1, X2, Y2, Z2, ω0, φ0, κ0). V primeru neznanke X0 se odvoda po x in y računata z enačbama: x  x  X  S   (46) X  X  X  0 S 0 y  y  X  S   . (47) X  X  X  0 S 0 Podobne enačbe veljajo za vse ostale neznanke. Parcialne odvode prvega člena glede na neznanke dobimo z (Kraus 2007): x  f   ( m Nx  m D)  b 2 13 11 1 , 1 X  D S x  f   ( m Nx  m D)  b 2 23 21 , 1 2 Y  D S 71 x  f   ( m Nx  m D)  b 2 33 31 , 1 3 Z  D S x  f  Nx    (( Y  Y ) m ( Z  Z ) m ) ( Y  Y ) m ( Z  Z ) m   b S 33 S 23 S 31 S 21 , 1 4   D  D  S x  f  Nx   ( Nx cos  Ny sin)  D cos  b ,15   D  D  S x  f   Ny  b (48) , 1 6   D S y  f   ( m Ny  m D)  b 2 13 12 2 1 , X  D S y  f   ( m Ny  m D)  b 2 23 22 2,2 Y  D S y  f   ( m Ny  m D)  b 2 33 32 2,3 Z  D S y  f  Ny    (( Y  Y ) m ( Z  Z ) m ) ( Y  Y ) m ( Z  Z ) m   b S 33 S 23 S 32 S 22 2,4   D  D  S y  f  Ny   ( Nx cos  Ny sin)  D sin  b 2,5   D  D  S y  f  Nx  b . 2,6   D S Parcialni odvodi drugega člena pa so (Poli 2005):  X  Y  Z       S  S  S  S  S  S 1  X  Y  Z       0 0 0 0 0 0 X  Y  Z  S S S    t (49) X  Y  Z  1 1 1 X  Y  Z  S S S 2    t . X  Y  Z  2 2 2 72 Enačbe (48) in (49) skupaj s približnimi vrednostmi neznank se uporabijo za izračun vrednosti matrike koeficientov. 5.3.2.3 Matrična oblika enačb popravkov Enačbe popravkov v izravnavi geometričnega modela zapišemo v matrični obliki z enačbo (22). Enačbe popravkov v grobem ločimo na dve skupini: prva obravnava oslonilne točke (OT), druga pa enačbe zveznosti tirnice (ZT). Matrično obliko enačb (22) lahko zapišemo kot: v  f  B  OT OT OT (50) v  f  B . ZT ZT ZT Matriko koeficientov B se polni glede na število točk, število neznank in število segmentov. Matrika je sestavljena iz zgornjega dela, ki je vezan na oslonilne točke in spodnjega, ki se nanaša na enačbe zveznosti tirnice. Če imamo n točk, je zgornji del matrike B (označimo ga z B1) dimenzije [2  n, 12  št. segmentov]. Vsaka oslonilna točka tvori 2 enačbi (zavzame 2 vrstici), ker pa imamo 2 segmenta, je število stolpcev 24. Vsaka točka zapolni le 12 stolpcev obeh vrstic, ki pripadajo prvemu ali drugemu segmentu, ostalo so ničle. Spodnji del matrike B je razdeljen na tri dele oz. podmatrike. Prva podmatrika ZT0 upošteva samo enakost vrednosti funkcij segmentov, ima dimenzije [6, 24] in jo zapišemo kot: 1 1 1 1 0 0   1 1 1    1 0 0  . (51) ZT   1 1 1 1 0 0  0  1 1   1 1     1   1 73 Druga in tretja podmatrika ( ZT1 in ZT2) upoštevata enakost prvih in drugih odvodov funkcij zunanjih parametrov. Obe imata dimenzije [3, 24] in jih zapišemo: 0 1 2 0  1 0  ZT   0 1 2 0   (52) 1  1 0   0 1 2 0  1 0  0 0 1 0 0  1  ZT   0 0 1 0 0   . (53) 2  1   0 0 1 0 0  1  Celotna matrika B ima torej dimenzije [2  n + 12, 24] in jo zapišemo:  B 1    ZT B   0  . (54)  ZT 1   ZT 2  Matrika oz. vektor neznank Δ ima velikost [24, 1], vektor popravkov v in vektor opazovanj f pa [2  n + 12]. V matrični obliki jih zapišemo:   v  fOT x 1  OT x 1  X 10         v f X 11   OT y 1   OT y 1  Δ   ...  , v   ...  , f   ...  . (55)        v f 20   ZT 2 y   ZT 2 y          v f 20   ZT 2 z   ZT 2 z  5.3.3 Stohastični model Stohastični model je del matematičnega modela, ki se nanaša na natančnost opazovanj. Vsa opravljena opazovanja se statistično spreminjajo in imajo lahko različne natančnosti. V geometričnem modelu se natančnosti opazovanj upoštevajo v obliki simetrične matrike uteži P, ki jo lahko zapišemo kot: 2 1  P     , (56) 0 kjer je Σ kovariančna matrika. 74 Matriko uteži, kjer ima vsaka od n točk in parametri zveznosti tirnice različno natančnost, zapišemo z: 1  2   x 1  2   y 1    ...  2     xn  2  2 P     yn  . (57) 0 2    ZT 1  2   ZT 2    ...  2    ZT 11   2   ZT 12 5.3.4 Reševanje sistema linearnih enačb Po nastavitvi sistema matrik in preureditvi v normalne enačbe oz. v linearen sistem oblike: N  t , (58) lahko izračunamo iskane količine (neznanke). Čeprav se zdi izračun trivialen, je numerično reševanje tako velikega sistema (zlasti, če imamo veliko število oslonilnih točk) zelo kompleksno in lahko vodi v napačne rezultate. Reševanje je oteženo zaradi velikosti in sestave matrike koeficientov normalnih enačb N. Težave se pojavijo pri množenju in inverziji matrike, ki je sestavljena iz veliko ničelnih elementov in elementov, ki zavzemajo zelo velike ali zelo majhne vrednosti. Za numerično reševanje takih sistemov po metodi najmanjših kvadratov se je uveljavilo več metod. Največkrat se uporabljajo metode, ki temeljijo na dekompoziciji LU (angl. lower upper; delitev matrike na spodnjetrikotno (L) in zgornjetrikotno (U) matriko), dekompoziciji QR (delitev matrike na ortogonalno (Q) in zgornjetrikotno (R) matriko), dekompoziciji Cholesky in dekompoziciji singularne vrednosti (angl. singular value). Vse omenjene metode delujejo na dekompoziciji matrike N na produkt več podmatrik z različnimi lastnostmi. Dekompozicija omogoči enostavnejše in robustnejše reševanje zapletenih sistemov enačb. Poleg omenjenih metod se v praksi 75 veliko uporablja tudi iterativna bikonjugirana gradientna metoda (angl. iterative biconjugate gradient method, IBGM), ki se je v našem primeru izkazala za najbolj primerno. Iterativna bikonjugirana gradientna metoda (Press in dr. 2007) je zelo uporabna pri delu s t.i. razpršenimi matrikami (angl. sparse matrix), ki vsebujejo veliko ničelnih elementov. Metoda je splošna verzija osnovne konjugirane gradientne metode, ki sloni na minimiziranju funkcije: 1 f ()   N   t  . (59) 2 Funkcija (44) je najmanjša, ko je njen gradient: f   N  t (60) enak 0. Minimizacija poteka z ustvarjanjem zaporedja iskalnih smeri pk in minimumov xk. V vsaki iteraciji se izračuna količino αk, ki minimizira funkcijo f( xk+αkpk), in vrednost izraza xk+αkpk se priredi naslednji spremenljivki xk+1. Po m iteracijah se dobi minimum za celotni vektorski prostor, ki poda tudi končno rešitev sistema. Osnovna metoda deluje samo, če je matrika N simetrična in pozitivno definitna. V vseh drugih primerih se uporablja bikonjugirano metodo, ki ima še dodatne korake in pogoje. Med drugimi mora zaporedje vektorjev pk zadovoljiti pogoj bikonjukcije in biortogonalnosti. Po več iteracijah metoda reševanja sistema linearnih enačb pripelje do rešitve (vrednosti neznank). Prva rešitev je v posredni izravnavi le prvi približek končne rešitve koraka izravnave. Po nadaljnjih iteracijah pa dobimo rešitev prve izravnave. Rezultati so lahko obremenjeni z grobimi napakami, ki so prisotne v opazovanjih, zato je potrebno izravnave ponavljati in napake iterativno izločati. 5.3.5 Izločanje grobih napak Zaradi metode določanja so v samodejno določenih koordinatah oslonilnih točk večkrat prisotne napake. Če je prišlo do slabe korelacije med referenčnimi podatki o cestah in cestami, ki so bile dobljene iz satelitskih 76 posnetkov, imajo lahko take točke grobe napake. Ker pri samodejni določitvi točk vizualna kontrola ustreznosti le-teh ni možna, moramo slabe točke izločiti samodejno v postopku izravnave. Poleg izločitve točk z grobimi napakami je za visoko natančnost rezultatov potrebna tudi določitev najboljših točk in njihova ustrezna uporaba na tak način, da imajo največji možni vpliv na končni rezultat. Zaradi takega večplastnega preizkušanja točk je geometrični model razdeljen na dva zaporedna dela: algoritem RANSAC in robustna ocena. 5.3.5.1 Algoritem RANSAC Prvi del je namenjen izločanju oslonilnih točk z grobimi napakami. Za izločanje smo uporabili algoritem RANSAC (RANdom SAmple Consensus), ki je že bil uporabljen za izločanje točk, na primer v programu BINGO (Kruck 2001), uporabili pa so ga tudi pri določitvi parametrov kamere (Zhou in dr. 2013). RANSAC (Fischler in Bolles 1981) je iterativna metoda za oceno parametrov matematičnega modela iz opazovanj, ki lahko vsebujejo tudi grobe napake. Algoritem določi vrednost parametrov z verjetnostjo, ki se veča z naraščanjem števila ponovitev iteracij. RANSAC je zanesljiv tudi pri velikem številu grobih napak in velikih napakah, kjer nekatere druge metode za izločanje grobih napak (npr. data snooping) niso uporabne. Čeprav je RANSAC zelo robustna metoda, je lahko neučinkovita, če je grobih napak več kot 50 %. Pri velikem številu grobih napak potrebujemo veliko iteracij, kar lahko vodi v počasno izvajanje. Nasploh pa je postopek težko avtomatizirati zaradi obvezne nastavitve treh začetnih parametrov: število iteracij algoritma, vrednost odstopanja opazovanja od modela in število opazovanj, ki se prilegajo modelu. Začetne vrednosti teh parametrov so odvisne predvsem od kakovosti in razporeditve točk ter uporabljenega senzorja. Zaradi tega so bile te vrednosti večinoma določene empirično in z iterativnimi postopki v modelu. V geometričnem modelu se RANSAC izvaja v več korakih. Najprej se izvede izravnava z najmanjšim možnim številom oslonilnih točk. Za določitev 24 parametrov zunanje orientacije z uporabljenim modelom vrstičnega senzorja je potrebnih najmanj šest točk. Izmed vseh samodejno določenih točk se te 77 točke določijo naključno. Ker se morajo izračunani parametri prilagajati celotnemu posnetku, je tudi določitev položaja osnovnih šestih točk pomembna. Zaradi tega izbira točk ni povsem naključna, ampak je podvržena posebnim pravilom. V tem koraku se iz vseh točk ustvari 100 šesteric točk, ki se jim izračunajo naslednje lastnosti: - standardni odklon koordinat točk, - seštevek razdalj med vsemi točkami, - odstopanje točk od premice skozi točke in - število točk v sredini posnetka. Odstopanje točk od premice smo računali s statistiko  2 (Harris Geospatial 2017): 2   a, b N   y  a  b x 2 , (61) i i  i 1 kjer je N število opazovanj, a in b sta koeficienta premice, x in y pa koordinati točke. Prve tri lastnosti se normalizirajo s povprečjem in standardnim odklonom rezultatov. Med vsemi skupinami šestih točk se izbere tisto, ki ima največji standardni odklon med točkami, največji seštevek razdalj med točkami, največje odstopanje  2 od premice in ima vsaj eno točko v sredini posnetka (kvadratno območje, ki zavzema 16 % posnetka). Poskusi so pokazali, da je povprečje med rezultati prvih treh lastnosti dovolj dober pokazatelj razporeditve in zato ne potrebujemo utežnih funkcij. Zadnja lastnost je pomembna, ker se le s prvimi tremi lastnostmi običajno izberejo točke na robovih posnetka. Če točk v sredini ni, algoritem vzame najboljši rezultat prvih treh lastnosti. Postopek je torej iterativen, vendar vseeno zelo hiter in nujno potreben pri iskanju pravilne rešitve modela, posebno, če je prisotnih več grobih napak. Prva iteracija algoritma RANSAC vzame najbolje razporejeno šesterico točk in prične z izravnavo. Izravnava poteka po modelu, ki je bil opisan v prejšnjih poglavjih in iz katerega dobimo parametre zunanje orientacije. Izračunani parametri se nato vstavijo v kolinearne enačbe, s katerimi preverimo prileganje vseh samodejno določenih točk modelu. Točke 78 preverimo z vstavljanjem slikovnih koordinat in računanjem prostorskih koordinat, ki se nato primerjajo s koordinatami samodejno pridobljenih točk. Če je razlika koordinat znotraj določenih meja, ki so običajno za vsak senzor različne, si algoritem točko zapomni. Ko se kontrola izvede za vse točke, se preveri odstotek vseh točk, ki so znotraj predhodno določene meje. V kolikor je teh točk enako ali več kot je bilo določeno v začetnih parametrih, se ponovno izvede izravnava samo s temi točkami. Po izravnavi se izračuna natančnost v obliki standardnega odklona, ki si ga algoritem zapomni. Izračunan standardni odklon bo merilo natančnosti izravnave in referenca za vse nadaljnje iteracije, saj bo iteracija z najmanjšim standardnim odklonom oziroma njene točke ohranjene za drugi del računanja modela. Zaradi specifičnosti senzorjev in spremenljive točnosti samodejno določenih točk je vse omenjene meje težko fiksno določiti. Odstotek dobrih točk smo določili empirično, medtem ko se meja, ki določa ustreznost točk, določa iterativno. Če po prvem izvajanju vseh iteracij število točk ne doseže ustreznega odstotka, se meja poveča in celotna obdelava ponovno izvede. Postopek se lahko tudi večkrat ponovi, kar pa da običajno slabše rezultate in kaže na slabo samodejno določitev točk. 5.3.5.2 Robustna ocena Drugi del modela uporablja samo točke, ki jih ni izločil algoritem RANSAC. Večina preostalih točk nima večjih grobih napak, ampak praviloma le manjše napake, ki jih je zelo težko izslediti. Čeprav so v večini primerov rezultati prvega dela natančni pod velikostjo dveh pikslov, se jih lahko dodatno izboljša z robustno oceno (angl. Robust estimation). Postopek robustne ocene so v fotogrametriji začeli uporabljati že v osemdesetih letih (Kubik in dr. 1987), sama ideja pa se je izoblikovala v šestdesetih letih (Huber 1964). Robustne cenilke so relativno neobčutljive na variacije porazdelitvene funkcije opazovanj in posledično na grobe in slučajne napake. V nasprotju z metodo najmanjših kvadratov robustna ocena ne minimizira kvadratov napak, ampak izbrano funkcijo napak. V fotogrametriji se izbrana funkcija uporablja za izračun uteži opazovanj glede na njihovo točnost, ki jo dobimo 79 po vsaki izravnavi. V praksi to pomeni, da se opazovanjem z večjimi popravki iterativno zmanjšuje utež in posledično vpliv na rezultate izravnave. V prvem delu (RANSAC) je bila matrika uteži opazovanj enotska, v drugem delu pa se sestavi glede na rezultate kriterijske funkcije (uspešnost korelacije točk) iz modula za samodejno določanje točk. Matrika predstavlja začetno stanje, ki se spreminja z vsako ponovitvijo izravnave. V primeru majhnih popravkov bodo končni parametri znani že po prvi izravnavi. Število iteracij je poleg natančnosti točk odvisno tudi od izbrane funkcije napak in meje, ki določa, kdaj je popravek dovolj velik za zmanjšanje njegove uteži v izravnavi. Izbira prave funkcije je odvisna od več faktorjev, med katerimi je najpomembnejša točnost oslonilnih točk. Ker se slednja spreminja, smo funkcijo znova izbrali z empiričnimi testi. Vrsto utežne funkcije in njeno obliko smo izbirali med tremi različicami. Preizkusili smo Dansko metodo (Kubik in dr. 1987), metodo po Kleinu (Klein in Förstner 1984) in metodo z osnovno hiperbolično funkcijo (poglavje 5.6). Danska metoda uporablja eksponentno funkcijo, nove uteži pa smo računali z (Kubik in dr. 1987): 2  v  P  exp i (62) i    ,  2 2  kjer je Pi nova utež, v popravek opazovanja in σ standardni odklon. Metoda po Kleinu deluje s hiperbolično funkcijo (Klein in Förstner 1984): 1 P  P (63) i j 1  a  v i  d i v kateri je: Pj a  i 4 4 .  r  82 d  5 . 3  (64) 4 81 Q  Q  ,  a priori 80 Pj je stara utež in r je število nadštevilnih opazovanj. Tudi tretja preizkušena metoda deluje s hiperbolično funkcijo: 1 P  . (65) i v 1 i  Mejo, ki določa, kdaj ima opazovanje preveliko napako, smo določili pri dveh standardnih odklonih (2 σ). Algoritem najprej preveri katera opazovanja imajo popravke večje od določene meje in nato le-tem spremeni utež glede na izbrano funkcijo. Če je popravek tudi v drugi iteraciji prevelik, se opazovanju ponovno izračuna utež. Postopek se ponavlja, dokler niso popravki vseh opazovanj znotraj meje. Na ta način se točke ne odstrani, ampak se jim le zmanjša vpliv pri izravnavi in izračunu končnih parametrov (neznank). S tem se ohrani razporeditev točk oz. večja pokritost posnetka s točkami, kar je zelo pomembno za doseganje visokih točnosti z rigoroznimi geometričnimi modeli. 5.4 Geometrični model polnoslikovnega senzorja Za ortorektifikacijo polnoslikovnih senzorjev je bil uporabljen geometrični model, ki ne vsebuje časovne komponente. Zaradi tega je model enostavnejši od vrstičnega in računa le šest parametrov zunanje orientacije, kar je enako kot pri letalski fotogrametriji. Tudi geometrični model polnoslikovnega senzorja uporablja kolinearne enačbe, ki pa imajo dodane še parametre notranje orientacije. 5.4.1 Funkcionalni model Funkcionalni model polnoslikovnega senzorja temelji na osnovnih kolinearnih enačbah. Razlika z vrstičnim modelom je predvsem v številu in vrsti neznank. V primeru polnoslikovnega senzorja se celoten posnetek zajame v trenutku, kar pomeni, da ima posnetek le en perspektivni center in en niz rotacij. Skupno je torej le šest parametrov zunanje orientacije – koordinate perspektivnega centra ( XS, YS, ZS) in rotacije ( ωS, φS, κS). 81 Ker imajo majhni sateliti omejen prostor za komponente, večinoma uporabljajo polnoslikovni senzor, ki zaradi enostavnosti delovanja zasede malo prostora. Ker imajo tovrstni sateliti manjši uporabni volumen, lahko pride do večjih temperaturnih razlik v ohišju in tudi slabšega odvajanja visokih temperatur. Zaradi termalne nestabilnosti se materiali krčijo in raztezajo, kar vodi tudi v spreminjanje oblike in položaja komponent. Te spremembe so sicer zelo majhne, vendar lahko vplivajo na opazovanja in kasneje na delovanje modela. Zaradi tega v modelu polnoslikovnega senzorja tudi parametre notranje orientacije obravnavamo kot neznanke. Za modeliranje notranje orientacije smo izbrali tri parametre – goriščno razdaljo f, in premik detektorja v slikovni ravnini, ki ga definirata količini x0 in y0. Dopolnjena oblika kolinearnih enačb je (Mikhail in dr. 2001): m  X  X  m  Y  Y  m  Z  Z 11  S  12  S  13  S  x  x   f 0 m  X  X  m  Y  Y  m  Z  Z 31  S  32  S  33  S  (66) m  X  X  m  Y  Y  m  Z  Z 21  S  22  S  23  S  y  y   f 0 m  X  X  m  Y  Y  m  Z  Z 31  S  32  S   S  , 33 pri čemer so parametri enaki kot pri enačbi (29). V okviru raziskave je bil izdelan splošen in rigorozen geometrični model za polnoslikovne senzorje. S pravilno predpripravo metapodatkov lahko model deluje praktično z vsemi snemalnimi sistemi, ki uporabljajo polnoslikovne senzorje. Poleg tega model rigorozno opiše razmerje med posnetkom in objektnim prostorom s kolinearnimi enačbami. Geometrični model polnoslikovnega senzorja je podrobno opisan v naslednjih podpoglavjih. 5.4.1.1 Linearizacija Linearizacija kolinearnih enačb poteka enako kot v primeru vrstičnega senzorja, le da imamo pri polnoslikovnem senzorju več vrst neznank (parametre notranje in zunanje orientacije). Tokrat se parcialni odvodi računajo za devet neznank ( XS, YS, ZS, ωS, φS, κS, f, x0, y0). Parcialni odvodi prvega člena glede na neznanke za x in y so (Kraus 2007): 82 x  f   ( m Nx  m D)  b 2 13 11 1 , 1 X  D S x  f   ( m Nx  m D)  b 2 23 21 , 1 2 Y  D S x  f   ( m Nx  m D)  b 2 33 31 , 1 3 Z  D S x  f  Nx    (( Y  Y ) m ( Z  Z ) m )  ( Y  Y ) m  ( Z  Z ) m   b S 33 S 23 S 31 S 21 , 1 4   D  D  S x  f  Nx   ( Nx cos  Ny sin)  D cos  b ,15   D  D  S x  f   Ny  b ,16   D S x  1 b ,17 x  0 x  0  b ,18 y  0 x  Nx    b ,19 f  D y  f   ( m Ny  m D)  b (67) 2 13 12 2 1 , X  D S y  f   ( m Ny  m D)  b 2 23 22 2,2 Y  D S y  f   ( m Ny  m D)  b 2 33 32 2,3 Z  D S y  f  Ny    (( Y  Y ) m  ( Z  Z ) m )  ( Y  Y ) m  ( Z  Z ) m   b S 33 S 23 S 32 S 22 2,4   D  D  S 83 y  f  Ny   ( Nx cos  Ny sin)  D sin  b 2,5   D  D  S y  f  Nx  b 2,6   D S y   0  b 2,7 x  0 y  1 b 2,8 y  0 y  Ny    b . 2,9 f  D Enačbe (67) skupaj s približnimi vrednostmi neznank se uporabijo za izračun vrednosti matrike koeficientov. 5.4.1.2 Matrična oblika enačb popravkov Enačbe popravkov v posredni izravnavi geometričnega modela zapišemo v matrični obliki z enačbo (22). Vse enačbe popravkov so vezane na uporabljene oslonilne točke, število členov enačb pa je vezano na število neznank. Ker vsaka oslonilna točka tvori 2 enačbi, imamo pri n točkah matriko koeficientov B velikosti [2  n, 9]. Matrika B ima obliko:   x  x  x  x 1 1 1 1   ...   X  Y  y  f S S 0    y  y  y  y 1 1 1 1   ...  X  Y  y  f   S S 0  B   ... ... ... ...  . (68)   x  x  x  x n n n n   ...  X  Y  y  f   S S 0   y  y  y   y n n n n ...   X  Y  y  f S S 0  Matrika oz. vektor neznank Δ ima velikost [9, 1], vektor popravkov v in vektor opazovanj f pa [2  n]. V matrični obliki jih zapišemo: 84    v  fOT x 1  OT x 1  X          v f Y   OT y 1   OT y 1  Δ   ...  , v   ...  , f   ...  . (69)        v f y 0   OTnx   OTnx        v   f f   OTny   OTny  5.4.2 Stohastični model Matrika uteži za geometrični model polnoslikovnega modela je simetrična in vsebuje samo natančnosti oslonilnih točk. Matriko uteži, kjer ima vsaka od n točk različno natančnost, zapišemo z: 1  2   x 1   2   y 1  2 P     ...  . (70) 0   2   xn   2    yn 5.4.3 Reševanje sistema linearnih enačb in izločanje grobih napak Reševanje sistema linearnih enačb in izločanje grobih napak poteka zelo podobno kot pri geometričnem modelu za vrstične senzorje. Največja razlika je v vrednosti najmanjšega možnega števila oslonilnih točk, s katerim se izvede izravnava, saj jih polnoslikovni model za določitev devetih parametrov zunanje in notranje orientacije potrebuje le pet. 5.5 Implementacija Geometrični modeli senzorjev so bili implementirani v samostojen procesni modul za izračun parametrov geometričnega modela, ki se izvede za modulom za samodejno določanje oslonilnih točk. 85 5.5.1 Modul za izračun parametrov geometričnega modela Modul za izračun parametrov geometričnega modela je najbolj zapleten del postopka in osrednji del knjige. Poleg izhodne podatkovne strukture modula za pripravo podatkov geometrični model uporablja še tekstovno datoteko s samodejno določenimi oslonilnimi točkami in DMR. V modelu se uporablja veliko iterativnih korakov in napredno določanje ter izločanje oslonilnih točk z grobimi napakami. Prav zaradi iterativnosti in iskanja konvergenc lahko postopek traja različno dolgo časa, odvisno od posnetka in števila dobljenih oslonilnih točk. V modulu se izvaja vrstični ali polnoslikovni model, odvisno od vrste posnetkov. V nekaterih segmentih modula se geometrična modela izvajata enako (npr. branje točk), v drugih (npr. izravnava) pa je postopek drugačen zaradi različnega načina snemanja in števila neznank. Polnoslikovni senzor je preprostejši in računsko manj zahteven (razlike so opisane v začetku poglavja 5). V modulu za izračun parametrov geometričnega modela se najprej prebereta DMR in datoteka s točkami. Iz datoteke se preberejo vsi ključni podatki za model: identifikacijska številka (ID) točke, X in Y koordinati točke iz referenčnih podatkov, x in y slikovni koordinati iz surovega posnetka in natančnost določitve točke (kriterijska funkcija oz. stopnja korelacije). Iz DMR- ja se za vsako točko pridobi še njena višina. Slikovne koordinate se preračunajo iz slikovnega koordinatnega sistema v koordinatni sistem kamere z:  N  v x   v   p  2 v  (71)  N  s y   s   p .  2 v  kjer je pv velikost piksla, Nv število pikslov v vrstici, Ns število pikslov v stolpcu, v in s pa sta slikovni koordinati točke. Pri vrstičnem skenerju nastavimo koordinato x na vrednost 0, saj ima vsaka vrstica svoje izhodišče. Nato se v modulu slikovne koordinate popravijo zaradi vpliva ukrivljenosti Zemlje in atmosferske refrakcije. Glede na višino leta in kot gledanja se koordinatam spremeni vrednosti za omenjena vpliva. Tako pripravljene 86 koordinate vstopijo v geometrični model. Podrobno delovanje geometričnega modela prikazuje slika 11. Iz nje je razvidno da, je celotna obdelava razdeljena na dva dela, ki potekata drug za drugim. Oba dela delujeta iterativno na več ravneh. Poleg tega se v trenutni implementaciji celoten postopek ponovi petkrat, na koncu pa obvelja najboljši rezultat. Zaradi slabo razporejenih in netočnih oslonilnih točk lahko postane izravnava nestabilna. V takih primerih so rezultati po robustni oceni preslabi za uporabo in kot končni rezultati obveljajo le parametri, dobljeni po izravnavi z algoritmom RANSAC, ki je na take razmere manj občutljiva. Rezultat izravnave so končni parametri zunanje in notranje orientacije, ki se zapišejo v podatkovno strukturo. Poleg parametrov se v podatkovno strukturo zabeležijo tudi osnovni podatki o posnetku, ki se uporabijo še v modulu za izvedbo ortorektifikacije (glej poglavje 6). Na koncu modula se preveri še notranja natančnost modela. Parametre se uporabi za izračun prostorskih koordinat oslonilnih točk iz pikselskih koordinat s kolinearnimi enačbami. Tako izračunane prostorske koordinate se odštejejo od koordinat, pridobljenih v modulu za samodejno določanje oslonilnih točk. Iz razlik se izračuna RMSE za obe koordinati in še skupni po enačbi: n r  2 i RMSE i   1 , (72) n kjer je n število opazovanj in r razlika med koordinatama. RMSE oslonilnih točk sicer ne poda absolutne natančnosti, vendar je vseeno dober pokazatelj prileganja modela vhodnim točkam. Ob odsotnosti kontrolnih točk je to edina zanesljiva mera natančnosti modela. 87 88 …nadaljevanje. Slika 11: Shematizacija poteka izračunov parametrov geometričnega modela. 5.6 Testiranje Uspešnost in robustnost delovanja modula za izračun parametrov geometričnega modela smo preverili z uporabo treh empiričnih poskusov za vsako vrsto podatkov. Poleg tega smo s poskusi izbrali tudi metode, ki bodo uporabljene v dejanski izvedbi postopka, najprimernejše vhodne podatke ter vrednosti (začetnih) parametrov. 5.6.1 Vpliv števila oslonilnih točk na izračun parametrov geometričnega modela Določanje točk poteka po izsekih posnetka. Na vsakem izseku smo določili več oslonilnih točk. Število izsekov in oslonilnih točk se spreminja glede na vrsto posnetka. Posnetki RapidEye, na primer, zajemajo večjo površino in so zaradi tega razdeljeni na več izsekov kot posnetki WorldView-2. Število izsekov in točk smo določili z empiričnimi poskusi na več posnetkih. Zaradi pomanjkanja cest na nekaterih izsekih, različne pokrovnosti in napačne detekcije, je število točk za vsak posnetek različno. S poskusi smo določili tudi 89 najvišje število točk na izsek, ki še znatno vpliva na natančnost delovanja postopka. Število teh točk je največje število (niz) točk za poskuse. Za vsak posnetek smo poleg največjega niza točk ustvarili še več nizov z izločanjem točk iz največjega niza glede na njihovo kakovost in razporeditev. Vsak niz je bil uporabljen v geometričnem modelu za izračun neznanih parametrov. Natančnost smo ocenili z računanjem napake RMSE, ki jo dobimo s primerjavo med ročno izmerjenimi koordinatami kontrolnih točk in koordinatami, ki jih dobimo z rezultati geometričnega modela. 5.6.1.1 Poskusi s posnetki RapidEye Za vsak posnetek smo izdelali več nizov oslonilnih točk s 15 do 500 točkami. Vsak niz smo uporabili v geometričnem modelu, ki je deloval s 100 iteracijami algoritma RANSAC in Kleinovo metodo za robustno oceno (glej str. 80), ki je pri opravljenih poskusih delovala zelo robustno in je dosegala najboljše rezultate. RANSAC smo ponavljali, dokler v novonastalem nizu ni ostalo vsaj 90 % točk prvotnega niza, torej dokler ni bilo izločenih manj kot 10 % prvotnih točk. Prag je bil določen glede na rezultate ocene točnosti določevanja točk (glej podpoglavje 4.3.1). Če je bilo v nizu več grobih napak, je njihov vpliv na rezultate zmanjšala metoda robustne ocene. Rezultati za vse posnetke so prikazani v preglednici 10. Iz preglednice je razvidno, da RANSAC običajno izloči okoli 10 % točk ne glede na število začetnih točk. Neizločene točke naj ne bi imele grobih napak oz. so te majhne. Preostale točke se uporabijo v izravnavi z robustno oceno. Končna točnost rezultatov, ki jo podaja RMSE na ročno merjenih kontrolnih točkah, je skoraj v vseh primerih manj kot piksel, ne glede na število uporabljenih oslonilnih točk (slika 12). Napaka RMSE je v primeru manjšega števila točk večja, vendar je razlika majhna. Rezultati kažejo na uspešno delovanje modela tudi z majhnim številom oslonilnih točk. Kljub temu v večini primerov večje število točk (100 in več) zagotavlja boljše rezultate, zato se v geometričnem modelu priporoča uporaba večjega števila oslonilnih točk. Ker je določanje točk samodejno, je določitev več sto točk hitro in enostavno ter pri izračunu parametrov občutno večje število točk geometričnega modela ne upočasni. 90 Preglednica 10: Točnost rezultatov izravnave, ki so bile dobljene z različnim številom oslonilnih in kontrolnih točk za posnetke RapidEye. Število Število Število RMSE OT RMSE KT Posnetek uporabljenih oslonilnih točk kontrolnih točk [piksel] [piksel] oslonilnih točk 501 460 22 0,59 0,82 303 273 22 0,62 0,85 207 196 22 0,58 0,83 Radgona 96 87 22 0,65 0,87 50 46 22 0,59 0,84 32 32 22 0,82 0,87 16 15 22 0,77 0,96 467 433 26 0,67 0,95 295 276 26 0,70 0,86 202 187 26 0,73 0,90 Koper 98 92 26 0,71 0,96 56 53 26 0,72 0,94 32 31 26 0,69 1,01 15 14 26 0,88 1,26 325 297 27 1,01 0,96 203 184 27 1,01 0,99 101 91 27 0,95 0,91 Bohinj 53 48 27 0,94 1,17 29 27 27 0,85 1,05 16 15 27 0,80 0,97 Slika 12: RMSE na kontrolnih točkah za vsak točkovni niz za posnetke RapidEye. 91 5.6.1.2 Poskusi s posnetki WorldView-2 Za vsak posnetek smo izdelali več nizov oslonilnih točk, ki so vsebovali od približno 15 do približno 200 točk. Vsak niz smo uporabili v geometričnem modelu, ki je zaradi nekoliko slabše točnosti točk deloval s 300 iteracijami algoritma RANSAC in Kleinovo metodo za robustno oceno. RANSAC smo ponavljali, dokler v novonastalem nizu ni ostalo vsaj 90 % točk prvotnega niza. Rezultati za vse posnetke so prikazani v preglednici 11, iz katere je razvidno, da tudi v teh primerih RANSAC večkrat izloči okoli 10 % točk ne glede na število začetnih točk. Končna točnost rezultatov, ki jo podaja RMSE na ročno merjenih kontrolnih točkah, je v vseh primerih nad velikostjo enega piksla (slika 13), na posnetku Ljubljana doseže celo vrednost nad velikostjo petih pikslov. Najboljše rezultate smo dobili s posnetkom Jesenice, ki pa je zelo specifičen glede razporeditve točk. Pri tem posnetku so tako oslonilne kot kontrolne točke prisotne samo v spodnjem delu posnetka, ki leži znotraj Slovenije (priloga A). Zaradi tega je točnost severnega dela posnetka neznana. Vrednosti RMSE na kontrolnih točkah so v rangu vrednosti RMSE glede na racionalne polinomske koeficiente RPC (preglednica 8), kar kaže na dobro delovanje geometričnega modela. Na dober model kažejo tudi RMSE oslonilnih točk. Do podobnih zaključkov lahko pridemo tudi pri posnetku Bled. Model se je dobro prilagodil točkam (majhen RMSE oslonilnih točk), večje napake na kontrolnih točkah pa so nastale zaradi slabe razporeditve oslonilnih točk po posnetku (slika 18). Zaradi pomanjkanja točk v hribovju (predvsem južni in zahodni del posnetka) imajo kontrolne točke na tem delu večja odstopanja. Najslabše in nepričakovane rezultate smo dobili pri posnetku Ljubljana. Kljub temu, da je bila razporejenost in natančnost točk dobra, je model dosegal zelo slabo notranjo natančnost (RMSE oslonilnih točk), kar je imelo za posledico tudi slabe rezultate na kontrolnih točkah. Pravega vzroka za slabše delovanje modela v tem primeru ne poznamo, sklepamo pa lahko, da je posledica načina snemanja (asinhrono, usmerjeno proti severu ipd.) in obdelave posnetka po snemanju, ki ju geometrični model trenutno še ne upošteva. Zaradi manjšega števila razpoložljivih posnetkov za preizkušanje je trenutno težko bolj natančno določiti vzroke slabših rezultatov. 92 Preglednica 11: Točnost rezultatov izravnave, ki so bile dobljene z različnim številom oslonilnih in kontrolnih točk za posnetke WorldView-2. Število Število Število RMSE OT RMSE KT Posnetek uporabljenih oslonilnih točk kontrolnih točk [piksel] [piksel] oslonilnih točk 75 70 24 0,92 1,54 50 48 24 1,17 1,44 Jesenice 31 28 24 0,98 1,58 15 14 24 0,74 1,78 158 149 30 0,99 2,30 102 97 30 1,11 3,17 Bled 50 47 30 1,15 3,70 30 29 30 1,35 3,80 15 14 30 0,88 4,88 206 191 30 2,61 3,94 152 145 30 2,85 3,94 101 92 30 3,05 5,20 Ljubljana 48 45 30 3,06 4,10 31 29 30 3,33 5,49 16 15 30 2,97 5,77 Slika 13: RMSE na kontrolnih točkah za vsak točkovni niz posnetkov WorldView-2. Če upoštevamo število oslonilnih točk, ugotovimo, da je RMSE večja v primeru manjšega števila točk in v vseh primerih narašča do najmanjšega niza, kar je v skladu s pričakovanji. Model sicer uspešno deluje tudi z majhnim 93 številom oslonilnih točk, vendar so v nekaterih primerih točnosti z manjšimi nizi točk precej slabše. Tudi zaradi slabše določitve točk je zelo priporočljiva uporaba večjega števila oslonilnih točk. Kot kažejo rezultati, je zelo pomembna enakomerna razporeditev točk, ki pa jo na nekaterih območjih, kjer ni cest oz. jih je malo ali so slabo vidne, težko dosežemo. 5.6.2 Metode za robustno oceno Za postopek je pomembna tudi izbrana metoda za robustno oceno. Najboljšo metodo smo izbrali s preizkušanjem in primerjavo treh različnih metod: Danske metode, metode po Kleinu in metode s hiperbolično funkcijo. Metoda, ki je dosegla najboljše rezultate pri uporabljenih posnetkih, je bila tudi implementirana v ortorektifikacijski postopek. 5.6.2.1 Poskusi s posnetki RapidEye Primerjavo metod za robustno oceno smo izvedli z nizom oslonilnih točk, ki je dosegel najboljše rezultate v prejšnjem poskusu. Uporabili smo 100 iteracij algoritma RANSAC; kot že omenjeno, smo jih ponavljali, dokler v novonastalem nizu ni ostalo vsaj 90 % točk prvotnega niza. V preglednici 12 so prikazane napake RMSE na oslonilnih in kontrolnih točkah po izravnavi z algoritmom RANSAC in po izravnavi z metodami robustne ocene, ki naj bi omejile vpliv preostalih oslonilnih točk z napakami. Iz rezultatov je razvidno, da robustna ocena le redko izboljša rezultate izravnave z algoritmom RANSAC. Rezultati večinoma ostanejo nespremenjeni, kar kaže na dobro delovanje iterativne metode RANSAC in na manjše število grobih napak, ki ne presegajo 10 % točk prvotnega niza. Robustna ocena bi prišla do izraza v primeru večjega števila grobih napak. Pri primerjanju metod robustne ocene (slika 14) najbolj izstopa slabše delovanje Danske metode, ki uporablja eksponentno utežno funkcijo. Metoda iz neznanih razlogov v nekaterih primerih ni zanesljiva, saj privede do velikih napak. Delovanje ostalih metod, ki uporabljajo hiperbolično funkcijo, je uspešno in zanesljivo. 94 Izbira najboljše metode za robustno oceno ni trivialna, saj lahko nekatere metode celo poslabšajo rezultate. Metode, ki uporabljajo hiperbolično funkcijo, so se v primeru posnetkov RapidEye izkazale za bolj zanesljive in robustne. Kljub temu so se rezultati izboljšali le v nekaterih primerih. Med mnogimi ponavljanji poskusov smo ugotovili, da je uporaba robustne ocene smiselna v primerih, ko je po izravnavi z algoritmom RANSAC med točkami prisotnih vsaj nekaj grobih napak. Zaradi tega je priporočljivo odstraniti največ 10 % oslonilnih točk, tudi če je med opazovanji prisotnih še več grobih napak. Preglednica 12: Primerjava metod za robustno oceno na posnetkih RapidEye. RANSAC Robustna ocena Posnetek RMSE OT RMSE KT Metoda RMSE OT RMSE KT [piksel] [piksel] [piksel] [piksel] 0,76 0,82 Danska 0,85 0,86 Radgona 0,66 0,83 Hiperbolična 0,66 0,83 0,76 0,82 Klein 0,77 0,83 0,81 0,93 Danska 2,64 3,29 Koper 0,75 0,87 Hiperbolična 0,71 0,84 0,75 0,88 Klein 0,74 0,88 1,01 0,98 Danska 0,98 0,89 Bohinj 0,97 0,94 Hiperbolična 0,97 0,94 0,98 0,96 Klein 0,98 0,96 Slika 14: Obnašanje preizkušenih metod robustne ocene za posnetke RapidEye. 95 5.6.2.2 Poskusi s posnetki WorldView-2 Primerjavo metod za robustno oceno smo izvedli z nizom oslonilnih točk, ki je dosegel najboljše rezultate v prejšnjem poskusu. Druge začetne nastavitve so ostale enake kot pri prejšnjem poskusu. V preglednici 13 so prikazane napake RMSE na oslonilnih in kontrolnih točkah po izravnavi z algoritmom RANSAC in na koncu, po izravnavi z metodami robustne ocene. Tudi pri tem poskusu so bili rezultati metod zelo različni po posnetkih (slika 15). Pri posnetkih Jesenice in Bled so vse metode dosegale približno enake rezultate. Nekoliko slabše je delovala le Danska metoda, ki pa je imela tudi najslabše natančnosti točk po izravnavi z algoritmom RANSAC. Iz preglednice je razvidno, da je popravek robustne ocene opazen le v primeru slabših rezultatov izravnave z algoritmom RANSAC, kar vendarle dokazuje primerno delovanje vseh metod. Zaradi slabega delovanje modela pri posnetku Ljubljana je bila v tem primeru ocena in primerjava metod robustne ocene težavna. Metode so imele precejšnje težave pri izračunu parametrov, Danska metoda je npr. le enkrat dosegla smiselne rezultate (v drugih primerih so obveljali rezultati po izravnavi z algoritmom RANSAC). Najbolj robustno je delovala Kleinova metoda. Preglednica 13: Primerjava metod za robustno oceno na posnetkih WorldView-2. RANSAC Robustna ocena Posnetek RMSE OT RMSE KT Metoda RMSE OT RMSE KT [piksel] [piksel] [piksel] [piksel] 1,69 2,30 Danska 0,97 1,62 Jesenice 0,89 1,44 Hiperbolična 0,90 1,47 0,87 1,33 Klein 0,92 1,54 2,01 4,71 Danska 1,03 3,71 Bled 1,45 4,09 Hiperbolična 1,01 2,92 1,02 1,89 Klein 0,99 2,30 2,43 3,11 Danska 2,51 3,21 Ljubljana 1,95 2,26 Hiperbolična 3,08 4,27 2,74 4,57 Klein 2,70 4,08 Tudi v primeru posnetkov WorldView-2 se je najbolje izkazala uporaba hiperbolične funkcije. Čeprav je v večini primerov tudi eksponentna funkcija dobro delovala, je njena zanesljivost v taki obliki vprašljiva. Tudi pri teh primerih se je pokazalo, da je uporaba robustne ocene smiselna le, ko je po 96 izravnavi z algoritmom RANSAC med točkami prisotnih vsaj nekaj grobih napak. Zaradi tega je priporočljivo odstraniti manjše število oslonilnih točk, tudi če je med opazovanji prisotnih še več grobih napak. Slika 15: Lastnosti preizkušenih metod robustne ocene za posnetke WorldView-2. 5.6.3 Detekcija in izločanje grobih napak Robustnost geometričnega modela je bila preizkušena z umetnim dodajanjem grobih napak oslonilnim točkam. Za vsak posnetek smo izbrali večji in manjši niz oslonilnih točk. Pri posnetkih RapidEye je bil večji niz tisti, ki je dosegel najboljše rezultate pri drugem poskusu, medtem ko je pri posnetkih WorldView-2, zaradi manjšega števila določenih točk, večji niz vseboval vse točke. Manjši niz je predstavljal približno četrtino večjega niza ali tretjino v primeru manjšega števila točk v večjem nizu. V večjem nizu smo izbrali deset natančnih in enakomerno razporejenih točk ter obema slikovnima koordinatama točk dodali napake velikosti od enega do deset pikslov. V manjšem nizu smo izbrali le pet točk, ki smo jim dodali napake velikosti od enega do pet pikslov. V kolikor je imel niz le okoli trideset točk, smo izbrali približno 10 % točk, ki smo jim dodali napake ali pa smo zmanjšali mejo za izločanje točk. Oba niza smo uporabili v geometričnem modelu ter analizirali natančnost rezultatov in delovanje samodejnega izločanja točk. 97 5.6.3.1 Poskusi s posnetki RapidEye Pri preizkušanju izločanja grobih napak smo uporabili enake začetne nastavitve in metode kot pri prejšnjih poskusih. Rezultati izravnave so prikazani v preglednicah 14 in 15. Pri posnetku Radgona je RANSAC v večjem nizu izločil 13 oslonilnih točk, med njimi sedem točk z umetno dodanimi napakami od štirih do desetih pikslov. Le tri točke z najmanjšimi dodanimi napakami so ostale v nizu, ki je bil nato uporabljen v izravnavi z robustno oceno. Natančnost končnih rezultatov je primerljiva z natančnostjo rezultatov brez dodanih napak (preglednica 10). Iz rezultatov lahko sklepamo, da robustna ocena v večini primerov popravi natančnost rezultatov izravnave (slika 16). Preglednica 15 prikazuje odstopanja točk z napakami, ki so ostale v nizu uporabljenih točk do konca geometričnega modela. Vrednosti so zelo podobne dodanim napakam, kar pomeni, da jih je geometrični model prepoznal in omejil njihov vpliv na rezultate. Uspešno izločanje je predvsem posledica večjega števila točk in robustne ocene, ki je zmanjšala vpliv točk z napakami. Podobne zaključke lahko postavimo tudi pri manjšem nizu. Preglednica 14: Rezultati izravnave s številom prisotnih in izločenih oslonilnih točk z dodanimi napakami, številom vseh in uporabljenih oslonilnih točk ter RMSE po izravnavi z algoritmom RANSAC in na koncu postopka (robustna ocena) za posnetke RapidEye. RANSAC Robustna ocena [%] ke Posnetek Niz e/ ene ločanja ablj OT KT OT KT isotne napa Vse točk Upor točke Meja iz Izločene napake/ Pr RMSE [piksel] RMSE [piksel] RMSE [piksel] RMSE [piksel] več. 207/194 90 7/10 0,75 0,85 0,76 0,82 Radgona man. 50/46 90 2/5 0,94 0,97 0,95 0,91 več. 202/187 90 6/10 0,95 0,95 0,88 0,86 Koper man. 56/52 90 3/5 0,92 0,99 0,85 1,01 več. 101/94 90 5/10 1,38 1,09 1,21 0,93 Bohinj man. 29/25 80 3/5 0,94 1,15 0,94 1,11 Pri posnetku Koper je RANSAC v večjem nizu odstranil šest točk, pri manjšem pa tri točke z napakami (slika 17). Med izločenimi točkami ni nobene od tistih, 98 ki smo jih dobili s podatki OpenStreetMap, kar kaže, da lahko tudi s temi podatki dobimo kakovostne oslonilne točke. Podobne ugotovitve veljajo tudi za druga dva posnetka (priloga A). Tudi v tem primeru na končno natančnost rezultatov dodane napake ne vplivajo in odstopanja pokvarjenih točk odražajo dodane grobe napake. Preglednica 15: Velikost odstopanj oslonilnih točk z dodanimi napakami, ki so ostale po izravnavi z algoritmom RANSAC za večji in manjši niz točk posnetkov RapidEye. Groba napaka Večji niz Manjši niz Posnetek [piksel] dx [piksel] dy [piksel] dx [piksel] dy [piksel] 1 0,05 1,07 -1,48 0,68 Radgona 2 -2,24 1,73 -1,96 1,49 3 -3,34 2,98 -2,72 2,24 1 -1,11 1,09 -0,86 0,52 2 -2,52 1,60 -2,52 1,09 Koper 3 -3,91 2,39 - - 4 -3,52 3,29 - - 1 -0,32 0,59 -1,69 1,07 2 -2,07 2,28 - - Bohinj 3 -2,78 2,44 -2,01 1,10 4 -4,63 2,96 - - 5 -4,56 3,15 - - Slika 16: Rezultati izravnave (RMSE na kontrolnih točkah) po izravnavi z algoritmom RANSAC (črtkana črta) in na koncu postopka (polna črta) za večji in manjši niz točk posnetkov RapidEye. Iz slike je razvidno, da v splošnem robustna ocena izboljša natančnost rezultatov izravnave. 99 100 Slika 17 (na prejšnji strani): Porazdelitev in odstopanja oslonilnih in kontrolnih točk. Na zgornji sliki je prikazan posnetek Koper (RapidEye). Oslonilne točke so prikazane z rumeno barvo, kontrolne pa z zeleno. Na spodnji sliki sta neodstranjeni točki z dodanimi napakami označeni z rdečima elipsama. Rezultati za posnetek Bohinj so podobni prejšnjim primerom, vendar so opazne nekatere razlike. Po izločanju v večjem nizu, ki je na začetku imel 101 točko, je ostalo pet oslonilnih točk z napakami. Njihova končna odstopanja odražajo dodane napake z nekaj manjšimi odkloni. Nepričakovani rezultati so se pojavili v manjšem nizu, kjer samo točki z enim in tremi piksli dodanih napak nista bili izločeni v prvem delu izravnave. Poleg tega imata točki manjša odstopanja od pričakovanih. Razlog za tako obnašanje modela je verjetno v manjšem številu točk v nizu (29) in zelo razgiban teren, ki otežuje enakomerno razporeditev oslonilnih točk. Kljub temu je natančnost rezultatov po izravnavi še vedno primerljiva z natančnostjo rezultatov brez dodanih napak. 5.6.3.2 Poskusi s posnetki WorldView-2 Pri preizkušanju izločanja grobih napak smo uporabili enake začetne nastavitve in metode kot pri prejšnjih poskusih. Rezultati izravnave so prikazani v preglednicah 16 in 17. Zaradi nezanesljivega delovanja geometričnega modela pri posnetku Ljubljana smo poskus izločanja grobih napak izvajali samo na posnetkih Jesenice in Bled. Ker je bilo za posnetek Jesenice na voljo manj točk, smo v večjem nizu dodali napake šestim točkam, v manjšem pa trem. S tem smo lahko ohranili mejo izločanja točk pri 10 %. V obeh nizih je algoritem RANSAC pustil le točki z dodano napako enega in dveh pikslov. Čeprav je bila RMSE po izravnavi z algoritmom RANSAC nekoliko večja, so se natančnosti po izravnavi z robustno oceno izboljšale (slika 19) in se lahko primerjajo z natančnostjo rezultatov brez dodanih napak (preglednica 11). Postopek robustne ocene je torej v teh primerih nujen za dosego visoke točnosti. Preglednica 17 prikazuje odstopanja točk z napakami, ki so ostale v nizu uporabljenih točk do konca delovanja geometričnega modela. Vrednosti so precej podobne dodanim napakam, kar pomeni, da jih je geometrični model prepoznal in omejil njihov 101 vpliv na rezultate. Uspešno prepoznavanje je predvsem posledica večjega števila točk in robustne ocene. Ugotovljeni zaključki veljajo tako za večji kot manjši niz. Slika 18: Porazdelitev in odstopanja oslonilnih in kontrolnih točk. Na zgornji sliki (posnetek Bled – WorldView-2) so oslonilne točke prikazane z rumeno barvo, kontrolne pa z zeleno. Na spodnji sliki so neodstranjene točke z dodanimi napakami označene z rdečimi elipsami. 102 Preglednica 16: Rezultati izravnave s številom prisotnih in izločenih oslonilnih točk z dodanimi napakami, številom vseh in uporabljenih oslonilnih točk ter RMSE po izravnavi z algoritmom RANSAC in na koncu postopka (robustna ocena) za posnetke WorldView-2. Zaradi slabega delovanja geometričnega modela poskusi s posnetkom Ljubljana niso bili izvedeni. RANSAC Robustna ocena [%] ke Posnetek Niz e/ ene ločanja ablj OT KT OT KT isotne napa Vse točk Upor točke Meja iz Izločene napake/ Pr RMSE [piksel] RMSE [piksel] RMSE [piksel] RMSE [piksel] več. 75/68 90 4/6 1,74 2,74 1,06 1,56 Jesenice man. 30/28 90 1/3 1,58 2,65 1,18 1,94 več. 158/145 90 7/10 1,25 2,19 1,10 2,81 Bled man. 50/47 90 1/5 1,80 4,81 1,40 3,97 - - - - - - - - Ljubljana - - - - - - - - Preglednica 17: Velikost odstopanj oslonilnih točk z dodanimi napakami, ki so ostale po izravnavi z algoritmom RANSAC za večji in manjši niz točk posnetkov WorldView-2. Zaradi slabega delovanja geometričnega modela poskusi s posnetkom Ljubljana niso bili izvedeni. Groba napaka Večji niz Manjši niz Posnetek [piksel] dx [piksel] dy [piksel] dx [piksel] dy [piksel] 1 -1,10 0,59 -0,66 -0,36 Jesenice 2 -2,20 2,83 -1,91 1,13 1 -1,64 0,97 -0,99 -0,84 2 -1,80 0,81 -1,59 2,35 Bled 3 -3,57 1,44 -1,31 2,78 4 - - -3,76 3,12 Ljubljana - - - - - Pri posnetku Bled smo imeli na voljo več točk in smo lahko večjemu nizu dodali napake desetim, manjšemu pa petim točkam. V obeh nizih je RANSAC v večini primerov izločil vse točke z napakami razen treh (slika 18). Pri manjšem nizu se je le enkrat zgodilo, da je bila izločena le ena točka z napakami, ostale so pa štiri. Tudi pri tem posnetku je robustna ocena dobro delovala in zmanjšala vpliv točk z dodanimi napakami. To je vidno tudi pri odstopanjih oslonilnih točk od končnega modela, ki odražajo velikost 103 dodanih grobih napak. Manjši odkloni odstopanj od pričakovanih so posledica manjših napak določanja točk in vpliva ostalih točk z napakami. K napakam prispevata tudi razgiban teren in slabša razporeditev točk. Kljub temu je točnost rezultatov po izravnavi še vedno primerljiva s točnostjo rezultatov brez dodanih napak. Slika 19: Rezultati izravnave (RMSE kontrolnih točk) po izravnavi z algoritmom RANSAC (črtkana črta) in na koncu postopka (polna črta) za večji in manjši niz točk posnetkov WorldView-2. Iz slike je razvidno, da v splošnem robustna ocena izboljša natančnost rezultatov izravnave. 104 6 ORTOREKTIFIKACIJA Izdelan postopek se konča z ortorektifikacijo posnetka. V tem poglavju so najprej predstavljene metode za izdelavo ortofotov, ki se delijo na dve skupini: direktne (neposredne) in indirektne (posredne). Opisani so postopki obdelave in rezultati posameznih metod. Sledi še podroben opis indirektne metode ortorektifikacije, ki je bila vključena v postopek samodejne ortorektifikacije. Ortorektifikacija je pretvorba posnetka v ortogonalno projekcijo (glej poglavje 2.3). Z ortorektifikacijo določimo, kateri slikovni element na surovem posnetku ustreza slikovnemu elementu na ortofotu oz. položaju v prostoru. Ker v samodejnem postopku obdelujemo le en posnetek naenkrat, ne moremo uporabiti standardne fotogrametrične metode prostorskega preseka žarkov (angl. space intersection), ki se uporablja v primeru stereoparov. V primeru enega posnetka se uporablja DMR in pa pogoj kolinearnosti, ki rigorozno poveže slikovno koordinato s koordinato v prostoru. 6.1 Metode izdelave ortofotov Običajno se pri ortorektifikaciji uporabljata dve skupini metod (Konecny 1979). V prvo skupino sodijo t.i. direktne metode, ki delujejo tako, da vsak piksel oz. slikovne koordinate projicirajo v trirazsežni prostor (3R). Druga skupina, t.i. indirektne metode, delujejo v obratni smeri. Metode začnejo v objektnem prostoru (angl. object space) in projicirajo položaj slikovnega elementa preko DMR-ja na surov posnetek oz. slikovne koordinate. Radiometrično vrednost na ortofotu dobimo s prevzorčenjem radiometričnih vrednosti okoliških pikslov na surovem posnetku. 105 6.1.1 Direktne metode Direktne metode se veliko uporabljajo za določitev posameznih prostorskih koordinat, redkeje za celotno ortorektifikacijo. Za ortorektifikacijo so služile predvsem pri vrstičnih senzorjih v osemdesetih letih (npr. za satelit SPOT), preden so se uveljavili indirektni modeli (Chen in Lee 1993). Tudi v primeru direktnih metod potrebujemo poleg parametrov zunanje orientacije še informacijo o višini, ki jo dobimo iz DMR-ja. Prostorske koordinate piksla dobimo s presekom med površjem in žarkom gledanja (angl. view ray), ki se začne v gorišču (perspektivnem centru) in gre skozi obravnavani piksel. Metode, s katerimi se na ta način določa koordinate preseka, imenujemo tudi metode za obrnjeno projiciranje posamičnega žarka (angl. single-ray backprojection) (Mikhail in dr. 2001). V splošnem poznamo dve metodi, s katerima rešujemo problem obrnjenega projiciranja (Sheng 2004): iterativna fotogrametrična metoda (angl. iterative photogrammetric method) in metoda sledenja žarku (angl. ray tracing method). Iterativna fotogrametrična metoda (IP) se najpogosteje uporablja v fotogrametriji, medtem ko se metoda sledenja žarku (RT) največ uporablja v računalniški grafiki. V praksi se je uveljavila tudi tretja metoda, imenovana iterativna metoda sledenja žarku (IRT), ki sicer temelji na metodi RT, vendar je od nje hitrejša in učinkovitejša. V vseh naštetih primerih se uporablja inverzne kolinearne enačbe, kjer iščemo prostorski koordinati X in Y (Mikhail in dr. 2001):         X  X   S  Z ZS m x x m y y m f 11  0  21  0  31   m  x  x  m  y  y  m   f 13  0  23  0  33           (73) Y  Y   S  Z ZS m x x m y y m f 12  0  22  0  32   m  x  x  m  y  y  m   f 13  0  23  0   , 33 pri čemer so parametri enaki kot pri enačbi (29). 6.1.1.1 Iterativna fotogrametrična metoda V fotogrametriji je najpogosteje uporabljena iterativna fotogrametrična metoda, ki deluje iterativno in je od vseh najhitrejša. Slika 20 prikazuje njeno delovanje, ki temelji na izračunu preseka slikovnega žarka z DMR-jem. 106 Postopek, s katerim metoda pridobi prostorske koordinate točke T, deluje v naslednjih korakih: 1. Z inverznimi kolinearnimi enačbami in povprečno višino modela reliefa ( Zp) se izračunajo koordinate začetne projicirane točke 1. 2. Iz točke 1 sledimo do točke 2 na površju modela reliefa, kjer izračunamo višino. 3. Izračunano višino uporabimo za izračun prostorskih koordinat točke 3. 4. Izračune iz prvih treh točk ponavljamo dokler koordinate ne konvergirajo. 5. Če sta zadnji dve projicirani točki znotraj predhodno določenih meja, se iteracije ustavijo in zadnje koordinate obveljajo kot končni rezultat. 6. Radiometrična vrednost iz posnetka se prepiše na izračunane koordinate v ortofotu. Slika 20: Iterativna fotogrametrična metoda (IP). Metoda IP je v primerjavi z drugimi metodami zelo hitra in dokaj zanesljiva, v kolikor pri ortorektifikaciji uporabljamo DMR srednje ali nizke ločljivosti. Pri zelo razgibanem terenu in visokoločljivem modelu reliefa pa lahko odpove ali 107 privede do napačnih rezultatov. Napake nastopijo zaradi divergence in zakrivanja (angl. occlusion), ki sta med drugim posledica velikega kota gledanja in naklona površja. V redkih primerih se lahko postopek ujame v neskončno zanko. Zaradi tega uporaba modelov reliefa iz lidarskih podatkov v polni ločljivosti, ki je običajno en meter in manj, ni priporočljiva. 6.1.1.2 Iterativna metoda sledenja žarku Ker je metoda sledenja žarku procesno zelo zahtevna in počasna, se je v fotogrametriji uveljavila njena iterativna različica, ki dosega enake rezultate v krajšem času. Iterativna metoda sledenja žarku, ki sta jo predlagala O'Neill in Dowman (1988), določi točko preseka z iterativnim podaljševanjem žarka gledanja, dokler ne doseže površja. Ključni parameter je dolžina koraka, s katerim podaljšamo žarek, saj je od njega odvisna hitrost izračuna koordinat in točnost rezultata. Če je korak velik, je metoda zelo hitra, vendar lahko pride do napak. V primeru kratkega koraka je postopek bolj zamuden, rezultati pa so zelo točni. Postopek, s katerim metoda pridobi prostorske koordinate točke T, deluje v naslednjih korakih (slika 21): 1. Z inverznimi kolinearnimi enačbami in največjo višino modela reliefa ( Zmaks) se izračunajo koordinate začetne točke. 2. Z izračunanimi koordinatami se pridobi pripadajočo višino iz modela reliefa. 3. Izračuna se razliko med zadnjimi višinami in se jo primerja s preddoločeno mejo. Če je razlika večja od meje, se žarek gledanja podaljša za korak (višini se vrednost koraka odšteje). 4. Izračune iz prvih treh točk ponavljamo, dokler ni razlika višin negativna ali pade znotraj meje. 5. V primeru negativne razlike se korak zmanjša ter prišteje zadnji višini. 6. Izračun se ponavlja, dokler razlika višin ne pade znotraj meje. V tem primeru se iteracije ustavijo in zadnje koordinate obveljajo kot končni rezultat. 7. Radiometrična vrednost iz posnetka se prepiše na izračunane koordinate na ortofotu. 108 Iterativna metoda sledenja žarku sicer ni tako hitra kot iterativna fotogrametrična metoda, vendar je veliko bolj zanesljiva in pri izbranem majhnem koraku tudi ne prihaja do napak. Zanesljiva je tudi pri modelih reliefa visoke ločljivosti in pri velikih kotih gledanja satelita ter se uporablja za izdelavo popolnega ortofota (angl. true orthophoto). Slika 21: Iterativna metoda sledenja žarku (IRT). Direktne metode so ob pravilni uporabi pri izdelavi ortofotov zelo zanesljive. Njihova največja pomanjkljivost je računska zahtevnost in trajanje obdelave. Končni ortofoto ima ponavadi tudi prazne piksle, saj se zaradi različnih velikosti pikslov in predvsem razgibanosti reliefa vsak piksel ne preslika vedno v svoj element ortofota (slika 22). V tem primeru je potrebno piksle naknadno zapolniti, metoda zapolnitve pa je odvisna od velikosti nastalih lukenj. Če so te manjše, se lahko uporabi prevzorčenje iz sosednjih pikslov, drugače potrebujemo dodatne posnetke, posnete iz drugačnih zornih kotov, s katerimi zapolnimo manjkajoče dele. Izdelava ortofotov z direktno metodo iz posnetka vrstičnega skenerja je prikazana na sliki 23. Zaradi reliefa, kota gledanja in orientacije senzorja med snemanjem je število pikslov in tudi njihove vrednosti nekoliko drugačno od surovega posnetka, saj 109 med transformacijo (ortorektifikacijo) vse piksle prevzorčimo, nekatere lahko celo združimo. V primeru lukenj moramo vrednost praznih pikslov pridobiti iz drugih virov. Slika 22: Ponazoritev preslikave pikslov pri izdelavi ortofota z direktnimi metodami. Zaradi preslikave pikslov (T6 in T8) v sosednje celice sta dve celici v ortofotu prazni (P1 in P2). Slika 23: Direktna metoda ortorektifikacije za vrstične skenerje. 110 6.1.2 Indirektne metode Indirektne metode se v fotogrametriji pojavijo že v šestdesetih letih prejšnjega stoletja (Smith 1995), digitalne oblike metod pa v sedemdesetih (Konecny 1979). Pred uvedbo vrstičnih skenerjev oz. visokoločljivih DMR-jev so bile tudi edine praktično uporabljene metode za ortorektifikacijo. Zaradi hitrosti obdelave in enostavnosti algoritmov se jih tudi danes veliko uporablja. Indirektne metode delujejo v obratni smeri kot direktne metode - s projiciranjem točke iz prostorskega sistema v točko na posnetku, ki je v 2R- sistemu. Za vrstične skenerje so se prve indirektne metode začele pojavljati razmeroma pozno. Chen in Lee (1993) sta predlagala rešitev, ki temelji na metodi Newton-Raphson. Čeprav je ta metoda uspešna, je zelo občutljiva na začetne pogoje in pravilno deluje le, kjer se enačba tirnice spreminja monotono. To sicer velikokrat drži za satelite, redko pa za letalska snemanja. Izboljšano verzijo indirektne metode so predlagali Kim in dr. (2001). Metoda je bolj robustna, saj deluje tudi pri nemonotonem premikanju in rotacijah platforme, na kateri je nameščena kamera. V naštetih primerih se uporablja enačbe (29), kjer iščemo slikovni koordinati x in y. Izdelava ortofota z indirektno metodo za polnoslikovni senzor je dokaj preprosta. Zaradi zaporednega snemanja vrstic ob različnih časih pa pri vrstičnih skenerjih transformacija slikovne koordinate iz prostorske ni reverzibilna. To pomeni, da lahko rešitev dobimo z iteracijami ali dodatnimi izračuni, odvisno od geometričnega modela. V splošnem pa indirektne metode delujejo v naslednjih korakih: 1. Izračuna se velikost in koordinate ortofota ter določi velikost piksla. 2. Za vsak piksel se pridobi pripadajočo višino iz modela reliefa. 3. S kolinearnimi enačbami, 2R-koordinatami pikslov in višinami modela reliefa na pikslih se izračunajo slikovne koordinate na surovem posnetku. 4. Z interpolacijo radiometričnih vrednosti v okolici izračunane koordinate se dobi vrednost piksla na ortofotu. 5. Na enak način se izračuna vrednosti vsakega piksla na ortofotu, kamor se vrednosti tudi zapišejo. 111 Indirektne metode imajo več prednosti v primerjavi z direktnimi. Dokazano je, da so indirektne metode boljše od direktnih v kakovosti ter hitrosti obdelave (Kim in dr. 2001). Kljub vsem prednostim pa so lahko indirektne metode neuporabne pri uporabi visokoločljivih modelov reliefa in/ali pri obdelavi posnetkov, ki so bili zajeti pod velikim kotom. Poleg tega indirektne metode ne puščajo praznih pikslov, zaradi česar dodatna obdelava ni potrebna. Pri posnetkih, ki so snemani v nadirju in v vrstičnem načinu, je uporaba indirektnih metod zanesljiva. Izdelava ortofota z indirektno metodo iz polnoslikovnega posnetka je prikazana na sliki 24. Tudi v primeru indirektnih metod je zaradi reliefa, kota gledanja in orientacije senzorja med snemanjem število pikslov in tudi njihove vrednosti nekoliko drugačno od surovega posnetka, saj med transformacijo (ortorektifikacijo) vse piksle prevzorčimo, nekatere lahko celo združimo. Slika 24: Ortorektifikacija polnoslikovnega posnetka z indirektno metodo. 112 6.2 Postopek izdelave ortofotov Indirektne metode so hitre in zanesljive, če ne gre za ekstremne primere snemanja. Zaradi tega smo to metodo uporabili tudi v okviru izdelanega postopka. V primeru polnoslikovnega sistema je uporaba enostavna, pri vrstičnih senzorjih pa je prisotna še časovna komponenta, ki zahteva dodatne izračune. Slika 25 prikazuje potek indirektne metode ortorektifikacije v primeru razvitega geometričnega modela za vrstični skener. Pri izračunu slikovnih koordinat je potrebno upoštevati delitev tirnice na dva dela ter parametre zunanje orientacije, ki jih dobimo po enačbah (34). Slika 25: Indirektna metoda ortorektifikacije za vrstični skener. Pri ortorektifikaciji je treba upoštevati geometrični model, v katerem je tirnica razdeljena na dva segmenta (S1 in S2). 113 Ker je pri vrstičnem skenerju slikovna koordinata x enaka 0 za vsako vrstico, lahko kolinearne enačbe zapišemo (Kim in dr. 2001): m  11  X  X       S  m 12  Y YS  m 13  Z ZS  0  f  m 31 X  X       S  m 32  Y YS  m 33  Z ZS  (74) m  21  X  X       S  m 22  Y YS  m 23  Z ZS  y  f  m 31 X  X       S  m 32  Y YS  m 33  Z S  , Z pri čemer so parametri enaki kot pri enačbi (29). Torej lahko problem inverzne metode rešimo, če izračunamo koordinato x iz enačbe: 0  m  . (75) 11  X  X  m  Y  Y  m  Z  Z S  12  S  13  S  Za razliko od perspektivnih posnetkov (en projekcijski center za celoten posnetek) pri vrstičnih posnetkih elementi rotacije in položaj satelita niso konstantni. Zaradi tega enačba (75) ni linearna in je ne moremo rešiti analitično (Kim in dr. 2001). Če v enačbo (75) vstavimo enačbe parametrov zunanje orientacije, dobimo: f ( t ) 2  t ( X  m  Y  m  Z  m )  t ( X  m  Y  m  Z  m )  2 11 2 12 2 13 1 11 1 12 1 13 (76)  m  X  X  m  Y  Y  m  Z  Z 11  0  12  0  13  0 . Z rešitvijo kvadratne enačbe dobimo čas snemanja t za posamezno točko, pri čemer je f(t) = 0. Celoten postopek izračuna slikovnih koordinat točke je naslednji: 1. Izberemo 3R-točko (koordinate X, Y, Z). Ker ni znano v katerem segmentu se nahaja projicirana slikovna točka, se čas snemanja t izračuna dvakrat, enkrat z vsakim nizom parametrov zunanje orientacije. 2. Glede na izračunana trenutka snemanja se lahko ugotovi, kateremu segmentu pripada točka, s čimer določimo čas snemanja točke. 3. Izračuna se pripadajoče parametre zunanje orientacije ob času snemanja točke. 4. Slikovno koordinato x dobimo po enačbi: 114 t x  , (77) tv kjer je tv čas snemanja ene vrstice. 5. Slikovno koordinato y v metrih dobimo najprej po kolinearni enačbi. Končno pa z: y n m s y   , (78) p 2 vel kjer je ym slikovna koordinata y v metrih, pvel velikost piksla v metrih, ns pa število stolpcev posnetka. 6. Z bilinearno interpolacijo radiometričnih vrednosti pikslov v okolici izračunane slikovne koordinate se dobi vrednost piksla na ortofotu. Celoten postopek se opravi za vsak piksel ortofota. Na koncu se rezultatu pripiše še informacija o položaju (koordinate zgornjega levega vogala ortofota) in ortofoto shrani v zahtevan zapis. 6.3 Implementacija Izdelava ortorektificiranega posnetka se izvaja z indirektno metodo v samostojnem procesnem modulu za izvedbo ortorektifikacije. Glavni izdelek zadnjega modula je ortofoto. 6.3.1 Modul za izvedbo ortorektifikacije Podatkovna struktura z izračunanimi parametri se prenese v modul za izvedbo ortorektifikacije, ki izdela končni ortofoto. Algoritem najprej prebere surov posnetek z vsemi kanali. Prvi korak ortorektifikacije je določitev velikosti matrike končnega ortofota. Ker ima večina satelitov v sončno-sinhroni tirnici inklinacijo okoli 98º, je večina posnetkov zasukanih v smeri urinega kazalca glede na geografski sever. Zaradi tega največkrat vsi štirje vogalni piksli posnetka določajo največji in najmanjši obseg ortofota (slika 26). Zato se v modulu najprej izračuna prostorske koordinate vogalov končne matrike. Pri tem se uporabi direktno metodo iterativnega sledenja žarku, ki sicer ni najhitrejša, je pa najbolj zanesljiva in splošno uporabna. Izračun poteka z 115 inverznimi kolinearnimi enačbami, začetna vrednost višine pa je največja vrednost uporabljenega modela reliefa. Iz izračunanih koordinat se z odštevanjem dobi velikost območja posnetka. Z deljenjem z ločljivostjo posnetka pa dobimo velikost matrike ortofota, izraženo v pikslih. Tako izračunane meje ne veljajo vedno za posnetke zajete pod velikim kotom, tiste, ki prikazujejo gorata območja, ali posnetke satelitov z zmožnostjo snemanja v smeri poldnevnika. Zaradi tega se velikost matrike poveča za 200 pikslov v obeh smereh. Na koncu postopka ortorektifikacije se prazne piksle odstrani in tako zmanjša velikost datoteke. Pri izdelavi ortofota se v trenutni različici modula uporablja DMR z ločljivostjo 12,5 m. Ker model reliefa ni visokoločljiv in večina satelitov snema pod majhnimi koti, se za ortorektifikacijo lahko uporabi veliko različnih metod. Uporaba drugega DMR-ja v razvitem postopku je preprosta, paziti je treba le na uporabljeno metodo glede na prostorsko ločljivost podatkov. Za izdelan postopek smo izbrali indirektno metodo ortorektifikacije (glej podpoglavje 6.1.2), saj ima nizke zahteve glede strojne opreme, je zelo hitra in ustreza uporabljenemu DMR-ju. Slika 26: Določanje velikosti ortofota (matrike) glede na surovi posnetek. Ortorektifikacija zaporedno obdeluje piksle za vsako vrstico vhodnega posnetka. Najprej obdela zgornji levi piksel in nadaljuje po vrstici do zadnjega piksla. Vsakemu pikslu oz. prostorski koordinati piksla ( X in Y) se najprej določi višina na DMR-ju. Višino se izračuna z bilinearnim prevzorčenjem sosednjih pikslov. Nato sledi postopek interpolacije, ki je opisan v poglavju 6.2. Rezultat so slikovne koordinate piksla. Če sta koordinati znotraj meja surovega posnetka, se radiometrično vrednost piksla izračuna z bilinearno interpolacijo sosednjih pikslov na surovem posnetku in nato zapiše v ortofoto. Če sta koordinati izven meja posnetka, ostane piksel prazen oziroma se mu pripiše 116 vrednost 0. Dobljeno matriko radiometričnih vrednosti, ki jo praviloma tvorimo v bralno-pisalnem pomnilniku (angl. random-access memory, RAM), se na koncu izpiše na trdi disk v zapisu GeoTIFF. Datoteke GeoTIFF se lahko bere v vseh geografskih informacijskih sistemih in v večini grafičnih programov. 6. 4 Testiranje Uspešnost in robustnost delovanja razvitega modula za ortorektifikacijo smo preverili z empiričnim poskusom na obeh vrstah satelitskih posnetkov. 6.4.1 Ocena položajne točnosti ortofotov Ortofote smo izdelali s parametri zunanje orientacije najbolj natančnega rezultata in modelom reliefa z ločljivostjo 12,5 m (GURS 2017a). Natančnost ortofotov smo preverili z ročno izmerjenimi kontrolnimi točkami. Na vsakem posnetku smo izmerili 30 kontrolnih točk, ki smo jih postavili na cestna križišča, enakomerno po celem posnetku. Koordinate točk na istih položajih smo izmerili še na državnem ortofotu z ločljivostjo 0,5 m in ocenjeno položajno točnostjo 1 m. S primerjavo koordinat smo lahko ocenili položajno točnost ortofotov. 6.4.1.1 Poskusi s posnetki RapidEye Ortofote RapidEye smo izdelali s parametri zunanje orientacije, ki smo jih dobili pri najnatančnejši rešitvi modela, in z DMR-jem z ločljivostjo 12,5 m. Ortorektificiran posnetek Koper je prikazan na sliki 28. Točnost ortofotov smo ocenili z ročno izmerjenimi kontrolnimi točkami. Na vsakem ortofotu smo izbrali 30 enakomerno porazdeljenih točk, ki so ležale na križiščih cest. Iz državnega ortofota z ločljivostjo 0,5 m in ocenjeno točnostjo pod 1 m smo pridobili položajne koordinate istih lokacij. Točnost (RMSE) ortorektificiranih posnetkov smo dobili s primerjavo obeh nizov koordinat. Preglednica 18 prikazuje položajno točnost (RMSE) ortofotov. Dobljene napake RMSE (slika 27) kažejo na visoko podobnost z rezultati, dobljenimi z geometričnim modelom na kontrolnih točkah (preglednica 10). Posnetek 117 Radgona, ki prikazuje večinoma ravninski teren, je dosegel najboljšo natančnost (0,85 piksla), medtem ko je RMSE preostalih posnetkov okoli piksla oz. 6,5 m. Izračunana natančnost je sicer odraz le izmerjenih točk z dobro vidnih križišč, ki pa se večinoma nahajajo na nižjem, ravninskem terenu. Zaradi tega položajne natančnosti gorskih območij (sredinski del posnetka Bohinj) nismo mogli ustrezno oceniti. Točnost takih območij bi lahko ocenili s postavljanjem točk na gorske grebene, kar pa je težko zaradi prisotnosti snega in slabše točnosti DMR-ja na gorskih območjih. Preglednica 18: Točnost izdelanih ortofotov RapidEye glede na državni ortofoto, ki je bil uporabljen kot referenca. RMSE [piksel] Posnetek Število kontrolnih točk X Y Skupni Radgona 30 0,69 0,49 0,85 Koper 30 0,73 0,77 1,06 Bohinj 30 0,60 0,80 1,01 Slika 27: RMSE na 30 ročno določenih kontrolnih točkah za posnetke RapidEye. 118 Slika 28: Ortorektificiran posnetek Koper (RapidEye). 6.4.1.2 Poskusi s posnetki WorldView-2 Ortofote WorldView-2 smo izdelali s parametri zunanje orientacije, ki smo jih dobili pri najnatančnejši rešitvi modela, in z digitalnim modelom reliefa z ločljivostjo 12,5 m. Ortorektificiran posnetek Bled je prikazan na sliki 29. Natančnost ortofotov smo preverili z ročno izmerjenimi kontrolnimi točkami. Na vsakem ortofotu smo izbrali približno 30 enakomerno porazdeljenih točk, ki so ležale na križiščih cest. Iz državnega ortofota z ločljivostjo 0,5 m smo pridobili položajne koordinate istih lokacij. Položajno točnost (RMSE) ortorektificiranih posnetkov smo dobili s primerjavo obeh nizov koordinat. Preglednica 19 prikazuje položajno točnost (RMSE) ortofotov. Dobljene napake RMSE (slika 30) so celo nižje od povprečnih vrednosti rezultatov, ki so 119 bili dobljeni z geometričnim modelom na kontrolnih točkah (preglednica 11). Pri tem moramo upoštevati, da smo pri ortorektifikaciji uporabili najboljše rezultate izravnav in ne povprečja rezultatov. Pričakovano je najbolj točen posnetek Jesenice z RMSE 1,16 piksla. Ostala posnetka sta nekoliko slabše položajne točnosti in imata RMSE okoli 2 piksla. Točnost zelo visokoločljivih ortofotov je odvisna tudi od ločljivosti in točnosti uporabljenega DMR-ja. Če bi imeli pri ortorektifikaciji na razpolago boljši DMR (npr. lidarski), bi lahko dosegli še višjo točnost ortopodob. Točnost je sicer odraz le izmerjenih točk, ki smo jih postavili na dobro vidna križišča, ki pa se večinoma nahajajo na nižjem, ravninskem terenu. Zaradi delne lege zunaj Slovenije, kjer nismo imeli kontrolnih podatkov, je bila točnost posnetka Jesenice ocenjena samo za spodnji del posnetka. Slika 29: Ortorektificiran posnetek Bled – WorldView-2. 120 Preglednica 19: Točnost izdelanih ortofotov WorldView-2 glede na državni ortofoto, ki je bil uporabljen kot referenca. RMSE [piksel] Posnetek Število kontrolnih točk X Y Skupni Jesenice 30 0,80 0,84 1,16 Bled 32 0,69 1,78 1,91 Ljubljana 32 1,03 1,75 2,03 Slika 30: RMSE na 30 ročno določenih kontrolnih točkah za posnetke WorldView-2. V primeru zelo visokoločljivih posnetkov je položajna točnost ortofotov tudi odraz točnosti referenčnih podatkov. Čeprav smo državni ortofoto uporabili za referenco pri določanju točk, lahko vsebuje naključne napake, ki so lahko velikosti 1 m ali več. Zaradi tega ima pri tako visoki ločljivosti, kjer je prisotno tudi veliko detajlov, geometrični model problem pri izračunu parametrov oz. modeliranju točne geometrije snemanja. Dodatni problem je tudi način snemanja (npr. asinhrono), ki v modelu ni upoštevan. Pri ločljivosti 2 m je sicer vpliv referenčnih podatkov na točnost manjši kot pri pankromatskih posnetkih (0,5 m), kjer bi bila za oceno točnosti nujna primerjava s terensko izmerjenimi točkami (npr. izmera z GPS). Za primere naših poskusov lahko položajno točnost izdelanih ortofotov generalno ocenimo okrog vrednosti dveh pikslov, kar je v skladu s pričakovanji. 121 7 ZAKLJUČKI V okviru raziskave smo razvili nov samodejen postopek za ortorektifikacijo optičnih satelitskih posnetkov ravni 1A. Postopek izdela satelitske ortofote v ortogonalni projekciji v približno eni uri in na povsem samodejen način, brez posredovanja operaterja. Ortorektificiran posnetek je v državnem koordinatnem sistemu in lahko služi za primerno podlago v prostorskih analizah in kot vir za kartiranje (npr. interpretacijo in vektorizacijo podatkov). Zaradi (običajno) večje spektralne in časovne ločljivosti lahko v določenih aplikacijah dopolnjuje ali celo nadomesti državni ortofoto. Večina dosedanjih postopkov ortorektifikacije je ročnih ali polsamodejnih, oziroma delujejo le za določen satelit ali le v določenih pogojih (glej podpoglavje 1.1.3). Naš izdelan postopek ortorektifikacije pa povezuje več različnih metod v enoten in robusten sistem za samodejno izdelavo ortofotov tako iz posnetkov pridobljenih z vrstičnimi kot polnoslikovnimi senzorji. Celoten postopek je sestavljen iz štirih osnovnih modulov: modula za branje in pripravo metapodatkov, modula za samodejno določanje oslonilnih točk, modula za izračun parametrov geometričnega modela in modula za izvedbo ortorektifikacije. Modul za branje in pripravo metapodatkov prebere metapodatkovno datoteko in izbere tiste parametre, ki se kasneje uporabijo v postopku za določanje točk in ortorektifikacijo posnetka. Tako pridobi podatke o geometriji senzorja, lastnostih posnetka in začetne približke parametrov zunanje in notranje orientacije. Drugi korak modula služi za preračun prebranih podatkov v koordinatni sistem, ki se uporablja v geometričnem modelu in v katerem bo tudi končni ortofoto. Pridobljeni podatki vstopajo v model kot približne vrednosti iskanih parametrov (neznanke). Tudi modul za samodejno določanje oslonilnih točk je razdeljen na dva zaporedna koraka. Prvi korak za določitev točk uporablja digitalne 122 podatke o cestah, ki so sestavljeni iz treh virov (podatke o cestah iz DTK5, iz ZKGJI in iz podatkov OpenStreetMap), in zaznane ceste iz surovega posnetka. Drugi korak položaj točk še izboljša s slikovnim ujemanjem med izseki surove podobe in izseki iz državnega ortofota. Pri visokoločljivih posnetkih že prvi korak poda zanesljive točke, medtem ko je za zelo visokoločljive posnetke potreben še drugi korak. Modul za izračun parametrov geometričnega modela je najbolj kompleksen del postopka. V njem izračunamo parametre notranje in zunanje orientacije posnetkov, katerih natančnost določitve vpliva na končno položajno točnost ortofota. Zaradi reševanja problema izločanja grobih napak, ki se jim v samodejnem postopku ne moremo izogniti, je postopek računanja parametrov geometričnega modela razdeljen na dva dela: prvi deluje z algoritmom RANSAC, drugi pa z robustno oceno. Zadnji modul izdeluje ortofote z uporabo digitalnega modela reliefa z ločljivostjo 12,5 m. Pri tem uporabi indirektni model ortorektifikacije, ki je bil izbran zaradi zanesljivosti in hitrosti delovanja. Realizirani postopek ortorektifikacije je torej popolnoma samodejen, vključuje obdelavo tako vzdolžnih kot polnoslikovnih senzorjev in pri tem uporablja le dostopne državne podatke in prostodostopne podatke o cestah. Glavni in osrednji del raziskave je postopek računanja parametrov geometričnega modela. Jedro modela je bilo izbrano iz že obstoječih modelov, ki so bili razviti za določene senzorje in so večinoma delovali v ročnih postopkih. V raziskavi smo jedro dodelali, dopolnili ter mu dodali funkcionalnosti za obdelavo podatkov v samodejnem postopku. Izdelali smo dve izvedbi modela: model vrstičnega senzorja in model polnoslikovnega senzorja. Izvedbi sta indirektni in splošni, kar pomeni, da delujeta le z oslonilnimi točkami in sta prilagojeni za modeliranje geometrije več snemalnih sistemov. Za dodajanje novih sistemov je potrebno prilagoditi le modul za branje in pripravo metapodatkov. Zelo pomemben del postopka je tudi modul za samodejno določanje oslonilnih točk. Modul je bil sicer izdelan za uporabo v samodejnem postopku ortorektifikacije, ni pa bil v celoti razvit v okviru te raziskave. Poleg uporabljene metode korelacije za določitev točk so pomembni tudi referenčni podatki, ki so bili v našem primeru sestavljeni iz podatkov o cestah (sestavljenega iz treh virov) in državnega ortofota. Raziskave so pokazale, da 123 so državni digitalni podatki o cestah z visoko položajno točnostjo (večina pod 1 m, ostalo med 1– 5 m) in digitalni podatki o cestah OpenStreetMap dovolj dobri za določevanje točk na visokoločljivih posnetkih (večinoma podpikselska točnost), ne zadoščajo pa za potrebe zelo visokoločljivih posnetkov. V teh primerih je nujna uporaba dodatnih virov, zato smo v našem primeru uporabili državni ortofoto z ločljivostjo 0,5 m, ki je sicer položajno dovolj točen za posnetke z ločljivostjo 2 m, problem pa je predvsem v različni radiometriji in časovni usklajenosti z vhodnim satelitskim posnetkom. Slikovno ujemanje z državnim ortofotom izboljša položajno točnost točk, zaradi omenjenih neskladij pa je odstopanje v nekaterih primerih še vedno večje od velikosti enega piksla. S poskusi na posnetkih RapidEye in WorldView-2 smo ovrednotili delovanje samodejnega postopka in ocenili položajno točnost ortofotov. Z vsemi posnetki RapidEye smo dosegli zelo dobre rezultate. Rezultati testov so dokazali, da lahko postopek izdela ortofote s položajno točnostjo (RMSE) pod velikostjo enega piksla, tudi če je med podatki prisotnih več grobih napak. Točnost glede na kontrolne točke je bila manj kot piksel tudi v primeru uporabe manj kot 20 oslonilnih točk. Enake rezultate smo dobili tako za ortofote kot za kontrolne točke v geometričnem modelu. Poskusi na posnetkih WorldView-2 kažejo na nekoliko slabšo samodejno določanje oslonilnih točk, kar je vplivalo tudi na končne rezultate. Pri zelo visokoločljivih posnetkih je geometrični model deloval dobro le v predelih z enakomerno razporeditvijo točk. Na območjih, kjer oslonilnih točk ni bilo, so bila odstopanja večja. Nepojasnjeno je slabše delovanje modela na enem posnetku, čeprav je končna točnost ortofota malo nad velikostjo dveh pikslov. V drugih dveh primerih je model pravilno izločal umetno dodane grobe napake in dosegal RMSE ortofotov pod velikostjo dveh pikslov. Rezultati izravnave so sicer nihali, boljši so bili v primeru uporabe več točk. V kolikor bi imeli na razpolago DMR z višjo ločljivostjo, bi verjetno dobili še boljše rezultate. Za robustno oceno so se za najboljše izkazale metode s hiperbolično utežno funkcijo. Čeprav se vrednosti med preizkušenimi metodami le malo razlikujejo, so metode s hiperbolično funkcijo dosegle boljše rezultate, predvsem pa delujejo bolj robustno. Izbira metode robustne 124 ocene sicer nima velikega vpliva na končne rezultate, je pa njena uporaba nujna v primeru, ko je v podatkih veliko grobih napak. Celoten razvit postopek ortorektifikacije je bil zapisan v programskem jeziku IDL. Izjema je prvi korak modula za določanje točk, ki je v jeziku C++. Uporaba razvitih algoritmov (programa) je mogoča na vsakem računalniku z zadostno količino bralno-pisalnega pomnilnika in povezavo do strežnika s pomožnimi podatki. Trenutno algoritmi delujejo le s slovenskim koordinatnim sistemom in podpirajo štiri vrstične senzorje (RapidEye, WorldView-2, Pleiades in SPOT 6) pri čemer dva (Pleiades in SPOT 6) zaradi pomanjkanja satelitskih posnetkov še nista bila temeljito preizkušena in tudi nista predstavljena v tej knjigi. Slabše delovanje postopka v primeru satelita WorldView-2 bi lahko bolje razumeli in izboljšali z dodatnimi posnetki. Za preizkus polnoslikovnega modela nismo uspeli pridobiti posnetkov, lidarski digitalni model reliefa, za preverjanje vpliva DMR-jev različnih ločljivosti na ortorektifikacijo, pa še ni bil na voljo. Izdelan postopek za ortorektifikacijo deluje popolnoma samodejno od sprejetega surovega posnetka do končnega ortofota, ki je umeščen v koordinatni sistem. Operater, razen po potrebi, vidi le končni izdelek s pripadajočimi podatki o obdelavi. Na postopek ni mogoče interaktivno vplivati, kar lahko v določenih primerih da tudi slabše rezultate. Natančnost delovanja geometričnega modela se samodejno preverja samo na koncu izravnav na osnovi RMSE oslonilnih točk. Ta podaja le notranjo natančnost, ki ni pravi pokazatelj uspešnosti delovanja postopka. Absolutno kontrolo modela in ortofotov lahko priskrbi le neodvisen postopek (neodvisno določene kontrolne točke), ki pa ga v samodejnih postopkih ni mogoče zagotoviti. Za skoraj vse aplikacije je pomemben tudi čas poteka celotnega postopka. Ortorektifikacija posameznega posnetka RapidEye (5 kanalov) velikosti 50 km  80 km traja približno eno uro s povprečnim osebnim računalnikom, ki vsebuje procesor Xeon Quad-Core s hitrostjo 2,5 GHz in 8 GB RAM-a. Korak branja metapodatkov traja le trenutek, medtem ko določanje oslonilnih točk poteka približno 38 % celotnega časa. Za posnetke RapidEye je izračun parametrov geometričnega modela ponavadi zelo hiter in traja le 1 % – 2 % 125 časa. Zadnji korak, ortorektifikacija, je računsko najbolj zahteven in zajema preostalih 60 %. V primeru večspektralnega (8 kanalov) posnetka WorldView- 2 velikosti 15 km  17 km obdelava traja približno 70 minut (odvisno predvsem od števila samodejno določenih točk). Korak branja metapodatkov traja le trenutek, medtem ko določanje oslonilnih točk poteka približno 28 % celotnega časa v primeru določitve približno 200 točk. Za posnetke WorldView-2 se izračun parametrov geometričnega modela ponavadi konča v približno 7 % časa. Trajanje izračuna parametrov geometričnega modela je zaradi manj natančnih točk daljše kot v primeru posnetkov RapidEye. Ortorektifikacija je računsko najbolj zahteven del obdelave in zajema preostalih 65 % časa. Samodejni postopek ortorektifikacije je del procesne verige STORM, ki je namenjena samodejni izdelavi izdelkov različnih ravni in kompleksnosti iz podatkov daljinskega zaznavanja. Veriga bo poleg posnetkov komercialnih satelitov podpirala tudi posnetke slovenskega satelita, katerega posnetki in izdelki bodo na voljo v skoraj realnem času. Pomembno vlogo pri obdelavi bo imel tudi algoritem za ortorektifikacijo, ki je že v trenutni različici dovolj zanesljiv in robusten. Kljub temu se lahko v prihodnje izdelan algoritem še dodatno izboljša in se mu doda druge funkcionalnosti. Možne izboljšave in nadgradnje so: - razširitev postopka za samodejno podporo drugim vrstičnim in polnoslikovnim snemalnim sistemom, - podpora tudi za posnetke z ločljivostjo 0,5 m oz. posnetke snemane v pankromatskem načinu, - podpora univerzalno uporabljenim projekcijam (npr. UTM) oz. drugim specifičnim (državnim) projekcijam in koordinatnim sistemom (npr. ECEF), - izboljšava algoritma za določanje oslonilnih točk na zelo visokoločljivih posnetkih, predvsem glede njihove razporeditve, - dodatna pospešitev celotnega postopka s pametnim programiranjem in uporabo porazdeljene obdelave (angl. distributed computing), - robustna in hitra implementacija direktne metode za ortorektifikacijo z lidarskim modelom reliefa ali modelom površja, 126 - izdelava geometričnega modela za obdelavo posnetkov geometrično že obdelanih stopenj (npr. ravni 1B), - razvoj in izdelava empiričnega geometričnega modela, ki deluje z racionalnimi polinomskimi koeficienti RPC. Poleg omenjenih izboljšav in nadgradenj je še veliko več možnosti razširitve postopka, ki jih ponujajo satelitski stereoposnetki (blokovna obdelava) in pa letalski posnetki oziroma posnetki brezpilotnih letalnikov (angl. unmanned aerial vehicles). Izdelana samodejna ortorektifikacija omogoča hitro pripravo satelitskih ortofotov, kar je pomembno tudi iz družbeno-ekonomskega vidika. V številnih aplikacijah, na primer pri opazovanju naravnih nesreč, nadzorovanju aktivnosti na morju, spremljanju vegetacije in podobno, je hitro kartiranje nujnega pomena. Pravočasna dostava ortorektificiranih satelitskih posnetkov lahko v primeru naravnih nesreč s hitrim posredovanjem prepreči večjo gmotno škodo ter celo rešuje človeška življenja. Hitro kartiranje je zelo pomembno tudi pri spremljanju kmetijskih zemljišč in različnih (legalnih in nelegalnih) posegov v okolje, kjer je hitra odzivnost nujna za preprečitev nastanka nepopravljive škode. Za pospešitev celotnega postopka je torej nujno potrebna visoka stopnja avtomatizacije geometrične obdelave posnetkov, ki jo lahko zagotovimo z izdelanim postopkom ortorektifikacije. 127 8 VIRI IN LITERATURA Airbus Defence and Space. 2017. Pléiades User Guide. http://www.intelligence-airbusds.com/en/4572-pleiades-technical- documents (3. 10. 2017). Allen, C. W. 1973. Astrophysical Quantities. London: The Athlone Press. http://archive.org/details/AstrophysicalQuantities. Bay, H., A. Ess, T. Tuytelaars in L. Van Gool. 2008. Speeded-Up Robust Features (SURF). Computer Vision and Image Understanding 110 (3): 346–359. Bigas, M., E. Cabruja, J. Forest in J. Salvi. 2006. Review of CMOS image sensors. Microelectronics Journal 37 (5): 433–451. Borčić, B. 1976. Gauss-Krügerova projekcija meridijanskih zona. Zagreb: Sveučilište, Geodetski fakultet. Brown, L. G. 1992. A survey of image registration techniques. ACM computing surveys (CSUR) 24 (4): 325–376. CEOS. 2017. The Committee on Earth Observation Satellite’s Earth Observation Handbook. http://www.eohandbook.com/ (3. 10. 2017). Chen, L.-C. in L.-H. Lee. 1993. Rigorous Generation of Digital Orthophotos from SPOT Images. Photogrammetric Engineering & Remote Sensing 59 (5): 655–661. Devaraj, C. in C. A. Shah. 2014. Automated geometric correction of multispectral images from High Resolution CCD Camera (HRCC) on- board CBERS-2 and CBERS-2B. ISPRS Journal of Photogrammetry and Remote Sensing 89: 13–24. DigitalGlobe. 2017a. Spletna stran podjetja DigitalGlobe. https://www.digitalglobe.com (3. 10. 2017). ———. 2017b. Accuracy of WorldView products. https://dg-cms-uploads- production.s3.amazonaws.com/uploads/document/file/38/DG_AC CURACY_WP_V3.pdf (3. 10. 2017). 128 Ebner, H., W. Konus in T. Ohlhof. 1992. A simulation study on point determination for the MOMS-02/D2 space project using an extended functional model. V Proceedings of the 17th ISPRS Congress, Commission IV, 29:458–464. International Archives of Photogrammetry and Remote Sensing. Washington, ZDA. e-GEOS. 2017. Price list. http://www.e-geos.it/products/pdf/prices.pdf (3. 10. 2017). El-Manadili, Y. in K. Novak. 1996. Precision rectification of SPOT imagery using the direct linear transformation model. Photogrammetric engineering and remote sensing 62 (1): 67–72. eoPortal. 2017a. RapidEye - eoPortal Directory - Satellite Missions. https://directory.eoportal.org/web/eoportal/satellite-missions/r/ rapideye (3. 10. 2017). ———. 2017b. WorldView-2 - eoPortal Directory - Satellite Missions. https://directory.eoportal.org/web/eoportal/satellite-missions/v-w- x-y-z/worldview-2 (3. 10. 2017). Eugenio, F. in F. Marques. 2003. Automatic satellite image georeferencing using a contour-matching approach. IEEE Transactions on Geoscience and Remote Sensing 41 (12): 2869–2880. Fischler, M. A. in R. C. Bolles. 1981. Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography. Communications of the ACM 24 (6): 381– 395. Fonseca, L. M. in B. S. Manjunath. 1996. Registration techniques for multisensor remotely sensed imagery. Photogrammetric Engineering & Remote Sensing 62 (9): 1049–1056. Gianinetto, M. in M. Scaioni. 2008. Automated geometric correction of high- resolution pushbroom satellite data. Photogrammetric Engineering & Remote Sensing 74 (1): 107–116. 129 Grocott, S. C. O., N. G. Orr, F. M. Pranajaya, R. E. Zee, T. Rodič, D. Matko, K. Oštir, M. Peljhan, A. Urbas, H. Frölich, S. Blažič in A. Marsetič. 2013. The NEMO-HD spacecraft with breadboard instrument test results. V Small satellites for Earth observation: digest of the 9th International Symposium of the International Academy of Astronautics, ur. R. Sandau, H. P. Röser in A. Valenzuela, 97–100. Berlin, Nemčija: Berlin: Wissenschaft und Technik. Grodecki, J. in G. Dial. 2003. Block Adjustment of High-Resolution Satellite Images Described by Rational Polynomials. Photogrammetric Engineering & Remote Sensing 69 (1): 59–68. Gugan, D. J. in I. Dowman. 1988. Topographic mapping from SPOT imagery. Photogrammetric Engineering and Remote Sensing 54 (10): 1409– 1414. GURS. 2017a. Centralna evidenca prostorskih podatkov. http://prostor3.gov.si/cepp/ (4. 10. 2017). ———. 2017b. E-prostor: Topografski podatki merila 1 : 5.000 (DTK 5). http://eprostor.agenda.si/si/zbirke_prostorskih_podatkov/topograf ski_in_kartografski_podatki/topografski_podatki_in_karte/topogra fski_podatki_merila_1_5000_dtk_5/ (4. 10. 2017). ———. 2017c. Zbirni kataster gospodarske javne infrastrukture – izdaja podatkov. http://www.e-prostor.gov.si/fileadmin/struktura/GJI_ izdaja_sifrant _in_struktura_2.pdf (4. 10. 2017). Gutman, G. in A. Ignatov. 1997. Towards a common language in satellite data management: a new processing level nomenclature. V Proceedings of the International Geoscience and Remote Sensing Symposium. IGARSS ’97, 3:1252–1254. Singapur. Harris Geospatial. 2017. Documentation Center. http://www.harrisgeospatial.com/docs/LINFIT.html (3. 10. 2017). Huber, P. J. 1964. Robust Estimation of a Location Parameter. The Annals of Mathematical Statistics 35 (1): 73–101. 130 Inglada, J. in A. Giros. 2004. On the possibility of automatic multisensor image registration. IEEE Transactions on Geoscience and Remote Sensing 42 (10): 2104–2120. Jovanovic, V., S. Lewicki, M. Smyth, J. Zong in R. Korechoff. 1999. Level 1 Georectification and Registration Algorithm Theoretical Basis. 20060034442. Pasadena, ZDA: Jet Propulsion Lab., California Inst. of Tech. Joyce, K. E., S. E. Belliss, S. V. Samsonov, S. J. McNeill in P. J. Glassey. 2009. A review of the status of satellite remote sensing and image processing techniques for mapping natural hazards and disasters. Progress in Physical Geography 33 (2): 183–207. Kim, T., D. Shin in Y.-R. Lee. 2001. Development of a Robust Algorithm for Transformation of a 3D Object Point onto a 2D Image Point for Linear Pushbroom Imagery. Photogrammetric Engineering and Remote Sensing 67 (4): 449–452. Klein, H. in W. Förstner. 1984. Realization of automatic error detection in the block adjustment program PAT-M43 using robust estimators. V Proceedings of the 15th ISPRS Congress, Commission III, 25:234–245. International Archives of Photogrammetry and Remote Sensing. Rio de Janeiro, Brazilija. Konecny, G. 1979. Methods and possibilities for digital differential rectification. Photogrammetric Engineering and Remote Sensing 45: 727–734. Konecny, G., P. Lohmann, H. Engel in E. Kruck. 1987. Evaluation of SPOT imagery on analytical photogrammetric instruments. Photogrammetric Engineering and Remote Sensing 53 (9): 1223– 1230. Kratky, V. 1989. Rigorous photogrammetric processing of SPOT images at CCM Canada. ISPRS Journal of Photogrammetry and Remote Sensing 44 (2): 53–71. Kraus, K. 1993. Photogrammetry, Vol. 1: Fundamentals and Standard Processes. Bonn: Dümmlers. 131 ———. 2007. Photogrammetry: Geometry from Images and Laser Scans. Prev. I. A. Harley in S. Kyle. 2 edition. Berlin: Walter de Gruyter. Kruck, E. 2001. Combined IMU and sensor calibration with BINGO-F. V Proceedings OEEPE Workshop on Integrated Sensor Orientation, Hannover, Germany, 1–3. Kubik, K., D. Merchant in T. Schenk. 1987. Robust estimation in photogrammetry. Photogrammetric Engineering and Remote Sensing 53 (2): 167–169. Kuipers, J. B. 2002. Quaternions and Rotation Sequences: A Primer with Applications to Orbits, Aerospace and Virtual Reality. Princeton, N.J: Princeton University Press. Le Moigne, J., W. J. Campbell in R. F. Cromp. 2002. An automated parallel image registration technique based on the correlation of wavelet features. IEEE Transactions on Geoscience and Remote Sensing 40 (8): 1849–1864. Leprince, S., S. Barbot, F. Ayoub in J. P. Avouac. 2007. Automatic and Precise Orthorectification, Coregistration, and Subpixel Correlation of Satellite Images, Application to Ground Deformation Measurements. IEEE Transactions on Geoscience and Remote Sensing 45 (6): 1529–1558. Liu, C.-C. in P.-L. Chen. 2009. Automatic extraction of ground control regions and orthorectification of remote sensing imagery. Optics Express 17 (10): 7970–7984. Lowe, D. G. 1999. Object Recognition from Local Scale-Invariant Features. V Proceedings of the International Conference on Computer Vision- Volume 2, 1150–1157. ICCV ’99. Washington, DC, USA: IEEE Computer Society. Marsetič, A. in K. Oštir. 2007. Uporaba satelitskih posnetkov SPOT za izdelavo ortopodob. Geodetski vestnik 51 (1): 69–84. Michalis, P. in I. Dowman. 2008. A Generic Model for Along-track Stereo Sensors Using Rigorous Orbit Mechanics. Photogrammetric Engineering & Remote Sensing 74 (3): 303–309. 132 Mikhail, E. M., J. S. Bethel in J. C. McGlone. 2001. Introduction to Modern Photogrammetry. 1 edition. New York : Chichester: Wiley. Netanyahu, N. S., J. Le Moigne in J. G. Masek. 2004. Georegistration of Landsat data via robust matching of multiresolution features. IEEE Transactions on Geoscience and Remote Sensing 42 (7): 1586–1600. Noerdlinger, P. D. 1999. Atmospheric refraction effects in Earth remote sensing. ISPRS Journal of Photogrammetry and Remote Sensing 54 (5): 360–373. NSSDC. 2017. The NASA Master Directory. https://nssdc.gsfc.nasa.gov/nmc/SpacecraftQuery.jsp (5. 10. 2017). Okamoto, A., C. Fraser, S. Hattori, H. Hasegawa in T. Ono. 1998. An alternative approach to the triangulation of spot imagery. V Proceedings of the ISPRS Commission IV Symposium, 32:457–462. International Archives of Photogrammetry and Remote Sensing. Stuttgart, Nemčija. O’Neill, M. in I. Dowman. 1988. The Generation of Epipolar Synthetic Stereo Mates for SPOT using a DEM. V Proceedings of the 16th ISPRS Congress, Commission III, 27:587–598. International Archives of Photogrammetry and Remote Sensing. Kyoto, Japonska. OpenStreetMap. 2017. Spletna stran organizacije OpenStreetMap. https://www.openstreetmap.org/ (5. 10. 2017). Oštir, K. 2006. Daljinsko zaznavanje. Ljubljana: Založba ZRC. Oštir, K., A. Marsetič, P. Pehani, M. Perše, K. Zakšek, J. Zaletelj in T. Rodič. 2014. Procesna veriga za samodejno obdelavo optičnih satelitskih posnetkov v skoraj realnem času. V GIS v Sloveniji 12: Digitalni prostor, ur. R. Ciglič, D. Perko in M. Zorn, 207–218. Ljubljana: Založba ZRC. Planet Labs. 2017. RapidEye Imagery Product Specifications. https://www.planet.com/products/satellite-imagery/files/160625- RapidEye%20Image-Product-Specifications.pdf (3. 10. 2017). Poli, D. 2005. Modelling of spaceborne linear array sensors. Doctoral Thesis, Zurich: ETH Zurich. 133 Press, W. H., S. A. Teukolsky, W. T. Vetterling in B. P. Flannery. 2007. Numerical Recipes 3rd Edition: The Art of Scientific Computing. 3rd edition. Cambridge ; New York: Cambridge University Press. Sandau, R. 2004. High resolution mapping with small satellites. V Proceedings of the 20th ISPRS Congress, Commission I, 35:108–113. International Archives of Photogrammetry and Remote Sensing. Istanbul, Turčija: ISPRS. Schenk, T. 1999. Digital Photogrammetry, Volume I: Background, Fundamentals, Automatic Orientation Procedures. 1st edition. Laurelville, OH: TerraScience. Sheng, Y. 2004. Comparative evaluation of iterative and non-iterative methods to ground coordinate determination from single aerial images. Computers & Geosciences 30 (3): 267–279. Sheng, Y. in D. E. Alsdorf. 2005. Automated georeferencing and orthorectification of Amazon basin-wide SAR mosaics using SRTM DEM data. IEEE Transactions on Geoscience and Remote Sensing 43 (8): 1929–1940. Slonecker, E. T., B. Johnson in J. McMahon. 2009. Automated imagery orthorectification pilot. Journal of Applied Remote Sensing 3 (1): 33552–33568. Smith, G. S. 1995. Digital Orthophotography and GIS. Esri International User Conference Proceedings. Redlands, ZDA. Tao, C. in Y. Hu. 2001. A Comprehensive study of the rational function model for photogrammetric processing. Photogrammetric Engineering and Remote Sensing 67 (12): 1347–1357. Toutin, T. 2003. Error Tracking in Ikonos Geometric Processing Using a 3D Parametric Model. Photogrammetric Engineering & Remote Sensing 69 (1): 43–51. ———. 2004. Review article: Geometric processing of remote sensing images: models, algorithms and methods. International Journal of Remote Sensing 25 (10): 1893–1924. 134 Veljanovski, T., P. Pehani, P. Lamovec in K. Oštir. 2012. Uporabnost podatkov satelitskega in letalskega daljinskega zaznavanja za opazovanje in kartiranje vodnih površin. Geodetski vestnik 56 (4): 786–801. Weser, T., F. Rottensteiner, J. Willneff, J. Poon in C. S. Fraser. 2008. Development and testing of a generic sensor model for pushbroom satellite imagery. The Photogrammetric Record 23 (123): 255–274. Westin, T. 1990. Precision rectification of SPOT imagery. Photogrammetric Engineering & Remote Sensing 56 (2): 247–253. Xiong, Z. in Y. Zhang. 2009. A Novel Interest-Point-Matching Algorithm for High-Resolution Satellite Images. IEEE Transactions on Geoscience and Remote Sensing 47 (12): 4189–4200. Zaletelj, J., U. Burnik in J. F. Tasic. 2013. Registration of satellite images based on road network map. V Proceedings of the 8th International Symposium on Image and Signal Processing and Analysis (ISPA), 49– 53. Trst, Italija: IEEE. Zhang, C. in C. S. Fraser. 2007. Automated registration of high-resolution satellite images. The Photogrammetric Record 22 (117): 75–87. Zhou, F., Y. Cui, Y. Wang, L. Liu in H. Gao. 2013. Accurate and robust estimation of camera parameters using RANSAC. Optics and Lasers in Engineering 51 (3): 197–212. Zhou, G. 2009. Near Real-Time Orthorectification and Mosaic of Small UAV Video Flow for Time-Critical Event Response. IEEE Transactions on Geoscience and Remote Sensing 47 (3): 739–747. Zitová, B. in J. Flusser. 2003. Image registration methods: a survey. Image and Vision Computing 21 (11): 977–1000. 135 Seznam slik Slika 1: Samodejni postopek ortorektifikacije. Za izdelavo satelitskega ortofota je potrebnih več korakov. ...................................................................................................... 13 Slika 2: Samodejni postopek za ortorektifikacijo je sestavljen iz štirih osnovnih modulov. ................................................................................................................................. 27 Slika 3: Uporabljeni surovi posnetki RapidEye (spodaj) in njihov položaj (zgoraj). .... 31 Slika 4: Uporabljeni surovi posnetki WorldView-2 (spodaj) in njihov položaj (zgoraj). .................................................................................................................................... 34 Slika 5: Potek modula za samodejno določanje oslonilnih točk. ....................................... 39 Slika 6: Položaj samodejno določene oslonilne točke. Točka se nahaja na križišču izseka zaznanih cest (siva), v ozadju pa je prikazan izsek podatkov oddaljenosti od cest. ........................................................................................................... 40 Slika 7: Zajem posnetka z vrstičnim sistemom. ........................................................................ 58 Slika 8: Zajem posnetka s polnoslikovnim sistemom. ........................................................... 60 Slika 9: Geometrični prikaz atmosferske refrakcije. ................................................................ 62 Slika 10: Modeliranje tirnice z dvema segmentoma. V praksi je prehod med segmentoma veliko bolj gladek, kot je prikazano na sliki, kar zagotovimo z enačenjem vrednosti funkcij segmentov ter njihovih prvih in drugih odvodov. ................................................................................................................................. 69 Slika 11: Shematizacija poteka izračunov parametrov geometričnega modela. ........... 89 Slika 12: RMSE na kontrolnih točkah za vsak točkovni niz za posnetke RapidEye. ........ 91 Slika 13: RMSE na kontrolnih točkah za vsak točkovni niz posnetkov WorldView-2. ... 93 Slika 14: Obnašanje preizkušenih metod robustne ocene za posnetke RapidEye. ....... 95 Slika 15: Lastnosti preizkušenih metod robustne ocene za posnetke WorldView-2. .. 97 Slika 16: Rezultati izravnave (RMSE na kontrolnih točkah) po izravnavi z algoritmom RANSAC (črtkana črta) in na koncu postopka (polna črta) za večji in manjši niz točk posnetkov RapidEye. Iz slike je razvidno, da v splošnem robustna ocena izboljša natančnost rezultatov izravnave. .............. 99 Slika 17: Porazdelitev in odstopanja oslonilnih in kontrolnih točk. Na zgornji sliki je prikazan posnetek Koper (RapidEye). Oslonilne točke so prikazane z 136 rumeno barvo, kontrolne pa z zeleno. Na spodnji sliki sta neodstranjeni točki z dodanimi napakami označeni z rdečima elipsama. ............................... 101 Slika 18: Porazdelitev in odstopanja oslonilnih in kontrolnih točk. Na zgornji sliki (posnetek Bled – WorldView-2) so oslonilne točke prikazane z rumeno barvo, kontrolne pa z zeleno. Na spodnji sliki so neodstranjene točke z dodanimi napakami označene z rdečimi elipsami................................................ 102 Slika 19: Rezultati izravnave (RMSE kontrolnih točk) po izravnavi z algoritmom RANSAC (črtkana črta) in na koncu postopka (polna črta) za večji in manjši niz točk posnetkov WorldView-2. Iz slike je razvidno, da v splošnem robustna ocena izboljša natančnost rezultatov izravnave. ........... 104 Slika 20: Iterativna fotogrametrična metoda (IP). ................................................................... 107 Slika 21: Iterativna metoda sledenja žarku (IRT). ..................................................................... 109 Slika 22: Ponazoritev preslikave pikslov pri izdelavi ortofota z direktnimi metodami. Zaradi preslikave pikslov (T6 in T8) v sosednje celice sta dve celici v ortofotu prazni (P1 in P2). .................................................................................. 110 Slika 23: Direktna metoda ortorektifikacije za vrstične skenerje. ...................................... 110 Slika 24: Ortorektifikacija polnoslikovnega posnetka z indirektno metodo. ................ 112 Slika 25: Indirektna metoda ortorektifikacije za vrstični skener. Pri ortorektifikaciji je treba upoštevati geometrični model, v katerem je tirnica razdeljena na dva segmenta (S1 in S2). ................................................................................................... 113 Slika 26: Določanje velikosti ortofota (matrike) glede na surovi posnetek. .................. 116 Slika 27: RMSE na 30 ročno določenih kontrolnih točkah za posnetke RapidEye....... 118 Slika 28: Ortorektificiran posnetek Koper (RapidEye). ........................................................... 119 Slika 29: Ortorektificiran posnetek Bled – WorldView-2. ...................................................... 120 Slika 30: RMSE na 30 ročno določenih kontrolnih točkah za posnetke WorldView-2. ....................................................................................................................... 121 137 Seznam preglednic Preglednica 1: Lastnosti satelitov RapidEye (eoPortal 2017a). ............................................ 30 Preglednica 2: Lastnosti uporabljenih posnetkov RapidEye. ............................................... 30 Preglednica 3: Lastnosti satelita WorldView-2 (eoPortal 2017b). ....................................... 33 Preglednica 4: Lastnosti uporabljenih posnetkov WorldView-2. ....................................... 33 Preglednica 5: Točnost ročno določenih oslonilnih točk glede na koordinate izračunane iz racionalnih polinomskih koeficientov RPC za posnetke RapidEye. ............................................................................................... 51 Preglednica 6: Točnost samodejno določenih oslonilnih točk glede na koordinate izračunane iz racionalnih polinomskih koeficientov RPC za posnetke RapidEye. ............................................................................................... 52 Preglednica 7: Točnost ročno določenih oslonilnih točk glede na koordinate izračunane iz racionalnih polinomskih koeficientov RPC za posnetke WorldView-2. ........................................................................................ 53 Preglednica 8: Točnost samodejno določenih oslonilnih točk glede na koordinate izračunane iz racionalnih polinomskih koeficientov RPC za posnetke WorldView-2. ........................................................................................ 54 Preglednica 9: Vrednosti atmosferske refrakcije na morski gladini (Noerdlinger 1999). ........................................................................................................................... 63 Preglednica 10: Točnost rezultatov izravnave, ki so bile dobljene z različnim številom oslonilnih in kontrolnih točk za posnetke RapidEye. .............. 91 Preglednica 11: Točnost rezultatov izravnave, ki so bile dobljene z različnim številom oslonilnih in kontrolnih točk za posnetke WorldView-2. ....... 93 Preglednica 12: Primerjava metod za robustno oceno na posnetkih RapidEye. ............ 95 Preglednica 13: Primerjava metod za robustno oceno na posnetkih WorldView-2. ..... 96 Preglednica 14: Rezultati izravnave s številom prisotnih in izločenih oslonilnih točk z dodanimi napakami, številom vseh in uporabljenih oslonilnih točk ter RMSE po izravnavi z algoritmom RANSAC in na koncu postopka (robustna ocena) za posnetke RapidEye. ................................... 98 Preglednica 15: Velikost odstopanj oslonilnih točk z dodanimi napakami, ki so ostale po izravnavi z algoritmom RANSAC za večji in manjši niz točk posnetkov RapidEye. ................................................................................... 99 138 Preglednica 16: Rezultati izravnave s številom prisotnih in izločenih oslonilnih točk z dodanimi napakami, številom vseh in uporabljenih oslonilnih točk ter RMSE po izravnavi z algoritmom RANSAC in na koncu postopka (robustna ocena) za posnetke WorldView-2. Zaradi slabega delovanja geometričnega modela poskusi s posnetkom Ljubljana niso bili izvedeni. .............................................................................. 103 Preglednica 17: Velikost odstopanj oslonilnih točk z dodanimi napakami, ki so ostale po izravnavi z algoritmom RANSAC za večji in manjši niz točk posnetkov WorldView-2. Zaradi slabega delovanja geometričnega modela poskusi s posnetkom Ljubljana niso bili izvedeni. .................................................................................................................. 103 Preglednica 18: Točnost izdelanih ortofotov RapidEye glede na državni ortofoto, ki je bil uporabljen kot referenca. ...................................................................... 118 Preglednica 19: Točnost izdelanih ortofotov WorldView-2 glede na državni ortofoto, ki je bil uporabljen kot referenca. ............................................... 121 139 Priloga PRILOGA A: PORAZDELITEV IN ODSTOPANJA OSLONILNIH IN KONTROLNIH TOČK V prilogi se nahajajo posnetki z označenimi položaji točk in odstopanji za vse preizkušene posnetke z nizom s ~50 točkami. Na zgornjem delu slik so oslonilne točke prikazane z rumeno barvo, kontrolne točke pa z zeleno. Oslonilne točke, ki jih RANSAC ni izločil, imajo rdečo obrobo, začetne RANSAC točke pa modro. Na spodnjem delu slik so prikazana odstopanja (položajne napake) vseh točk. 140 RapidEye Radgona 141 RapidEye Koper 142 RapidEye Bohinj 143 WorldView-2 Jesenice 144 WorldView-2 Bled 145 WorldView-2 Ljubljana 146 Urednika zbirke: Nataša Gregorič Bon in Žiga Kokalj, ZRC SAZU Prostor Zbirka je namenjena objavi krajših tematsko zaokroženih znanstvenih kraj raziskav s področja sodobnega merjenja prostora, ki temeljijo na geografskih informacijskih sistemih in daljinskem zaznavanju, kot tudi na družbenih in čas 15 kulturnih konstrukcijah prostora in časa: kako ljudje v različnih obdobjih in pokrajinah mislimo prostor in čas, kako ju živimo, čutimo, ustvarjamo, spreminjamo in uporabljamo. GEOMETRIČNI POPRAVKI OPTIČNIH SATELITSKIH POSNETKOV Aleš Marsetič, Krištof Oštir in Mojca Kosmatin Fras Dr. Aleš Marsetič je zaposlen kot znanstveni sodelavec na Inštitutu za antropološke in prostorske študije ZRC SAZU. V svojem raziskovalnem delu se posveča predvsem daljinskemu zaznavanju in fotogrametriji. Veliko se ukvarja z avtomatizacijo predobdelave podatkov in izdelavo produktov. Že vrsto let sodeluje pri projektiranju in izdelavi majhnih satelitov. Izvaja tudi terensko delo in obdelavo ter vizualizacijo podatkov na arheoloških projektih v Mehiki. Dr. Krištof Oštir je redni profesor na Univerzi v Ljubljani in znanstveni svetnik na Znanstvenoraziskovalnem centru SAZU. Ukvarja se z raziskavami in uporabo optičnih in radarskih satelitskih posnetkov ter razvojem tehnologije malih satelitov za opazovanje Zemlje. Dr. Mojca Kosmatin Fras je doktorirala leta 2003 iz geodetskih znanosti na Politehniki v Milanu, Italija. Od tega leta naprej je zaposlena kot docentka na Univerzi v Ljubljani, Fakulteti za gradbeništvo in geodezijo. Njeno raziskovalno delo obsega različne vidike zajema in ocene kakovosti topografskih podatkov, pridobljenih iz aeroposnetkov, visokoločljivih satelitskih posnetkov in lidarskih podatkov. S temo izdelave ortofota iz različnih virov se ukvarja že od začetka svoje službene poti pred več kot tridesetimi leti. Knjiga je prosto dostopna v e-obliki (pdf) na: http:\\zalozba.zrc-sazu.si http://zalozba.zrc-sazu.si/p/P15 ZRC Publishing