4.1.1 Popis řešení projektu InGeoCalcza rok 2008

ÚvodV tomto popisu je obsažen postup a dosažené výsledky za rok 2008. Projekt je rozdělen jakmetodicky, tak i časově na tři hlavní fáze: teoretickou, implementační a ověřovací. Rok 2008 bylposledním rokem teoretické fáze. Řešily se v něm teoreticky nejnáročnější úkoly, které přineslydosud nejvýznamnější přidané hodnoty ke studovaným tématům. Současně byla s předstihempřípravována implementace navržených algoritmů, zejména u bayesovské klasifikace a analýzyposunů a deformací. Dále popíšeme podrobněji řešení jednotlivých aktivit vprůběhu roku 2008.

Bayesovská klasifikace digitálních obrazů

Řešení bayesovské klasifikace navržené v rámci aktivity A07-02 bylo dále rozvíjeno v souladus plánem. Byly zobecněny některé zjednodušující předpoklady stanovené ve zprávě o řešeníaktivity A07-02, zejména ty, které se týkají prvního kroku bayesovské klasifikace, tj. modelovánínepřesné hranice oblastí pomocí fuzzy množin.Popis celkového řešení modifikované bayesovské klasifikace obsahuje zpráva předložená za

rok 2007. Obecný matematický popis řešení byl také publikován v [4]. Proto se v této zprávěomezíme jen na první krok bayesovské klasifikace, kterým je modelování nepřesné hranice oblastípomocí fuzzy množin.Současně s teoretickým výzkumem modelování nepřesné hranice oblastí pomocí fuzzy mno-

žin byly podniknuty přípravné práce pro implementaci bayesovské klasifikace jako webové služy.Při tom byly vyvinuty základní operace pro interaktivní vkládání dat uživatelem a pro grafickouprezentaci výsledků bayesovské klasifikace. O této části se zmíníme v odstavci (1).Předpoklady o rozdělení pravděpodobnosti vrcholů polygonu určujícího tvar klasifikované

oblasti byly při aktivitě A07-02 v loňském roce stanoveny takto:

1. Rozdělení pravděpodobnosti vrcholů polygonů jsou isotropní neboli kruhová, tzn. majístejný rozptyl ve všech směrech.

2. Hustotou pravděpodobnosti vrcholů je polynomiální funkce odpovídající Pearsonovu roz-dělení typu II.

3. Topologie rozložení oblastí v obraze (tj. sousednost jednotlivých oblastí nezávisle na jejichtvaru) je známa předem, např. na základě předběžné, přibližné klasifikace některou zestandardních metod.

V letošním roce byly první dva předpoklady modifikovány tak, aby lépe vyhovovaly praktic-kým požadavkům a přitom byla zvýšena efektivita výpočtu parametrů tohoto rozdělení prav-děpodobnosti. Proto byla znovu přezkoumána možnost použití normálního rozdělení pravděpo-dobnosti.První předpoklad byl zobecněn tak, že již není nutné uvažovat pouze isotropní rozdělení

pravděpodobnosti vrcholů polygonů. Přesnost polohy jednotlivých vrcholů může tudíž být dánachybovými elipsami. Druhý předpoklad byl nahrazen obecnějším předpokladem o rozdělenípravděpodobnosti vrcholů polygonu. Nyní jím může být vícerozměrné normální rozdělení nenutně nezávislých veličin, tj. s plnou (nediagonální) kovarianční maticí. Před použitím tohotorozdělení pravděpodobnosti je však nutné provést transformaci k hlavním osám, aby sdruženáhustota pravděpodobnosti měla diagonální kovarianční matici.Původně zvolené Pearsonovo rozdělení typu II mělo sice výhodu v tom, že jeho hustota

je ve tvaru polynomu a dá se tudíž analyticky integrovat i při jeho vícenásobném použití,avšak výsledný výpočet je přesto velmi náročný. Příčinou této výpočetní náročnosti je postupnéhromadění dílčích podpůrných dat o mnoharozměrné množině, přes kterou se integruje. Tatomnoharozměrná množina je totiž při integraci rozčleněna na velké množství podoblastí, jejichžpříspěvky k celkovému integrálu jsou všechny stejně významné bez ohledu na jejich velikost.To znemožňuje zvyšovat efektivitu výpočtu při zvětšování počtu vrcholů polygonu. Tato potížs Pearsonovým rozdělením nebyla zpočátku patrná, neboť byly zkoumány jen obecné vlastnostinavrhované metody a ověřovány na jednoduchých polygonech.Zmíněná skrytá výpočetní potíž Pearsonova rozdělení se u normálního rozdělení tolik ne-

projevuje, neboť je definováno na nekonečném oboru a není tudíž nutno starat se o omezenýsuport, který v případě použití Pearsonova rozdělení způsobuje ono rozčlenění integrační ob-lasti na velké množství podoblastí. Místo toho však nastupuje problém s analytickou integracísdružené hustoty normálního rozdělení.

Formulace problému modelování nepřesné hranice oblasti

Nad rastrovým digitálním obrazem je dána oblast ve tvaru uzavřeného polygonu, jehož vrcholymají dvojrozměrná normální rozdělení pravděpodobnosti. Je třeba sestrojit fuzzy množinu,jejíž funkce příslušnosti udává pro každý pixel pravděpodobnost, že tento pixel náleží do danépolygonální oblasti.

Řešení problému modelování nepřesné hranice oblasti

Nejprve byl odvozen obecný vztah pro výpočet funkce příslušnosti fuzzy množiny. Tento vztahmá tvar zlomku, v jehož čitatel i jmenovateli je mnohonásobný určitý integrál sdružené hustotynormálního rozdělení. O takovýchto integrálech je známo, že nejsou řešitelné pomocí obvyklýchmatematických funkcí, a proto se řeší zásadně numericky. V našem případě však numerickéřešení není přípustné, neboť v dalším kroku bude třeba řešit inverzní úlohu, jak k danýmpravděpodobnostem nalézt rozdělení pravděpodobnosti vrcholů polygonální oblasti. Bylo proto

nutné vyjádřit oba integrály analyticky pomocí speciálních funkcí a případný aproximativnívýpočet provést až nakonec.

1 Přípravna implementace bayesovské klasifikace

Základní přípravné implementační práce zahrnují volbu trénovacích množin a nejjednoduššíverzi bayesovské klasifikace (viz. např. [1]). Na příkladě družicového snímku (viz. obr. 1), naněmž jsou dobře patrny dvě třídy: les a moře, bylo vybráno několik trénovacích množin protyto dvě třídy a provedena klasifikace. Na obrázku 2 jsou plně barevnými obdélníky vyznačenyzvolené trénovacích množiny.

Obrázek 1: Vstupní snímek pro bayesovskou klasifikaci

Obrázek 2: Výsledek bayesovské klasifikace

Závěr k bayesovské klasifikaci

Byly zavedeny dvě nové speciální funkce, které spolu s exponenciální funkcí a distribuční funkcíjednorozměrného normálního rozdělení umožňují vyjádřit vícenásobný integrál sdružené hus-toty normálního rozdělení analyticky. Tyto dvě pomocné speciální funkce byly definovány jakonekonečné polynomické řady ve dvou proměnných. Podmínky konvergence těchto nekonečnýchřad nebyly zatím teoreticky studovány.

