DMT3 Prostorová interpolace dat DjcjjííjJ/jj jr/jodeJy iará/ju (3) Prostorová interpolace dat Ing. Martin KLIMÁNEK, Ph.D. 411 Ústav geoinformačních technologií Lesnická a dřevařská fakulta, Mendelova zemědělská a lesnická univerzita v Brně L* *'LV'C/,? lil DMT3 Prostorová interpolace dat Interpolace - úvod • procedura odhadu neznámých hodnot ze známých hodnot v okolí • interpolace • extrapolace • úkolem interpolace je určit vhodnou funkci y=f(x), která v daných (tzv. uzlových) bodech nabývá známých hodnot interpolační funkce f(x) - shoda s původními hodnotami v uzlových bodech aproximující funkce f(x) - nahrazuje původní funkci f0(x) s jistou přesností Dodržují se známé hodnoty? • interpolační (exaktní) metody (deterministické) • aproximující metody (stochastické) DMT3 Prostorová interpolace dat Typy odhadů - odhad veličiny v bodě, kde nebyla zjištěna - gridování - odhad pro body pravidelné sítě Y 16, DMT3 Prostorová interpolace dat Typy odhadů blokový globální rozmístění měření - některé metody selhávají pro silně nepravidelné sítě měření - zvláště problém hodnot na profilech DMT3 Zdrojová data rozmístěni zdrojových dat Prostorová interpolace dat v ► • • • • •Ä pravidelné rozmístění náhodné rozmístění transekty (profily) # # • # • i # • • # • • • # • • • • # # # • • # — \ •• % " stratifikované náhodné shluky (skupiny) vrstevnice (izolinie) DMT3 Prostorová interpolace dat Metody interpolace Vážený průměr DMT3 Prostorová interpolace dat Metoda inverzních vzdáleností (ID, IDW) - princip Velikost příspěvku je přímo úměrná velikosti hodnoty a na druhé straně nepřímo úměrná vzdálenosti. r=dX2+dY limMy =M „M," je známá hodnota v i-tém místě, „r" vzdálenost i-tého místa od místa X a „k" je vhodná mocnina vzdálenosti (např. 1 nebo 2). DMT3 Prostorová interpolace dat Metoda inverzních vzdáleností (ID) - varianty: • s výběrem z kvadrantů nebo oktantů - rozdělení prohledávané oblasti do 4 nebo 8 sektorů • se směrovým vážením - rozdělení prohledávané oblasti do nepravidelných sektorů. Vstupní parametry: • úhlové rozpětí sektorů, váhy přiřazené jednotlivým sektorům = (zavádění anizotropie pole) • z každého sektoru je vybrán 1 bod (nebo „k" bodů) a na tuto množinu vybraných bodů je aplikována metoda ID Shepardova metoda (SH) - metoda inverzních vzdáleností + vyrovnání metodou nejmenších čtverců DMT3 Prostorová interpolace dat Metoda triangulace s lineární interpolací (TR) • poskytuje odhad neznámé hodnoty pomocí lineární závislosti. Lineárním útvarem je tedy rovina. Rovnice roviny obsahuje tři koeficienty - to znamená, že pro její určení jsou zapotřebí tři známé body: P=[x1,y1,z1], Q=[x2,y2,z2], R=[x3,y3,z3]. • tři body P, Q a R jsou vybrány tak, aby bod X ležel uvnitř trojúhelníka získaného jako průmět trojúhelníka PQR. Těmito body je určena rovina, jejíž koeficienty lze získat řešením soustavy rovnic. Hledaný odhad Zx je dán vztahem: DMT3 Prostorová interpolace dat Metoda triangulace s lineární interpolací dodržuje naměřené hodnoty • výhodná je pro modelování nespojitostí v poli (zlom, hrana) nevýhodná pro nespojitost vypočtených hodnot (kostrbatý průběh izolinií) - různé stupně vyhlazení 7M A \ ■ÜUI ■ \ > / aoi P\ \ r\ \ v \ \ a \ 0 ílÁt ?Sa??i 3 ?SSäi *J V3 DMT3 Prostorová interpolace dat Metoda triangulace s lineární interpolací interpolace izolinií DMT3 Prostorová interpolace dat Metoda triangulace s lineární interpolací • Skrývá jistý prvek náhodnosti. Uvažujme následující data: P=[1,1,1], 0=13,1,10], S=[3,1,10], R=[3,3,1]. Odhad proX=[2,2,?]. • Náhodnost spočívá ve volbě trojúhelníka: Volbou PRS, je odhadem hodnota 1. Volbou QRS, je odhadem hodnota 10. Problém je tedy ve velkém rozdílu dat mezi „sousedními" hodnotami. Bohužel právě takový charakter má mnoho geoprostorových dat, a proto zvláště na taková data není vhodné metodu aplikovat. S dobrými výsledky ji lze použít v případech, kdy hodnoty nemají příliš velký rozptyl. DMT3 Thiessenovy (Dirichlet, Voronoi) polygc Představují přesnou metodu interpolace že neznámé hodnoty bodů odpovídají bodů. Zahrnuje šíření teritoria sdružení tak dlouho, až se narazí na obdol sousedního bodu. Je-li rozmístění bodů * • • / • \ Prostorová interpolace dat >ny (TP) , která vychází z předpokladu, hodnotě nejbližších známých 3ho s bodem, které pokračuje bně zpracovávané teritorium nepravidelné, výsledkem bude / / Thiessen polygons (thick lines) Delaunay triangulation (thin lines) DMT3 Prostorová interpolace dat Thiessenovy (Dirichlet, Voronoi) polygony (TP) • Toblerova pyknofylaktická metoda: l jf(x,y)dxdy = H „H" jsou hodnoty v regionu „R" pro všechna „i" Administrativní celky Pyknofylaktická interpolace DMT3 Prostorová interpolace dat Metoda minimální křivosti (MC) • spline-funkce - vychází z interpolace pomocí (nejčastěji kubických) funkcí. „Spojuje" tedy dvojice daných bodů segmenty kubické křivky (ten je dán čtyřmi body). Z prvních čtyř bodů se spočte kubická křivka a první dva body se spojí jejím segmentem. Pak se z druhého až pátého bodu spočte kubická křivka a druhé dva body se spojí jejím segmentem, atd. • polynomické funkce, na hranách spojité (spojitost interpolující funkce a stanoveného počtu jejich derivací) • povrch je interpolován po částech hladké povrchy • míra aproximace je dána stanovením tolerance a počtem iterací DMT3 Prostorová interpolace dat Metoda radiálních funkcí (RF) • hladký povrch • dodržuje naměřené hodnot; • m u Iti kvádri ková metoda D|(x,y) je radiální funkce vzdálenosti dj(x,y) dj(x,y) je relativní, anizotropní vzdálenost mezi místem měření (x, .y,) a místem odhadu (x, y) R2 je vyhlazovací faktor V každém místě jsou stanoveny optimální váhy řešením soustavy lineárních rovnic. DMT3 Prostorová interpolace dat Approximation Based On Smoothing (SurGe, M. Dressler, http://mujweb.cz/www/surge/) • pridelení nejbližší známe hodnoty každému nodu mrizky, (per partes konstantní interpolace, Thiessenovy polygony) • průměrování hodnot mřížky (vyrovnávání rozdílu): vyhlazování hodnot mřížky (s ohledem na požadované parametry) několik iterací výsledný povrch se vypočte bihnearni interpolaci z rohových nodu: DMT3 Prostorová interpolace dat mm* DMT3 Prostorová interpolace dat DMT 3 Prostorová interpolace dat Fourierova analýza (FA) • obecný interpolační postup, který se používá v geostatistických metodách pro odhad průběhu povrchu pomocí série funkcí sinus a kosinus • nejvhodnější pro datové soubory, které vykazují periodicitu - analýza signálů, modelování klimatických změn, pohyb (vln v oceánu, sněhových závějí, pískových dun apod.) a zpracování obrazu (odstranění šumu, korekce); redukce variability (kontrastů) v DMT pro snadnější identifikaci tvarů terénu (hřbety, údolí) frekvence (spektrum) = amplituda k = koeficient harmonické funkce a, b = koeficienty Fourierovy transformace Nyquistova frekvence - nejvyšší frekvence = nejkratší vlna = 2-pixel DMT3 Prostorová interpolace dat Krigování (KR) - úvod • „teorie regionalizovaných proměnných", která byla vytvořena G. Matheronem a D. G. Krigem (metoda interpolace dat získávaných při vyhledávání a průzkumu ložisek nerostných surovin) • vychází ze zjištění, že prostorová variabilita řady geoprostorových prvků je příliš nepravidelná než, aby mohla být modelována pomocí vyhlazovacích matematických funkcí - optimalizovaná soustava vah v každém kroku s požadavkem minimálního rozptylu •základ spočívá v nalezení průměrné hodnoty změn v závislosti na změně vzdálenosti mezi měřenými body - vyhodnocování strukturálních funkcí • optimalizuje prostorovou interpolaci tak, že rozděluje prostorovou variabilitu do tří komponent: • deterministickou variabilitu (různé trendy) • prostorově autokorelovanou variabilitu ■ nekorelovaný šum • charakter prostorově korelované variability lze znázornit jako semivariogram, který poskytuje informaci pro optimalizaci interpolačních vah a prohledávacích poloměrů DMT3 Prostorová interpolace dat Krigování (KR) - úvod Výchozí předpoklady pro strukturální analýzu a krigování: normální distribuce • transformace nebo použití nelineárních technik »růměrná hodnota konstantní • trend, univerzální krigování prostorová autokorelace konstantní DMT3 Prostorová interpolace dat Krigování (KR) - experimentální semivariogram Semivariogram je základní geostatistický nástroj pro vizualizaci, modelování a využití prostorové autokorelace regionalizované proměnné. Semivariance je míra variance. kde y(h) je semivariance proměnné z pro vzdálenost h Základní schéma výpočtu pro vzdálenost h sestává z následujících kroků: • v úvahu se vezmou všechny takové páry zx a zjf jejichž vzdálenost padne do třídy pro h • vypočtou se rozdíly hodnot těchto párů • sečtou se kvadráty těchto rozdílů • součet se vydělí dvojnásobkem počtu párů Tak se pro všechny hodnoty h získá řada hodnot nazývaných experimentální semivariance. DMT3 Prostorová interpolace dat Krigování (KR) - teoretický semivariogram • Funkční vztah, který pokud možno dobře sleduje (modeluje) experimentální semivariogram, se nazývá teoretický semivariogram 1 Chování semivariogramu je možno (víceméně) intuitivně popsat takto: • velmi blízká data mají velmi malou odchylku data ve větších vzdálenostech mají větší odchylky, avšak odchylky pro velmi vzdálená data a velmi velice vzdálená data se už příliš neliší • od jisté vzdálenosti už vzájemné odchylky nerostou (např. také proto, že vzdálenost překračuje rozměry zkoumané plochy nebo tělesa) DMT3 Prostorová interpolace dat Krigování (KR) - strukturní analýza • proces hledání teoretického semivariogramu pro daný experimentální semivariogram je nazýván strukturní analýzou. Model nalezený pro danou množinu dat závisí jak na experimentálních, tak teoretických předpokladech. Vlastnosti, které prakticky vedou k určení konkrétního teoretického modelu, jsou: • přítomnost nebo absence --------------------------------------------------- „ploché části" semivariogramu, to znamená, že při zvětšující se vzdálenosti se hodnoty variance nemění - existuje tzv. práh (sill); je dán konstantou • vzdálenost, ve které semi-variance dosáhne prahové hodnoty - tzv. dosah (range); je dán konstantou „a" • chování v počátku (tj. semi-variance mezi velmi blízkými body) i k dosah (a) (range) 3 \ Ö u > B f práh (C / (si») v 1 ' ) \ model semivariogramu cx> i L zbytkový rozptyl (nugget) ' ' vzdálenost (h) DMT3 Prostorová interpolace dat Krigování (KR) - strukturní analýza • dosah (range) je mírou korelace uvnitř množiny dat; „dlouhý" dosah indikuje vysokou korelaci, „krátký" korelaci nízkou hodnota prahu (sill) je rovna celkovému rozptylu • velmi často nenabývají experimentální semivariogra-my v počátku nulové hodnoty; protínají osu y v nenulové hodnotě, která je nazývána zbytkový rozptyl (nugget effect). To může ukazovat na rozptyl menší než je „vzorkovací" vzdálenost nebo na malou přesnost měření (např. jsou v datech obsaženy dva vzorky ze stejného místa, pokaždé s jinou hodnotou). i k dosah (a) (range) 3 \ Ö u > B f práh (C / (si») v 1 ' ) \ model semivariogramu cx> i L zbytkový rozptyl (nugget) ' ' vzdálenost (h) DMT3 Prostorová interpolace dat Krigování (KR) - modely semivariogramů Sférický model Exponenciální model Gaussův model • Složené modely - každý zdroj variability je popsán vlastní strukturní funkcí DMT3 Prostorová interpolace dat Krigování (KR) - modely semivariogramů KW C J b) .B--"^""^ Sill c1 range a nugget a) sférický b) exponenciální c) lineární en aaussovskv Lag(h)- Laq(h> Lag(h)- Lag(h)- DMT3 Prostorová interpolace dat Krigování (KR) - směrovost semivariancí • mnohé fenomény nevykazují prostorovou isotropii rozptylu - naopak se chovají anisotropně - rozsah influence je různý v různých směrech • pro data, vykazující anisotropii, se využívá postupu spočívajícího v konstrukci dílčích semivariogramů pro stanovené směrové tolerance V nejjednodušším případě rozdělíme všesměrové pole rovnoměrně: • například na směry sever ± 22.5°, jih ± 22.5°, východ ± 22.5°, západ ± 22.5°, severovýchod ± 22.5°, jihovýchod ± 22.5°, jihozápad + 22.5°, severozápad + 22.5° a zkonstruujeme osm dílčích směrových semivariogramů tak, že do každého z nich zahrneme pouze ty dvojice daných bodů, jejichž směrový vektor padne do intervalu směrů daného dílčího semivariogramů. DMT3 Prostorová interpolace dat Krigování (KR) - směrovost semivariancí • Při směrově symetrických datech stačí zkonstruovat pouze polovinu počtu semivariogramů. Po zjištění jejich dosahů se tyto dosahy zanesou do růžicového diagramu obsahujícího ty směry, pro něž byly dílčí semivariogramy sestaveny. • Anizotropická množina dat je charakterizována směrem maximální variance a směrem minimální variance. Tyto směry jsou směry hlavní a vedlejší poloosy tzv. elipsy anizotropie. Elipsa anizotropie je pak zjistitelná jako elipsa, která aproximuje dosahy vynesené do shora zmíněného růžicového diagramu. dosah semivariogramů pro úhel 180° ±22,5° osy anizotropie dosah semivariogramů pro úhel - 45° ± 22,5° DMT3 Prostorová interpolace dat Krigování (KR) - varianty Základní krigování: bodový odhad, (blokový odhad) • krigování proto stanoví odhad v místě X jako součet vážených známých hodnot •minimalizace rozptylů odhadů v místě X vede k soustavě (n+1) lineárních rovnic, přičemž Kr} jsou kovariance mezi i-tým a j-tým známým bodem vázané se semivariahcemi vztahem: • Variance tohoto odhadu je dána vztahem (kde X je Lagrangeuv multiplikátor): DMT3 Prostorová interpolace dat Krigování (KR) - základní krigování (příklad) • určit koeficienty \iv \i2, \i3, tyto koeficienty jsou výsledkem řešení soustavy lineárních rovnic •jako prvky rozšířené matice soustavy vystupují hodnoty semivariance. Protože jde o aproximované hodnoty, je prvním úkolem stanovení modelu (= teoretického semivariogramu) • model = lineární semivariogram bez zbytkového rozptylu se směrnicí přímky rovné 4,0 m2/km =^> y(h)=4h • Pro zjištění prvků rozšířené matice soustavy je tedy zapotřebí: • zjistit souřadnice jednotlivých bodů • zjistit vzájemné vzdálenosti mezi nimi zjistit hodnoty semivariancí pro zjištěné vzdálenosti DMT3 Prostorová interpolace dat Krigování (KR) - základní krigování (příklad) . Řešíme tedy soustavu: Bod Souřadnice X Souřadnice Y Hladina vody 1 3.0 4.0 120 2 6.3 3.4 103 3 2.0 1.3 142 X 3.0 3.0 ??? z - do [krn 1 2 3 X 1 0.000 3.354 2.879 1.000 2 0.000 4.735 3.324 3 0.000 1.972 X 0.000 Semivariance 1 2 3 X 1 0.000 13.416 11.517 4.000 2 0.000 19.142 13.297 3 0.000 7.339 X 0.000 --------------1| /'i A'2 Sh Á 1 ° 13.416 11.517 1 4.000 J13.416 0 19.142 1 13.297 li 1.517 19.142 0 1 7.889 1 1 1 0 1 Řešením je čtveřice [nv |i2, Ha, X] = [0,5954; 0,0975; 0,3071;-0.7298]. Odhad v bodě X je tedy: Zx = 0,5954 -120 + 0,0975 ■ 103 + 0,3071-142 = 125,1 [m] s variancí chyby odhadu: s2 = 0,5954.4 + 0,0975-13,416 + 0,3071-7,889-0,7298 1 = 5,3826 [m2] X= 125,1+2,3 m DMT3 Prostorová interpolace dat Krigování (KR) - varianty • simple kriging - jednoduché krigování - známá průměrná hodnota veličiny ^ """"I • cokriging - odhad proměnné na základě hodnot korelovaných veličin, vztahy popsané pomocí vzájemného semivariogramu ■ lognormální krigování • indikátorové krigování - transformace kvantitativních údajů na kvalitativní (ano/ne) - výsledkem lokálního odhadu je pravděpodobnost; (absolutní nebo relativní indikátory - problém konzistence) • probability kriging - kombinace indikátorů s původními údaji • soft kriging - kombinace hodnoty s dalšími údaji DMT3 Prostorová interpolace dat Krigování (KR) - varianty (ukázky) ■ a) Zdrojová data b) Třídy četnosti záplav c) Koncentrace Zn dle četnosti záplav d) Indikátorové krigování: pravděpodobnost překročení hodnoty 1000ppmZn e) Soft krigování: pravděpodobnost překročení hodnoty 500ppm s využitím pravděpodobné hodnoty obsahu Zn v jednotlivých třídách četnosti záplav f) Kokrigování In Zn a In vzdálenosti od řeky DMT3 Prostorová interpolace dat Podmíněná stochastická simulace (SS) 1. Vyberte náhodně dosud nezpracovaný bod (neznámé hodnoty) 2. Vypočítejte pro něj pomocí jednoduchého krigování odhad hodnoty a rozptyl odhadu s využitím známých (nebo již vypočítaných) hodnot. 3. Stanovte náhodnou hodnotu z distribuce pravděpodobnosti určené zjištěným průměrem (zjištěný odhad hodnoty) a rozptylem (zjištěný rozptyl odhadu). Výsledek představuje simulovanou hodnotu pro dané místo. Bod se stává místem se známou hodnotou a vstupuje do dalších výpočtů. 4. Opakujte první 3 kroky, až je pokryto celé území. 5. Opakujte první 4 kroky tolikrát, kolik simulací je potřebné provést pro dostatečně věrohodný odhad modelu (např. 100 krát). 6. Ze simulovaných hodnot vypočítejte průměrnou hodnotu a rozptyl (resp. směrodatnou odchylku) pro každé místo. DMT3 Podmíněná stochastická simulace (ukázky) Podmíněná simulace obsahu Zn v půdě po 100 iteracích: • vlevo očekávaná hodnota Zn vpravo chyba odhadu Podmíněná simulace obsahu Zn v půdě po 100 iteracích s využitím informace o pravděpodobné hodnotě obsahu Zn v jednotlivých třídách četnosti záplav: • vlevo očekávaná hodnota Z • vpravo chyba odhadu Prostorová interpolace dat DMT3 Volba metody a parametrů Prostorová interpolace dat • Průzkumová analýza dat • zkoumání distribuce hodnot, • statistická (normalita, extrémní hodnoty), • prostorová (existence trendu, anizotropie), • uspořádání sítě měření. • Transformace hodnot měření - normální distribuce (případně použití nelineárních technik) • Eliminace trendu, model trendu, výpočet reziduálních hodnot (zvlášť výpočet základní hodnoty z průběhu trendové funkce a interpolace odchylky) • Návrh vhodných parametrů pro vyhledávání s ohledem na existující síť (u některých metod nutné pro získání korektního výsledku, jinde pro zrychlení výpočtu): • volba sektorů • počet bodů ■ dosah hledání • tvar prohledávané oblasti (kruh, elipsa) DMT3 Volba metody a parametrů Prostorová interpolace dat • Ověření výsledku interpolace (verifikace) - hodnocení rozptylu, srovnání s jinými metodami, s původně naměřenými hodnotami: ručně provedená interpolace • bumerangová metoda (cross validation) itftiiMHBrararaim iHiaunt WMWm*] mapa chyb (střední chyba) změna průměrné hodnoty (rozptylu) velikost rozptylu u výsledku porovnání s původně naměřenými hodnotami • (požadavek exaktní interpolace) průběh pole DMT3 Prostorová interpolace dat Volba metody a parametrů Na čem tedy závisí výsledek ? • cíl interpolace (exaktní interpolace - aproximace - průběh pole) • počet a rozmístění známých měření • statistických charakteristikách zkoumaného souboru měření • vazbě principu zvolené metody a vývoje sledovaného pole Nelze spoléhat na 1 metodu í K DMT3 Prostorová interpolace dat DMT3 Kriging Unearvariogram model Minimum curvature Tension 0.0 Lv ABOS Sharp interpolation Normal interpolation Prostorová interpolace dat b-Efhhf Í1ÍA**T V ru 4tu.pt trnmlan (I.H t-------n-------n-------n- í--------Tí---------Tí--------T*--------O--------ÍT--------Tí--------7^~ HKTE Cn«ml lTit*rjM lat Infi) í----------TB----------T5----------O----------O----------37E----------ĎTB- AHS , 2008. Dressler, M. SurGe, gridding and mapping software. Manuál k programu [online]. Internet, < http://mujweb.cz/www/surge/ >, 2008. Golden software. Domovské stránky organizace [online]. Program Surfer, verze 8. Internet, < http://www.goldensoftware.com/products/surfer/surfer.shtml >, 2008. Homola, V. Sylaby geostatistiky a geoinformatiky [online]. Internet, < http://homel.vsb.cz/~hom50/ >, 2004. Horák, J. Prostorová analýza dat [online]. Internet, < http://gis.vsb.cz/pad/ >, 2002. Horák, J., Staněk, F. Interpolace prostorových dat. Seminář při mezinárodním sympoziu GIS Ostrava 2004. VŠB TU Ostrava 28.1.2004. DMT3 Prostorová interpolace dat Děkuji za pozornost. klimanek@mendelu.cz http://mapserver.mendelu.cz