Teoretické řešení analýzy posunů a deformací

K analýze posunů a deformací se přistupuje z hlediska metody kolokace. Při aplikaci metodykolokace se vychází z posunů v poli deformací, které byly určeny dvojím zaměřením souřadnicdaných identických bodů. Z nich se pak predikují posuny i v bodech, kde nebylo měřeno, tj.mimo identické body.Při tomto přístupu má klíčový význam stanovení statistické závislosti posunů v různých

místech deformačního pole. Tuto závislost zachycuje kovarianční matice, jejíž prvky je třebaodhadnout na základě zadaných posunů a dodatečných předpokladů. Základním předpoklademje pokles statistické závislosti dvou posunů se zvětšováním vzdálenosti příslušných identickýchbodů. Rychlost tohoto poklesu se při tradiční aplikaci metody kolokace udává pomocí tzv.kovarianční funkce. Tento tradiční přístup vychází z teorie gaussovských procesů a je v součas-nosti velmi podrobně rozpracován (viz. např. [2]). Přesto má použití kovarianční funkce některáomezení. Jedním z těchto omezení je nutnost předpokládat stejnou variabilitu všech identickýchbodů. Dalším nedostatkem kovarianční funkce jsou umělé předpoklady o závislosti vzdálenýchbodů skryté ve způsobu výpočtu nediagonálních prvků kovarianční matice i v samotné volbětvaru kovarianční funkce.Hlavním úkolem této aktivity je proto navržení alternativního postupu stanovení kovarianční

matice, který nahradí tradiční metodiku založenou na kovarianční funkci.

Formulace problému stanovení kovarianční matice

Jsou dány posuny v poli deformací, které byly určeny dvojím zaměřením souřadnic danýchidentických bodů. Je třeba odhadnout statistickou závislost zadaných posunů za předpokladu,že tato závislost klesá s rostoucí vzdáleností příslušných identických bodů. Míra poklesu tétozávislosti přitom musí být stanovena tak, aby se kovarianční matice daných i predikovanýchposunů přirozeně, plynule měnila v poli deformací.Požadavek přirozenosti a plynulosti změny kovarianční matice nelze předem jednoznačně

formulovat. Přesná matematická formulace tohoto požadavku proto vykrystalizuje až v průběhuřešení problému.

Řešení problému stanovení kovarianční matice

Kovarianční matici identických bodů označíme Cn, kde n značí počet identických bodů. Přidá-ním jednoho predikovaného bodu o souřadnicích x vznikne nová kovarianční Matice Cn+1(x).Matice Cn je řádu 2n, matice Cn+1(x) řádu 2n + 2. Souřadnice identických bodů určené vedvou etapách označíme po řadě xi, x′

i, kde i ∈ {1, . . . , n}. Posun na i–tém identickém bodě pakbude x′

i − xi.Kovarianční matici Cn+1(x) přísluší 2n+ 2 rozměrný elipsoid, jehož poloosy mají velikosti

rovné odmocninám vlastních čísel matice Cn+1(x) a směry poloos odpovídají směrům vlastníchvektorů matice Cn+1(x). Platí tedy

Cn+1(x) = U(x)V(x)U(x)T , (1)

kdeV(x) . . . diagonální matice vlastních čísel

U(x) . . . matice vlastních vektorů (sloupcových)

Pokud se predikovaný bod x ztotožní s některým identickým bodem xi, tzn. pokud ∃ i ∈{1, . . . , n} ; x = xi, pak se dvě poloosy elipsoidu příslušného kovarianční matici Cn+1(x)anulují. Tím vznikne degenerovaný elipsoid, který má o dva rozměry méně než při x 6= xi. Pokudse predikovaný bod x pouze přiblíží k některému identickému bodu xi, pak se příslušný elipsoidpouze zploští v oněch dvou poloosách, ale nezdegeneruje. To plyne z požadavku přirozenéa plynulé změny kovarianční matice. Tento požadavek lze geometricky interpretovat tak, abyse příslušný elipsoid plynule zplošťoval, nafukoval a natáčel.Na základě studia několika obvykle používaných kovariančních funkcí bylo vypozorováno

zplošťování a natáčení elipsoidu při pohybu predikovaného bodu x mezi identickými body.Přitom bylo zjištěno několik nápadných odchylek od plynulosti změny rozměrů a orientaceelipsoidu v prostoru. Korekce těchto odchylek umožnila zavést dvě vektorové funkce u, v tak,aby platilo

u(q,g(x)) = h(U(x)) , (2)

v(q,g(x)) = diag(V(x)) ,

g(x) . . . vektor vzdáleností predikovaného bodu x od identických bodů,

tj. g(x) = [|x− x1|, |x− x2|, . . . |x− xn|] ,h . . . vektorová funkce přiřazující vlastním vektorům kovarianční matice Cn+1(x)

(tj. sloupcům matice U(x)) jejich úhly se souřadnicovými osami,q . . . n-tice neurčených parametrů.

Současně musí platit, že elipsoid definovaný poloosami o velikostechV(xi) ve směrechU(xi)prostřednictvím funkcí u, v a rovností (2) je pro ∀ i ∈ {1, . . . , n} stejný a regulární (nedegene-rovaný). Matice jemu přidružené kvadratické formy je pak kovarianční maticí 2n-rozměrnéhonormálního rozdělení pravděpodobnosti. Tato kovarianční matice je tedy jednoznačně určenaaž na parametry q. Označíme ji proto C̃n(q).Zmíněné 2n-rozměrné normální rozdělení pravděpodobnosti dané kovarianční maticí C̃n(q)

je vlastně rozdělením pravděpodobnosti posunů x′i − xi. Toto rozdělení pravděpodobnosti má

nulové střední hodnoty všech jeho složek. Jeho hustota existuje, neboť kovarianční matice C̃n(q)má díky regularitě jejího elipsoid plnou hodnost. Tuto hustotu označíme fq. Parametry qmůžeme tedy bayesovsky odhadnout (viz např. [3]) pomocí zadaných posunů [x′

i − xi | i ∈{1, . . . , n} ], které mají sdružené normální rozdělení s hustotou fq.Protože parametry q byly bayesovsky odhadnuty, lze odvodit i jejich rozdělení pravděpo-

dobnosti. Jeho prostřednictvím pak s využitím rovností (2) a (1) odvodíme i rozdělení pravdě-podobnosti matice Cn+1(x) pro libovolný predikovaný bod x. Toto rozdělení pravděpodobnostije podobné Wishartovu rozdělení a je hledaným teoretickým výsledkem aktivity A08-02.

Odhad přesnosti tenzoru deformace simulací

Souběžně s teoretickými pracemi na stanovení kovarianční matice Cn+1(x) byl řešen alter-nativní způsob odhadu přesnosti tenzoru deformace pomocí počítačové simulace. Pro zadané

posuny byly pseudonáhodně generovány jejich chyby a registrován jejich vliv na výsledné tenzordeformace. Přitom bylo předpokládáno normální rozdělení pravděpodobnosti chyb posunů.Výsledky těchto simulací budou použity k ověření analytického způsobu odhadu přesnosti

tenzorů deformace. Tyto výsledky lze již nyní graficky prezentovat prostřednictvím webovéaplikace postupně vyvíjené v rámci tohoto projektu.Pro zadané posuny byly nejprve spočteny tenzory deformace pomocí bikubické interpolace

— viz. obr. 3.

Obrázek 3: Pole tenzoru deformace

Poté byly simulací napočítány střední chyby složek tenzoru deformace. Jejich grafické zná-zornění ve tvaru malých žlutých lichoběžníků ukazuje obrázek 4.

Obrázek 4: Přesnost tenzorů deformace

Jednorozměrná charakteristika tenzorového pole (tzv. čistý střih) je barevně znázorněna naobrázku 5.

Obrázek 5: Čistý střih tenzoru deformace

Závěr k analýze posunů a deformací

Byly navrženy dvě funkce u, v (viz (2) a následující postup), které řídí chování kovariančníhoelipsoidu a tím i jemu příslušné kovarianční matice Cn+1(x). Jedna z těchto řídicích funkcí(u) udává závislost rozměrů elipsoidu na poloze predikovaného bodu, druhá řídicí funkce (v)pak ovládá orientaci elipsoidu vůči souřadnicovým osám. Parametry obou dvou řídicích funkcíjsou odvozeny pouze od variability daných identických bodů, takže se do celkového postupu

nepromítají žádné umělé předpoklady, které by postrádaly statistickou interpretaci. K výslednékovarianční matici Cn+1(x) lze odvodit i její rozdělení pravděpodobnosti, což vede k většíobjektivitě při odhadu přesnosti metody kolokace.Byl připraven nástroj k ověření analytického způsobu odhadu přesnosti tenzorů deformace

založený na počítačové simulaci. Jeho grafická prezentace (viz. obrázky 3, 4, 5) bude využitai k prezentaci analytického odhadu přesnosti.

Literatura

[1] D. G. T. Denison, C. C. Holmes, B. K. Mallick, and A. F. M. Smith. Bayesian methods fornonlinear classification and regression. John Wiley Sons, 2002.

[2] T. Hida and M. Hitsuda. Gaussian processes. American Mathematical Society, 1993.

[3] Karl-Rudolf Koch. Bayesian Inference with Geodetic Applications, volume 31 of LectureNotes in Earth Sciences. Springer–Verlag, 1990.

[4] L. Soukup. Overall uncertainty of georeferencing and classification. In A. Stein, editor,International Symposium on Spatial Data Quality. ISSDQ ’07, Enschede, June 2007. ITCEnschede.

Tíhové redukce geodetických měření 

1 Vliv tíhového pole Země na měřené veličiny v geodézii Výsledky pozemní geodetická měření jsou ovlivněny proměnlivostí zemského tíhového pole. Této problematice byla v minulosti věnována značná pozornost, protože přesnost klasicky budovaných polohových základů silně závisela právě na opravách měřených veličin z vlivů tíhového pole Země. S nástupem technik kosmické geodézie do budování polohových základů již tíhové korekce nehrají v geodézii tak klíčovou roli, jako dříve. Stále je jejich použití nutné při některých speciálních pracích, například v důlním měřictví [1]. Přestože cílem projektu je výpočet pouze redukcí z tíhového pole Země, aplikace bude umožňovat úplnou transformaci uvedených měřených veličin na výpočetní plochu, kterou je elipsoid, aby nebylo nutné počítat další nezbytné korekce zvlášť v jiném programu. Aplikace tedy počítá všechny „matematické“ korekce k níže uvedeným veličinám. Program nepočítá žádné korekce označované jako „fyzikální“, jako je vliv teploty, refrakce apod. O tyto korekce by měla být měřené opravena před zavedením korekcí matematických.

1.1 Korekce směrníku Opravy směru měřeného z bodu Pi[φi, λi, hi] na bod Pj[φj, λj, hj] z vlivu tížnicových odchylek lze vyjádřit podle [3] jako

 kde ηi příčná složka tížnicové odchylky v bodu Pi, ξi meridiánová složka tížnicové odchylky v bodu Pi, αij azimut z bodu Pi na bod Pj, zij zenitová vzdálenost z bodu Pi na bod Pj.

stačí podle [4] určit s dostatečnou přesností jako

kde sij prostorová (šikmá) vzdálenost bodů Pi, Pj. Další zavedené korekce (které již nesouvisí přímo s tíhovým polem Země) jsou korekce z výšky cíle nad elipsoidem C3 a směrníku normálového řezu na směrník geodetiky C4, které se spočítají podle [3] jako

 kde Mi meridiánový poloměr křivosti v bodě Pi, Ni příčný poloměr křivosti v bodě Pi, e excentricita použitého elipsoidu φm střední šířka (= 1/2(φi + φj)) Nm střední příčný poloměr křivosti (= 1/2(Ni + Nj))  Směrník na elipsoidu σE se pak spočítá z měřeného směrníku σ‘ jako

σE = σ‘ + C1 + C2 + C3 + C4

1.2 Korekce horizontálních úhlů Redukci horizontálního úhlu snadno odvodíme z vzorce pro redukci směru. Od redukce směru na záměru vpřed odečteme redukci směru spočítanou na záměru vzad. Člen C1 se tím eliminuje a zůstanou rozdíly členů C2, C3 a C4. Úhel na elipsoidu ωE tedy spočítáme z měřeného úhlu ω‘ podle vzorce

ωE = ω‘ + FC2 + FC3 + FC4 - BC2 + BC3 + BC4

kde index F označuje korekce na směr záměry vpřed a index B korekce spočítané pro směr na záměru vzad.

1.3 Korekce vzdáleností Měřená vzdálenost mez dvěma body nijak nezávisí na tíhovém poli Země. Žádné redukce v souvislosti s tíhovým polem není třeba zadávat. Aplikace však přesto bude umožňovat redukci měřených šikmých i vodorovných délek na elipsoid, aby měl uživatel jistotu, že všechny nutné redukce všech měřených veličin budou provedeny touto aplikací a vyhnul se tak nejistotě, jaké korekce byly zavedeny a které musí případně dopočítávat v jiném programu. Je-li měřena prostorová (šikmá) vzdálenost sij, nejprve se převede na vodorovnou vzdálenost lij podle vztahu

Nísledně se tato vzdálenost převede na vzdálenost v nulové hladině [3]

kde

Ze vzdálenosti pak spočítáme délku příslušného normálového řezu na elipsoidu vzorcem

Rozdíl mezi délkou normálového řezu a délkou geodetiky lze rovněž vyjádřit (např. [3]), ale projeví se až u délek přesahujících 1000 km a na území ČR může být tato korekce spolehlivě zanedbána, a proto ji v této aplikaci neuvažujeme.

1.4 Vliv na zenitovou vzdálenost (úhel) Korekce měřené zenitové vzdálenosti zij‘ z bodu Pi na bod Pj se určí jako [3]

a zenitová vzdálenost odpovídající normále k elipsoidu je pak

2 Formát a vlastnosti vstupních dat Pro výpočet jsou nutné tři druhy vstupních dat: 1) Geodetické souřadnice bodů, z nichž bylo měřeno (a na něž bylo

měřeno). 2) Měřené veličiny, kterými mohou být:

a) směrníky, b) úhly, c) vodorovné či šikmé (prostorové) vzdálenosti, d) zenitové úhly.

3) Obě složky tížnicové odchylky (meridiánová a příčná). Pro vstup souřadnic bodů a měřených veličin bylo nutno zvolit vhodný formát. Formátů pro geodetická měření existuje celá řada. Pro účely projektu se nejvíce hodí otevřený formát, který bude možné snadno sestavit i bez použití specializovaného software a který bude snadno čitelný. Zároveň by měl formát umožňovat i zpracování i většího množství dat z výstupu nějakého běžně používaného geodetického programu. Všechny tyto podmínky splňuje XML formát GKF používaný v projektu GNU Gama [2]. Formát je jednoduchý, čitelný a existuje podpora pro transformaci z mnoha jiných formátů (např. prostřednictvím SW Kokeš).

2.1 Geodetické souřadnice bodů Souřadnice bodů budou zadávány v S-JTSK (závazný systém na území ČR) a v XML souboru jsou uloženy v tagu <point>, který má strukturu: <point id=STR, x=NUM, y=NUM, z=NUM, fix=STR, adj=STR>

(Podrobný popis viz [2].) Pro naší aplikaci není rozhodné, které body jsou považovány za pevné, opěrné či neznámé, proto atributy fix a adj nevyužijeme. Použití vzorců však předpokládá znalost geodetických souřadnic, souřadnice tedy budou nejprve transformovány. Znalost geodetických souřadnic však stačí pouze přibližná (tížnicové odchylky máme v rastru 30“x30“ a jejich přesnost a rozlišení je zde limitujícím faktorem), a k transformaci tak postačí použít vhodný globální transformační klíč platný pro území ČR. Důležité je, aby všechny body měly v tagu point uvedeny všechny souřadnice, tj. x, y i z, tedy včetně výšky, jinak nemůže výpočet proběhnout. I u výšky však platí, že pro výpočet korekcí stačí znát její hodnotu přibližně. V praxi bude tento postup nutno provádět iteračně: nejprve se vypočtou přibližné souřadnice bodů bez zavedení tíhových oprav (např. programem Gama), následně se měření opraví s využitím přibližných souřadnic pomocí naší aplikace a pak se zopakuje vyrovnání s opravenými měřenými hodnotami. Zadané normální nadmořské výšky budou transformovány na elipsoidální pomocí lokálního modelu kvazigeoidu. (Podle [4] by bylo možné dokonce rozdíl mezi nadmořskými a elipsoidálními výškami zanedbat).

2.2 Měřené geodetické veličiny XML formát GKF umožňuje zadat všechny měřené veličiny, pro něž jsme uvedli korekce z tíhového pole i další matematické korekce. Přehled XML tagů odpovídající jednotlivým veličinám je uveden v tab. 1. veličina XML tag důležité atributy tagu směr direction from, to, val vodorovný úhel angle from, bs, fs, val zenitová vzdálenost

z-angle from, to, val

vodorovná délka distance from, to, val šikmá délka s-distance from, to, val

tab. 1 - Přehled označení veličin v XML

U některých korekcí je nutno znát azimut ze stanoviště na měřený bod. V těchto případech se nebudou používat měřené hodnoty, ale azimut se spočítá ze zadaných souřadnic bodů. Azimut totiž nebývá většinou přímo měřenou veličinou, při určování azimutu z měření by bylo nutné měření určitým způsobem vyrovnávat, což není účelem této aplikace.

obr. 1 - Příčná složka tížnicové odchylky η

2.3 Složky tížnicové odchylky Obě složky tížnicové odchylky budou počítány z přednostně z lokálního, případně z globálního modelu (způsob výpočtu byl vyřešen v rámci předchozích aktivit) a nebude zde znovu podrobně zmiňován. Hodnoty složek tížnicových odchylek uložené v lokálním modelu jsou znázorněny na obr. 1 a obr. 2.

3 Formát výstupních dat Výstupem z programu bude: 1. XML soubor GKF, který bude shodný se vstupním souborem pouze s tím

rozdílem, že příslušné měřené veličiny budou opraveny o korekce. 2. Protokol o výpočtu korekcí. Protože použitý XML formát neobsahuje tagy

nebo atributy, do nichž by bylo možné uvést přímo hodnoty jednotlivých vypočtených korekcí a není žádoucí tento XML formát z důvodů kompatibility s jinými programy rozšiřovat, bude výstupem dále protokol, ve kterém budou uvedeny podrobnosti o vypočtených korekcích. Protokol bude obsahovat u každého měření vyčíslené jednotlivé korekce podle vztahů uvedených v úvodní kapitole této zprávy. Vhodný formát protokolu bude upřesněn až při implementaci.

obr. 2 - Meridiánová složka tížnicové odchylky ξ

4 Výpočetní postup Výpočet bude probíhat podle níže uvedeného schématu. Jednotlivé kroky podle popisu a výpočetních vztahů uvedených v tomto textu.

5 Literatura [1] MUČKOVÁ, J. Vliv redukcí z tížnicových odchylek na přesnost základní orientační úsečky při připojovacím a usměrňovacím měření gyroteodolitem v dole. Sborník vědeckých prací Vysoké školy báňské – Technické univerzity Ostrava, řada hornicko-geologická, Volume L (2004), No.2, p. 33-38, ISSN 0474-8476. [2] ČEPEK, Aleš. GNU Gama 1.9.06a [online]. [2009] [cit. 2009-02-05]. Dostupný z WWW: <http://www.gnu.org/software/gama/manual/gama.pdf>. [3] VANÍČEK P., KRAKIWSKY E. Geodesy: the concepts. Second edition. Elsevier Publishers, Amsterdam. 1986. [4] CIMBÁLNÍK, M., MERVART, L. Vyšší geodézie 1. Geometrická. Vydavatelství ČVUT, Praha. 2002.

Harmonická syntéza do velmi vysokého řádu a stupně 

1 Úvod V dubnu 2009 byl vydán nový model globálního tíhového pole Země EGM08 s koeficienty do stupně a řádu 2190, resp. 2160, jehož přesnost na území ČR je díky použití pozemních tíhových dat s rozlišením 5’x5’ velmi vysoká (viz aktivita „Srovnání kvazigeoidu EGM08 a lokálního modelu“). Z tohoto modelu je možné generovat celou řadu parametrů tíhového pole Země, jak to bylo plánováno v návrhu projektu, a to s poměrně vysokou přesností. Model je volně k dispozici, je uložen v textovém souboru, který obsahuje koeficienty rozvoje Cnm, Snm včetně jejich středních chyb. Z modelu lze spočítat celou řadu parametrů tíhového pole Země, přesný výpočet však předpokládá hlubší znalosti fyzikální geodézie (je nutno transformovat souřadnice do různých systémů, zohledňovat různé parametry GM, a modelu EGM08 a použitého referenčního elipsoidu, vybrat mezi mean-tide nebo zero-tide systémem atd.). Navíc se jedná o výpočetně náročnou úlohu, a to jak z hlediska přesnosti (při generování Legendreových funkcí hrozí podtečení či přetečení reálných čísel ve formátu double), tak z hlediska rychlosti a paměťové náročnosti. Jenom velikost souboru s koeficienty přesahuje velikost 200 MB.

V současnosti navíc autorům projektu není znám spolehlivý a uživatelsky jednoduchý program, který by umožnil výpočet parametrů z tohoto modelu. Program Synth [1] volně poskytovaný spolu s daty je k dispozici pouze ve formě zdrojového kódu ve Fortranu 77. Zadávání dat do programu se provádí modifikací zdrojového kódu a překladem programu. Program má navíc omezený počet bodů, které je možné najednou zpracovat, a nezapisuje pod OS Windows spolehlivě výsledky do výstupních souborů. Program Gravsoft [2] funguje spolehlivě, ale je velmi pomalý a má velmi omezené možnosti nastavení.

Popisovaná aplikace umožní provést výpočty požadovaných parametrů tíhového pole Země i těm uživatelům, kteří problematice generování parametrů z modelu podrobně nerozumí, či těm, kteří potřebují data pouze menšího rozsahu. Protože aplikace poběží na straně serveru, nebude umožňovat výpočty ve vysokém rozlišení pro velká území, která jsou časově velice náročná a neúměrně by přetěžovala server. Přesné limity budou stanoveny ve fázi testování.

2 Volba algoritmu Rešerší dostupné literatury zabývající se generováním přidružených Legendreových funkcí byl zvolen jako nejvhodnější algoritmus „Forward column“ (FC) algoritmus [3] s numerickým řešením problému podtečení rozsahu pomocí škálování. Pro stupeň a řád použitý v současném modelu EGM08 stačí konstantní škálování (přibližně do stupně a řádu 2800). Zároveň byla testována metoda dynamického škálování podle [4], kterou je možné užít i pro případné budoucí modely tíhového pole Země do vyšších řádů a stupňů (funguje spolehlivě nejméně do stupně a řádu 20000). Popis obou algoritmů je uveden v citovaných publikacích a nebude zde podrobně uváděn.

Pro účely aplikace byl však algoritmus mírně modifikován následujícím způsobem: Podle [4] mají normované Legendreovy funkce tvar (pro bod o sférických souřadnicích P[r, ])

Tento tvar je výhodný při výpočtech vysoko nad zemským povrchem, kdy může hodnota

dosáhnout velmi vysokých hodnot a hrozí přetečení. Naše aplikace bude sloužit pro výpočty na povrchu Země či v jejím blízkém okolí, kde lze ukázat, že přetečení z tohoto důvodu nehrozí. Proto postačí tvar normovaných ALF

který má tu výhodu, že normované ALF jsou závislé pouze na polární vzdálenosti a nikoliv na velikosti průvodiče r, což výrazně redukuje počet nutných výpočtů normovaných ALF při výpočtech v rastru na povrchu Země, kdy se pro danou zeměpisnou šířku mění r ale zůstává konstantní.

3 Implementace Protože pro uvedený algoritmus neexistuje vhodná knihovna, bylo přistoupeno k jeho implementaci. Algoritmus byl implementován nejprve v jazyce Matlab, ve kterém byla ověřena jeho spolehlivost a přesnost. Metoda dynamického škálování byla rovněž implementována. Protože testy potvrdily vysokou přesnost a efektivnost metody konstantního škálování pro model do stupně a řádu 2190, bylo přistoupeno k její finální implementaci v jazyce C++.

Po úspěšném testování přesnosti implementace v C++ byla věnována značná pozornost časové optimalizaci programu. V neoptimalizovaném kódu trvá výpočet ALF pro jeden bod přibližně 2 s, po optimalizaci se podařilo dosáhnout času přibližně 0,07 s / bod na běžném PC. Na výkonném serveru, kde aplikace poběží, lze předpokládat ještě výrazně vyšší rychlost.

Hlavními metodami optimalizace bylo:

1. Důsledné využívání všech jednou vypočítaných hodnot i ve všech následujících výpočtech (např. normovací koeficienty).

2. Sekvenční přístup do polí v paměti pro lepší využití vyrovnávací paměti procesoru.

3. Úprava některých vzorců do rekurzivní podoby (např. časově náročné umocňování je nahrazeno postupným násobením).

Výsledkem optimalizace je velice rychlý program, což ukázaly i testy (viz dále).

Detaily implementace jsou patrné z přiložených zdrojových kódů a podrobné dokumentace programu.

4 Testy Protože není možné testovat přímo hodnoty ALF (chybí spolehlivý srovnávací SW), byl program rozšířen a nyní již umožňuje nejen spočítat hodnoty ALF, ale i provést syntézu a spočítat některé parametry tíhového pole Země (zatím však pouze ve speciálních případech pro účely testování s nutností modifikace zdrojového kódu, obecné řešení bude řešeno v dalších aktivitách). Níže uvedené testy tak testují nejen generování samotných ALF, ale i celou řadu výpočetních postupů použitých v aplikaci (transformace souřadnic, korekce C20 členu, výpočty parametrů normálního pole GRS80 elipsoidu…).

4.1 Testy přesnosti výpočtu K porovnání byly použity programy Synth [1] a Gravsoft [2]. Jak již bylo zmíněno, ani jeden z těchto programů neumožňuje přímý výstup samotných ALF, proto byly porovnány výsledky vhodné z nich z nich odvozené veličiny. Takovou veličinou je výšková anomálie ζ, protože její výstup umožňují oba programy a zároveň je výpočet z ALF a koeficientů modelu Cnm, Snm poměrně jednoduchý. Navíc má tato veličina rozměr [m], což usnadní kvantifikaci případných rozdílů ve výsledcích.

Oba programy dávají různé výsledky, lišící se přibližně o 1 m. Je to způsobeno použitím jiného systému pro Stokesův koeficient modelu C20 (zero-tide, resp. mean-tide) a také faktem, že Synth při výpočtu zanedbává vliv členu nultého řádu v rozvoji (kvůli neznalosti přesného parametru GM jej nelze spočítat přesně), kdežto Gravsoft tento člen počítá (i když pouze přibližně). Oba výsledky lze tedy považovat za správné. Naše aplikace umožňuje oba režimy výpočtu a je tak možné srovnání s oběma programy. Testování bylo provedeno na 10 bodech o různé zeměpisné šířce, délce a výšce rovnoměrně rozložených na území ČR a v blízkém okolí. Výsledky shrnuje tab. 1. Protože výsledky vychází absolutně shodně, testování na dalších bodech není považováno za nutné. Některé chyby způsobené podtečením a přetečením se projevují pouze v některých zeměpisných šířkách. Proto byly udělány další testy, které ukazují opět na absolutně shodné výsledky i pro tyto

problematické šířky, pro účely naší aplikace však nemají význam a nejsou zde uvedeny.

souřadnice bodu (GRS80) ζ ζ' φ λ h Gravsoft GeoCalc rozdíl Synth Geocalc rozdíl 49 13 300 46.279 46.279 0.000 47.195 47.195 0.000 50 13 500 46.378 46.378 0.000 47.292 47.292 0.000 51 13 800 43.985 43.985 0.000 44.898 44.898 0.000 49 15 300 45.601 45.601 0.000 46.517 46.517 0.000 50 15 500 43.932 43.932 0.000 44.846 44.846 0.000 51 15 800 42.023 42.023 0.000 42.936 42.936 0.000 49 17 300 42.807 42.807 0.000 43.722 43.722 0.000 50 17 500 43.471 43.471 0.000 44.386 44.386 0.000 51 17 800 40.261 40.261 0.000 41.173 41.173 0.000 50 18 500 42.026 42.026 0.000 42.941 42.941 0.000

tab. 1 - Porovnání přesnosti výpočtu

4.2  Testy rychlosti výpočtu Neoptimalizovaný výpočet parametrů tíhového pole Země trvá pro jeden bod přibližně 2 s (viz výše). Rychlost výpočtu je proto velmi důležitá. Provedli jsme proto srovnání rychlosti všech tří použitých programů. Výsledky shrnuje tab. 2 a graf 1.

program / počet bodů 1 10 20 50 200 500 GeoCalc 0.8 1.4 2.0 4.1 14.9 35.5 Gravsoft 7 32 55 115 Synth 35 36 37 42 64 109

Tab. 2 - Porovnání rychlosti výpočtu v [s]

Graf 1 – Závislost výpočetního času jednotlivých programů na počtu bodů v [s]

Výpočet byl prováděn pro náhodně rozmístěné body, kdy jednotlivé programy nemohly využít předchozích výsledků pro další body. Pokud se provádí výpočet v rastru geodetických souřadnic, stačí počítat časově nejnáročnější ALF pouze pro každou zeměpisnou šířku zvlášť. Neplatí to však již pro rastr v souřadnicích S-JTSK a program Synth této vlastnosti navíc dokáže využít pouze v případě, že jsou všechny body stejně vysoko nad referenční koulí, což výrazně limituje možnost využití tohoto způsobu zrychlení výpočtu. Srovnání výpočtu pro jednotlivé body lze proto považovat za nejobjektivnější porovnání rychlosti výpočtu.

Výsledky jednoznačně prokazují efektivitu implementovaného algoritmu, který je více než 2x rychlejší než druhý nejrychlejší program Synth. Synth by navíc bylo v aplikaci obtížné použít kvůli velmi dlouhému spouštění programu, který trvá přes 30 s i při výpočtu jediného bodu. To je způsobeno patrně pomalým načítáním souboru Stokesových koeficientů, přestože byl celý tento soubor při všech testech uložen v paměťové cache a nebylo jej nutno načítat z pevného disku.

4.3 Shrnutí testů Provedené testy ukazují na vynikající přesnost i rychlost výpočtu poruchového potenciálu a výškové anomálie. Výsledky přesně odpovídají srovnatelným programům (Gravsoft a Synth) a to při několikanásobně (až řádově) kratším výpočetním čase.

5 Závěr pro výpočet normovaných ALF byl zvolen a implementován algoritmus, jehož vysoká kvalita a efektivita byla potvrzena srovnávacími testy s obdobnými programy. Program je předán ve formě zdrojových kódu včetně obsáhlé dokumentace i ve formě spustitelného souboru. Testy jsou předány ve formě textových dokumentů včetně tabulek a grafů.

6 Literatura [1] Holmes S. A., Pavlis N. K. A Fortran program for very-high-degree harmonic synthesis. Version 05/01/2006. [online]. [2007] [cit. 2009-02-05]. Dostupný z WWW <http://earth.esa.int/gut/prototype/Harmonic_Synth/>

[2] Forsberg R., Tscherning, C. C. An overview manual for the GRAVSOFT. [online]. [2003] [cit. 2009-02-05]. Dostupný z WWW <http://cct.gfy.ku.dk/publ_cct/cct1792.pdf>

[3] Bethencourt A. , Wang J., Rizos C., Kearsley A. H. W. 2005. Using personal computers in spherical harmonic synthesis of high degree earth geopotential models. Dynamic Planet 2005, Cairns, Australia.

[4] Nesvadba O. Towards the numerical evaluation of high degree and order associated Legendre functions as in EGM08. 2008. IAG International Symposium "Gravity, Geoid and Earth Observation", Chania, červen 2008.

[5] PAVLIS N. K., HOLMES S. A., KENYON S. C., FACTOR J. K. An Earth Gravitational model to degree 2160: EGM2008. Preseneted at 2008 General Assembly of European Geosciences Union, Vienna, Austria, April 13-18, 2008.

Srovnání kvazigeoidu EGM08 a lokálního modelu 

1 Lokální a globální modely kvazigeoidu V roce 2008 byl představen nový globální model tíhového pole Země EGM08 [1] ve formě Stokesových koeficientů do stupně a řádu 2160. Nahradil předchozí oficiální model EGM96, který obsahoval koeficienty pouze do stupně a řádu 360. Prostorové rozlišení modelu EGM96 je přibližně 30’ × 30’, prostorové rozlišení modelu EGM08 je již 5’ × 5’. Pro sestavení takto podrobného modelu byla použita mimo jiné i pozemní tíhová data v odpovídajícím rastru (tj. 5’ × 5’) [1].

Rozlišení a přesnost modelu EGM96 byly pro přesný výpočet kvazigeoidu nedostatečné, proto byly počítány přesnější lokální modely kvazigeoidu z místně dostupných tíhových dat s vyšším rozlišením (pro území ČR např. [2], [3]). Z testování předběžné verze nového modelu EGM08 porovnáním rozdílů výškových anomálií na přibližně 1000 bodech výběrové údržby na území ČR Δζ = ζ – ζ’ (1) kde ζ jsou výškové anomálie vypočítané z modelu EGM08 a ζ’ výškové anomálie ζ’ = h – H (2) kde h je geodetická výška (změřená GPS) a H je normální výška (změřená nivelací) vyšla podle [4] střední kvadratická chyba RMSE Δζ = 3,6 cm (3) což ukazuje poměrně vysokou shodu mezi hodnotami ζ, ζ’. Naskýtá se otázka, zda lze takto dobrou shodu v ČR vůbec ještě vylepšit pomocí lokálních modelů kvazigeoidu, protože testy lokálních modelů kvazigeoidu pomocí obdobné metodiky (kde ovšem hodnota ζ není počítána z globálního modelu, ale z přesného lokálního modelu kvazigeoidu) vykazují střední kvadratickou chybu na území ČR (např. [3], [6]) rovněž kolem 3 cm.

Tento článek ukazuje odhad velikosti vlivu reziduálních tíhových anomálií na kvazigeoid, kde reziduální tíhové anomálie Δgres počítáme jako Δgres = Δgloc - ΔgEGM (4) kde Δgloc je měřená tíhová anomálie interpolovaná do lokálního modelu (ve formě rastru) a ΔgEGM je tíhová anomálie spočítaná z globálního modelu (ve formě Stokesových koeficientů).

2 Použitá data Na obrázku 1 je model tíhových anomálií Δgloc interpolovaný z pozemních tíhových měření v rastru s rozlišením 30˝ × 30˝ (lokální model VÚGTK [5]). Na dalších obrázcích jsou pro srovnání vykresleny tíhové anomálie jak pro nejnovější model

EGM08, tak pro starší EGM96. Rozdíl v podrobnosti původního modelu EGM96 a nového EGM08 na území České republiky je patrný z map tíhových anomálií na obrázcích 2 a 3, na kterých jsou znázorněny jak celé hodnoty tíhových anomálií spočítaných z příslušných modelů (ΔgEGM96 a ΔgEGM08), tak reziduální tíhové anomálie spočítané pro oba modely (Δgres96 a Δgres08). Všechny tři modely tíhových anomálií (Δgloc, ΔgEGM96 i ΔgEGM08 ) jsou spočítané na povrchu země v rastru s rozlišením 30˝ × 30˝. Všechny modely jsou tedy spočítané ve stejných bodech, a přestože globální modely nemají takto vysoké rozlišení, je tento postup vhodný kvůli závislosti tíhových anomálií na nadmořské výšce. Přehled hodnot je uveden v tabulce 1.

Obr. 1 Tíhové anomálie z lokálního modelu (Δgloc) [mGal]

Obr. 2 Tíhové anomálie z EGM08 ΔgEGM08 (vlevo) a reziduální tíhové anomálie Δgres08 (vpravo),

oba [mGal]

Obr. 3 Tíhové anomálie z EGM96 ΔgEGM96 (vlevo) a reziduální tíhové anomálie Δgres96 (vpravo),

oba [mGal]

Δgloc ΔgEGM08 ΔgEGM96 Δgres08 Δgres96

min -53,5 -50,0 -45,1 -88,7 -59,9

max 209,9 159,4 69,8 83,5 159,6

mean 13,1 15,1 16,3 -1,9 -3,2

std 24,6 23,8 20,5 6,2 13,4

Tab. 1 Statistika tíhových anomálií z různých modelů a jejich rozdílů na území zobrazeném na obr. 1 – 3 (celkem 720000 bodů), hodnoty v [mGal].

3 Vliv reziduálních tíhových anomálií na kvazigeoid K odhadu vlivu reziduálních pozemních tíhových dat nového modelu Δgres08 využijeme tzv. „remove-restore“ princip: výškovou anomálii ζ rozdělíme na dvě části ζ = ζEGM + ζres (5) kde ζEGM je část generovaná tíhovými anomáliemi ΔgEGM08 a ζres je část generovaná reziduálními tíhovými anomáliemi Δgres08 a měla by představovat vysokofrekvenční část řešení, která není obsažena v EGM. Výškové anomálie počítám pouze pro globální model EGM08, a proto je budu označovat pouze ζEGM a ζres bez uvedení přípony „08“, kterou jsem používal u tíhových anomálií k odlišení anomálií z modelu EGM08 od anomálií z modelu EGM96.

Pro testování ovšem nebylo z výpočetních důvodů použito území celé ČR, ale pouze území omezené území 14˚ – 17˚ východní délky a 49,5˚ - 51˚ severní šířky. Jedná se přibližně o území středních, severních a východních Čech a části severní Moravy.

3.1 Výpočet kvazigeoidu z globálního modelu Výškovou anomálii ζEGM z globálního modelu určíme podle vztahu

),,(),,(),,(

hhThEGM

λϕγλϕλϕζ = (6)

kde ϕ, λ, h jsou geodetické souřadnice, ),,( hλϕγ je velikost normálního zrychlení, T(ϕ, λ, h) je poruchový potenciál spočítaný řadou (např. [7])

∑=

⎟⎠⎞

⎜⎝⎛=

max

2),(),,(

N

nn

n

Tra

rGMrT λθλθ , (7)

kde λθ ,,r jsou sférické souřadnice (odpovídající bodu o geodetických souřadnicích ϕ, λ, h), a je velikost hlavní poloosy použitého elipsoidu a GM je geocentrická gravitační konstanta, Nmax je v našem případě 2160, dále

( )

( ) nPmSmCPJCT

nPmSmCTn

mnmnmnmnnnn

n

mnmnmnmn

sudé pro)(sinsincos)(sin)(),(

liché pro)(sinsincos),(

10

0

∑

∑

=

=

⋅++−=

⋅+=

θλλθλθ

θλλλθ (8)

nmC , nmS jsou (plně normalizované) Stokesovy koeficienty modelu EGM08, nmP jsou plně normované přidružené Legendrerovy funkce, nP plně normalizované Legendrerovy polynomy a

nJ2jsou plně normované Stokesovy koeficienty

normálního pole (např. [7]).

Výšková anomálie ζEGM je znázorněna na obr. 4.

Obr. 4 Výšková anomálie z EGM08 ζEGM [m]

Z globálního modelu musíme kromě výškové anomálie spočítat i tíhovou anomálií ΔgEGM08, kterou budeme potřebovat pro výpočet reziduální tíhové anomálie Δgres08. Tíhovou anomálii spočítáme řadou

),()1(),,(

22 λθλθ n

nNmax

nT

ran

rGMrg ∑

=⎟⎠⎞

⎜⎝⎛−=Δ

(9)

3.2 Výpočet kvazigeoidu z lokálních dat K dispozici máme hodnoty tíže g na povrchu země popsané v [5]. Z těchto hodnot nejprve spočítáme tíhové anomálie na volném vzduchu Δgloc podle vztahu (10) ),(),,(),,( HHgHg loc ϕγλϕλϕ −=Δ

kde ϕ, λ jsou geodetické souřadnice a H je normální výška. Tyto tíhové anomálie jsou redukovány o tíhové anomálie z globálního modelu ΔgEGM08 na reziduální tíhové anomálie Δgres. Z hodnot tíhových anomálií Δgres a normálních výšek H lze spočítat odpovídající výškovou anomálii ζres na povrchu Země různými postupy (např. [7]).

Výšková anomálie ζres je znázorněna na obr. 5. Obrázek výsledného modelu ζ = ζEGM + ζres není samostatně uveden, protože je velmi podobný obr. 4. Na obr. 7 je ovšem model ζ zobrazen, a to spolu s body výběrové údržby na nichž byl testován.

Obr. 5 Reziduální výšková anomálie (vzhledem k modelu EGM08) ζres [m]

Z obrázku a ze statistických parametrů souboru ζres (střední hodnota 0,000 m; směrodatná odchylka 0,028 m, minimum -0,101 m a maximum 0,136 m) je zřejmé, že tyto hodnoty ovlivňují průběh kvazigeoidu v řádu centimetrů, přičemž extrémní hodnoty na sledovaném území dosahují až velikosti 14 cm. Pokud požadujeme model kvazigeoidu s přesností v řádu centimetrů, příspěvek ζres nelze zanedbat.

3.3 Porovnání kvazigeoidů V předchozím odstavci jsme ukázali, že rozdíl mezi kvazigeoidem spočteným pouze z globálního modelu a kvazigeoidem doplněným o složku ζres dosahuje hodnot až 14 cm. V kapitole 1 však citacemi z různých zdrojů ukazuji, že při testování kvazigeoidu ζEGM spočítaného pouze z EGM08 na bodech výběrové údržby se dosáhlo podobné shody, jako při testování lokálních modelů vypočtených z lokálních dat (a tedy zahrnujících i složku ζres).

Provedl jsem proto rovněž testování jak tíhových anomálií ζEGM , tak tíhových anomálií ζ = ζEGM + ζres na bodech výběrové údržby. Do testovacího území padlo celkem 359 bodů výběrové údržby. Na těchto bodech lze spočítat výškovou anomálii ζ’ jako rozdíl geodetické a normální výšky. Hodnoty ζ na konkrétním bodě lze určit interpolací z rastru hodnot ζ a hodnoty ζEGM výpočtem řady z koeficientů modelu EGM08. Ta získáme pro každý bod výběrové údržby tři hodnoty výškové anomálie ζ, ζEGM a ζ’. Lze očekávat, že soubor výškových anomálií ζ’ bude systematicky posunut o jistou hodnotu vůči souborům ζEGM a ζ

z důvodu rozdílnosti nulové hladinové plochy W0 modelu EGM08 a počátku výškového systému použitých normálních výšek.

Ukázalo se, že na několika bodech výběrové údržby hodnoty ζ’ hrubě nesouhlasí s hodnotami ζEGM a ζ. Tyto body ze souboru odstraněny (odstraněno bylo celkem 6 bodů) Tabulka 2 ukazuje jak statistiku rozdílů mezi hodnotami ζ’ - ζEGM a mezi hodnotami ζ’ – ζ a to jak pro soubor všech bodů výběrové údržby, tak pro soubor s odstraněnými body (odlehlé hodnoty).

včetně odlehlých hodnot bez odlehlých hodnot

ζ’ – ζ ζ’ - ζEGM ζ’ – ζ ζ’ - ζEGM

střední hodnota -0,427 -0,419 -0,427 -0,419

směrodatná odchylka 0,070 0,072 0,032 0,036

minimum -0,921 -0,921 -0,538 -0,519

maximum 0,534 0,521 -0,335 -0,326

Tab. 2 Porovnání rozdílů výškových anomálií spočítanách na bodech výběrové údržby. Všechny hodnoty v [m].

Tabulka 2 ukazuje stejné výsledky, jaké byly citovány v kap. 1, tedy že ačkoliv se průběhy kvazigeoidu liší i o 1 dm, směrodatné odchylky, kterými lze přibližně vyjadřit přesnost modelů ζ a ζEGM jsou 3,2 cm pro model ζ (včetně reziduálních hodnot kvazigeoidu), resp. 3,6 cm pro model ζEGM (vypočítaný pouze z EGM08). Z tohoto pohledu tedy přináší lokální model zlepšení směrodatné odchylky pouze o 0,4 cm.

Takto malý rozdíl ve směrodatné odchylce při poměrně významné rozdílnosti obou modelů kvazigeoidu může znamenat pouze to, že na některých bodech výběrové údržby je přesnější model ζ a na jiných model ζEGM . Následovat tedy bude analýza, na kterých bodech je přesnější který model, případně proč. Na obr. 6 je vykreslen model ζEGM a kroužky představující body výběrové údržby. Střed kroužku vyjadřuje souřadnice bodu výběrové údržby a poloměr kroužku je přímo úměrný hodnotě |ζ’ - ζEGM na daném bodu. Na obr. 7 je znázorněn model ζ a kroužky mají obdobný význam, nyní je jejich velikost ovšem úměrná hodnotě |ζ’ – ζ|. Z těchto obrázků je zřejmé, že v některých místech má model ζ skutečně chyby menší, a tedy došlo ke zpřesnění globálního modelu, na jiných místech má však chyby větší a doplnění informace z lokálních tíhových dat vedlo ve skutečnosti ke snížení přesnosti globálního modelu kvazigeoidu.

Ještě lépe je tento výsledek patrný z obr. 8. Zde jsou zeleně vyznačeny body, kde je model ζ významně přesnější než model ζEGM a červeně vyznačeny body, kde je model ζ významně méně přesný než model ζEGM Hranice významně vyšší resp. nižší přesnosti byla zvolena jako 2 cm, což přibližně odpovídá předpokládané přesnosti hodnoty ζ’ na bodech výběrové údržby. Body, na nichž se hodnoty ζ a ζEGM liší o méně než 2 cm tedy nejsou na obr. 8 znázorněny. Pozadí obr. 8 znázorňuje ζres, aby bylo patrné, jak zlepšení či zhoršení souvisí s reziduálním kvazigeoidem. Poloměr znázorněných kroužků je opět přímo úměrný zlepšení, resp. zhoršení přesnosti vyjádřené snížením resp. zvýšením velikosti rozdílu |ζ’ – ζ| oproti rozdílu |ζ’ - ζEGM |v cm, neboli vyjadřuje velikost |ζ - ζEGM| .

Obr. 6 Výšková anomálie z EGM08 ζEGM [m] a body výběrové údržby.

Obr. 7 Výšková anomálie ζ (kombinace ζEGM a ζres) v [m] a body výběrové údržby.

Obr. 8 Reziduální výšková anomálie ζres [m] a body výběrové údržby.

4 Závěr Odpověď na otázku, zda jsou na území ČR významné rozdíly mezi modelem kvazigeoidu vypočítaným pouze z EGM08 a modelem zahrnujícím rovněž složku ζres vypočítanou z podrobnějších pozemních dat, je kladná. Nejednoznačná odpověď je však na otázku, který z modelů je přesnější. Na některých bodech výběrové údržby je to globální model ζEGM na jiných model ζ. Na 66 bodech došlo k významnému zlepšení a naopak na 43 bodech došlo ke zhoršení přesnosti modelu ζ vůči modelu ζEGM . Zdá se však, že většina bodů, na kterých došlo ke zhoršení přesnosti, leží v oblasti o velikosti 1˚×1˚ přibližně kolem Pardubic, případně pak v pohraničních oblastech, kde je přesnost modelu ζ limitována nedostatkem přesných tíhových dat za hranicemi ČR. Lze se domnívat, že zhoršení je způsobeno nepřesností tíhových nebo výškových lokálních dat, které způsobuje systematický lokální posun reziduálního kvazigeoidu až o několik cm.

Závěrem lze říci, že současný globální model kvazigeoidu EGM08 umožňuje na území ČR vypočítat kvazigeoid se střední chybou přibližně 3,6 cm. Tuto chybu je možné dále snížit pomocí lokálního modelu pozemních tíhových dat, ovšem na některých místech ČR je tento model nepřesný. Tyto lokality ale lze kombinací globálního modelu a bodů výběrové údržby odhalit. Možnost odstranění předpokládaného systematického posunu v těchto lokalitách tak, aby i zde došlo ke zlepšení, nebo aspoň nedošlo ke zhoršení, bude předmětem dalšího studia.

5 Literatura [1] PAVLIS N. K., HOLMES S. A., KENYON S. C., FACTOR J. K. An Earth Gravitational

model to degree 2160: EGM2008. Preseneted at 2008 General Assembly of European Geosciences Union, Vienna, Austria, April 13-18, 2008.

[2] NOVÁK P. Evaluation of local gravity field parameters from high resolution gravity and elevation data. Contributions to Geophysics and Geodesy 36/2006: 1-33.

[3] KOSTELECKÝ J., KOSTELECKÝ J. Jr., PEŠEK I., ŠIMEK J., ŠVÁBENSKÝ O., WEIGEL J., ZEMAN A. Quasigeoids for the territory of the Czech Republic. Studia Geophysica et Geodaetica, č. 48, 3/2004, s. 503-518. ISSN 0039-3169.

[4] NOVÁK P, KLOKOČNÍK J., KOSTELECKÝ J. Testing the global geopotential model EGM07a with GPS/leveling data over the territory of the Czech Republic. Report for the IAG Special Working Group 2.2: Evaluation of Global Earth Gravity Models. 2007.

[5] KADLEC M., KOSTELECKÝ J, NOVÁK P. Database for evaluation of gravity field parameters in Central Europe. Geodeticky a kartograficky obzor 12/2008: 282-288.

[6] KADLEC M. Odstranění systematické složky chyby modelu kvazigeoidu a odhad chyby transformace výšek. Zpráva k řešení projektu InGeoCalc, VÚGTK, 2007.

[7] HEISKANEN W. A., MORITZ H. Physical geodesy. Freeman and Co., San Francisco 1967.