Software pro interpolaci hodnot slunečního záření

ČESKÁ ZEMĚDĚLSKÁ UNIVERZITA V PRAZE Fakulta životního prostředí Katedra aplikované geoinformatiky a územního plánování

Software pro interpolaci hodnot slunečního záření Diplomová práce

Jan Hnilica

Vedoucí diplomové práce: RNDr. Vladimír Puš, CSc.

2010

Poděkování: Rád bych na tomto místě poděkoval vedoucímu diplomové práce RNDr. Vladimíru Pušovi CSc. za cenné rady a trpělivý přístup a konzultantovi Ing. Vojtěchovi Bartákovi za věcné připomínky k problematice. Dále bych rád poděkoval RNDr. Petru Mayerovi za přínosné konzultace numerických problémů. Nakonec bych chtěl poděkovat Soně za trpělivost a pochopení.

Prohlášení: Prohlašuji, že jsem diplomovou práci s názvem „Software pro interpolaci hodnot slunečního záření“ vypracoval samostatně pod vedením RNDr. Vladimíra Puše CSc. a s využitím uvedené literatury. V Praze dne 29. dubna 2010

……………………………….

ABSTRAKT Diplomová práce měla dva cíle. Prvním cílem bylo naprogramovat software umožňující interpolace odhadů globální radiace metodou kriging. Druhým cílem bylo zjistit, zda kriging může poskytovat kvalitní odhady globální radiace, pokud samotnými vstupními daty jsou odhady této veličiny z náhradních meteorologických prvků. Nejprve byla provedena literární rešerše uvedených problematik, na jejímž základě byl vytvořen původní software. Ten je koncipován do dvou oddělených částí. První část umožňuje vytvářet časové řady odhadů globální radiace a může nalézt uplatnění při tvorbě vstupních dat pro hydrologické modely. Druhá část implementuje metodu kriging ve dvou variantách (ordinary kriging a universal kriging), obsahuje široký výběr modelů variogramu a rovněž nástroje pro analýzu prostorových dat. Praktickým použitím softwaru na reálných meteorologických datech bylo zjištěno, že interpolace radiačních odhadů metodou kriging mohou poskytovat věrohodné údaje o globální radiaci. Klíčová slova: odhad globální radiace, kriging, variogram

ABSTRACT The aim of this diploma thesis was to programme an author´s software tool that would be able to calculate spatial interpolations of estimated global radiation data through kriging method. The developed software is divided into two modules. The first one computes time series of estimated global radiation applying five different methods using different meteorological quantities. The second module implements spatial interpolation method kriging in two variants (ordinary kriging and universal kriging) and includes tools for exploratory data analysis. The software was used for interpolation of real meteorological data. This study resulted into the finding that kriging of the estimated radiation data can produce plausible values of global radiation. Key words: global radiation, kriging, variogram

OBSAH 1. ÚVOD ...................................................................................................................... 6 2. METODIKA ............................................................................................................ 8 2.1 Metody odhadu radiace ...................................................................................... 8 2.1.1 Extraterestrická radiace............................................................................. 10 2.1.2 metoda Angstromova - Prescottova .......................................................... 11 2.1.3 Metoda Klabzubova .................................................................................. 13 2.1.4 Metoda Hargreavesova ............................................................................. 13 2.1.5 Metoda Supitova ....................................................................................... 14 2.1.6 Metoda Winslowova ................................................................................. 16 2.2 Interpolační metoda kriging ............................................................................. 19 2.2.1 Vyjádření prostorové závislosti................................................................. 19 2.2.1.1 Variogram .......................................................................................... 20 2.2.1.2 Odhady variogramu............................................................................ 21 2.2.1.3 Modely variogramu............................................................................ 23 2.2.2 Analýza vstupního datového souboru ....................................................... 27 2.2.2.1 Detekce nestacionarity ....................................................................... 27 2.2.2.2 Izotropie ............................................................................................. 28 2.2.2.3 Trend .................................................................................................. 29 2.2.3 Výpočetní schéma krigingu....................................................................... 30 2.2.3.1 Ordinary kriging................................................................................. 32 2.2.3.2 Universal kriging................................................................................ 34 3. POPIS AUTORSKÉHO SOFTWARE .................................................................. 38 3.1 Odhady radiace................................................................................................. 40 3.2 Bodové pole ..................................................................................................... 42 3.3 Kriging ............................................................................................................. 43 3.3.1 Odhady variogramu a analýza dat............................................................. 44 3.3.2 Prokládání odhadů variogramu ................................................................. 46 3.3.3 Kriging ...................................................................................................... 51 3.3.3.1 Ordinary kriging................................................................................. 51 3.3.3.2 Universal kriging................................................................................ 53 4. VÝSLEDKY – INTERPOLACE RADIAČNÍCH ODHADŮ .............................. 55 4.1 Porovnání odhadů radiace ................................................................................ 55 4.2 Porovnání variant krigingu............................................................................... 58 5. DISKUZE............................................................................................................... 61 6. ZÁVĚR .................................................................................................................. 63 Seznam literatury ....................................................................................................... 64 Seznam obrázků ......................................................................................................... 67

5

1. ÚVOD Primárním cílem předložené diplomové práce bylo vytvoření původního autorského software na bázi geografických informačních systémů, umožňujícího statistickou interpolaci bodových odhadů slunečního záření. Hodnoty dopadajícího slunečního záření (tzv. globální radiace) na zájmové území jsou jedním ze zásadních faktorů, majících vliv na redistribuci vody v systému a jako takové jsou nepostradatelnou vstupní veličinou do mnoha modelů, ve kterých je řešena otázka odhadu potenciální nebo aktuální evapotranspirace. To se týká jednak některých modelů hydrologických, dále například modelů půdní vlhkosti nebo růstových modelů plodin (Trnka et al. 2005). Problémem, který je v odborné literatuře často zmiňován, např.(Liu, Scott 2001), je nízká dostupnost měřených dat globální radiace. Ta vyplývá z toho, že měření této veličiny vyžaduje relativně náročnou instalaci a obsluhu nákladných zařízení, jakými jsou aktinometry či pyranometry, přičemž většina klimatických a meteorologických stanic není těmito přístroji vybavena. Jak uvádí Trnka et al. (2005), podíl počtu stanic, registrujících denní sumy globální radiace ku počtu stanic, provádějících pouze měření srážek a teplot, se v České Republice a Rakousku pohybuje kolem 1:20, v globálním měřítku je tento podíl podle Thornton el al. (1997) dokonce 1:500. Z výše uvedených důvodů bylo vyvinuto množství postupů, umožňujících hodnoty globální radiace v zájmovém bodě odhadovat. Jejich výčet zahrnuje odhady pomocí metod dálkového průzkumu Země, substituce s využitím dat z blízkých stanic, využití stochastických generátorů počasí, interpolace měřených hodnot, aplikace neuronových sítí a také empirické vztahy, umožňující odhadovat globální radiaci z jiných meteorologických prvků, které jsou obvykle běžně (nebo alespoň častěji) dostupné, jako jsou například hodnoty teplot, srážek či oblačnosti (Trnka et al. 2005). Posledně jmenovaná skupina, tedy metody odhadu z jiných meteorologických prvků, je předmětem zájmu této diplomové práce. Jak již bylo uvedeno výše, výstupem je software umožňující interpolaci těchto odhadů. Výsledkem jeho spuštění je (v optimálním případě) rastr prostorového rozložení hodnot globální radiace. Po konzultaci s vedoucím diplomové práce byl z mnoha možných interpolačních metod vybrán kriging.

6

Od doby, kdy byla tato metoda poprvé popsána (Matheron 1963), se těší značnému zájmu odborné veřejnosti a prodělala dlouhý vývoj. To je důvodem skutečnosti, že termín kriging v sobě dnes zahrnuje celou škálu interpolačních technik, jejichž společným jmenovatelem je snaha o optimální lineární predikci (Cressie 1993). Diplomová práce je rozdělena do tří základních částí. V první, metodické části, jsou popsány jednak metody odhadu globální radiace, které byly pro účely diplomové práce vybrány, dále jsou popsány principy a jednotlivé techniky interpolační metody kriging. Druhá část je zaměřena na představení vlastního software,

čtenář

je

seznámen

s celkovou

architekturou

i

konkrétním

programátorským řešením jednotlivých problémů. Ve třetí části je předloženo a diskutováno praktické použití software.Především je řešena otázka, zda (a pokud ano, tak při jakém nastavení parametrů) je kriging schopen poskytovat kvalitní odhady globální radiace v situaci, kdy samotnými vstupními hodnotami jsou odhady této veličiny.

7

2. METODIKA V této části jsou nejprve představeny jednotlivé metody odhadu globální radiace z náhradních meteorologických prvků. Je zdůvodněn jejich výběr a popsán jejich algoritmus. Dále je řešena otázka jejich aplikovatelnosti vzhledem k náročnosti na vstupní data a také vhodnost výběru jednotlivých metod vzhledem ohledem na klimatické podmínky zájmové oblasti. Druhý oddíl je věnován interpolační metodě kriging. Jsou vysvětleny základní principy prostorové statistky, které jsou s krigingem neodmyslitelně spjaty, dále je popsáno schéma jednotlivých technik krigingu. Metodická část nepojednává o problematice převedení těchto metod do softwarové podoby. Specifické problémy, které z tohoto úkolu plynou, jsou popsány v další části věnované představení autorského softwaru.

2.1 Metody odhadu radiace Před vlastním rozborem metod odhadu radiace je vhodné stručně definovat pojmy, které budou v následujících odstavcích používány. Sluneční záření je souhrnným pojmem, zahrnujícím v sobě širokou škálu částic a fotonů, dopadajících na Zemi v celém spektru vlnových délek od krátkých vln kosmického záření až po vlny velmi dlouhé. Po zanedbání okrajových částí tohoto spektra, které tvoří pouhé jedno procento energie, dostáváme záření o vlnových délkách od 170 do 4000 nanometrů, které je v meteorologické literatuře označováno pojmem krátkovlnné záření (Klabzuba 2001). Během překonávání dráhy mezi Sluncem a Zemí je krátkovlnné záření nepatrně ovlivňováno interakcemi s částicemi v meziplanetárním prostoru, ale toto ovlivnění je nevýznamné. Intenzita slunečního záření na horní hranici atmosféry je proměnlivá pouze díky měnící se vzdálenosti Země a Slunce, která je způsobena eliptickou trajektorií při oběhu Země. Průměrná intenzita slunečního záření na horní hranici atmosféry, vztažená ke střední vzdálenosti Země od Slunce, je nazývána solární konstantou, a její hodnota je stanovena na 1367 W ⋅ m −2 (Sobíšek 1993). Při průchodu atmosférou je část záření odražena a rozptýlena na prachových částečkách, molekulách plynů a kapičkách vody. Záření se proto na zemský povrch dostává ve dvou formách: jednak jako tzv. přímé záření tvořené (vzhledem k velké 8

vzdálenosti Slunce a Země) prakticky rovnoběžnými paprsky, a jednak jako tzv. difúzní záření, dopadající na povrch Země ze všech směrů po předchozím rozptylu a odrazu v atmosféře. Součet přímého a difúzního záření se nazývá globální záření (globální radiace) a představuje energetickou dotaci zemskému povrchu ve formě slunečního záření (Klabzuba 1991). Toto záření je předmětem této diplomové práce. Některé meteorologické prvky, běžně registrované na většině meteorologických stanic, např. teplotní extrémy, oblačnost, trvání slunečního svitu atp., jsou s množstvím globálního záření více či méně těsně spjaty. Jak již bylo řečeno v úvodu, předmětem této práce jsou empirické vztahy, které tuto závislost využívají a dovolují odhadovat denní sumu globální radiace z měřených hodnot těchto veličin. Za tímto účelem bylo různými autory vyvinuto značné množství metod, které se vzájemně liší jednak požadovanými vstupními údaji (tedy volbou meteorologických prvků, ze kterých je radiace odhadována), dále složitostí výpočtu, aplikovatelností s ohledem na klimatické podmínky zájmové oblasti a v neposlední řadě také přesností výsledných odhadů. Lepší výsledky je možno očekávat při použití prvků, které jsou s globálním zářením v přímém vztahu (obecně nejlepší veličinou je z tohoto pohledu trvání slunečního svitu), naopak prvky, jako např. teplotní extrémy, jsou sice běžněji dostupné, ale jejich hodnoty jsou ovlivněny dalšími faktory, takže výsledný odhad může být od skutečnosti výrazněji odchýlen (Trnka et al. 2005). Metody pro účely této diplomové práce byly vybrány na základě dvou kriterií. Zaprvé bylo požadováno, aby vstupní údaje jednotlivých metod pokrývaly jako celek co možná nejširší škálu běžně měřených meteorologických prvků. Druhým kriteriem byla použitelnost metod pro klimatické podmínky České republiky, potažmo středoevropského regionu. Po zvážení těchto kriterií bylo vybráno následujících pět metod:


Angstromova - Prescottova




Klabzubova




Hargreavesova




Supitova




Winslowova

(Martínez - Lozano et al. 1994) (Klabzuba et al. 1999) (Hargreaves et al. 1985) (Supit et van Kappel 1998) (Winslow et al. 2001)

S výjimkou metody Klabzubovy vyžadují všechny ostatní metody jako jeden ze vstupních údajů hodnotu extraterestrické radiace. 9

2.1.1 Extraterestrická radiace Jak již bylo uvedeno výše, je hodnota intenzity slunečního záření dopadajícího na vnější povrch atmosféry proměnlivá v důsledku měnící se vzdálenosti Země od Slunce v průběhu roku, přičemž střední hodnotu této intenzity udává solární konstanta. Aktuální intenzitu záření pro daný den v roce a dané místo na Zemi je tedy možno spočítat jako hodnotu solární konstanty násobenou bezrozměrným koeficientem, reflektujícím polohu bodu a čas. Tato aktuální intenzita je nazývána extraterestrickou radiací. Pro její výpočet byl použit způsob, který uvádí Allen et al. (1998):

Re =

1440

π

⋅ Ro ⋅ dr ⋅ [ω ⋅ sin (ϕ ) ⋅ sin (δ ) + cos(ϕ ) ⋅ cos(δ ) ⋅ sin (ω )]




Re – extraterestrická radiace




Ro – solární konstanta




dr – inverzní vzdálenost Země – Slunce





 2π  dr = 1 + 0.033 ⋅ cos ⋅J  365  J – číslo dne (1 - 366) ω – hodinový úhel západu slunce [rad]




(2.1)

ω = arccos[− tg (ϕ ) ⋅ tg (δ )]

ϕ – zeměpisná šířka bodu [rad]

δ – solární deklinace [rad]  2π  δ = 0.409 ⋅ sin ⋅ J − 1.39   365  Obecně platí, že radiace dopadající na zemský povrch je pouze částí radiace


extraterestrické. Záření je při průchodu atmosférou zeslabeno, což je způsobeno interakcemi s částicemi v atmosféře. Tyto interakce jsou velmi složité, souhrnně se vliv složek atmosféry na oslabení dopadající radiace vyjadřuje tzv. atmosférickou transmitancí (propustností), což je bezrozměrný koeficient z intervalu <0,1>, kterým se násobí hodnota extraterestrické radiace, čímž je obdržena aktuální hodnota dopadajícího záření (Sobíšek 1993). Některé z níže popsaných metod se soustředí právě na vyčíslení hodnoty transmitance. V dalších oddílech této části následuje popis jednotlivých metod odhadu globální radiace.

10

2.1.2 metoda Angstromova - Prescottova Metoda je modifikací původní práce Angstromovy (Martínez - Lozano et al. 1994 ex Angstrom 1924). Ten se zaměřil na dobu trvání slunečního svitu jako na nejvhodnější meteorologickou charakteristiku pro účely odhadu globální radiace. Původní formule pro odhad byla tohoto tvaru: Rcelk n =a +b⋅ Rmax N

(2.2)


Rcelk – celková suma globální radiace
Rmax – viz níže
n – skutečné trvání slunečního svitu v daném dni
N – maximální astronomicky možné trvání slunečního svitu v daném dni
a, b – empirické koeficienty. Vztah (2.2) staví do souvislosti hodnotu relativního trvání slunečního svitu, vyjádřeného jako poměr skutečné a maximální možné hodnoty této veličiny v daném dni a relativního množství příchozí globální radiace. Člen Rmax představuje v rovnici hypotetické maximum příchozí radiace a může být reprezentován různými způsoby. Nejčastěji je za Rmax dosazována hodnota extraterestrické radiace ( Re ), ale v literatuře je možno nalézt i práce, uvažující hypotetickou hodnotu radiace pod reálnou atmosférou za zcela bezmračného dne ( Rc ). S volbou Rmax samozřejmě souvisí odlišný náhled na hodnotu a význam koeficientů a, b. Použití Rc (např. Angstrom 1956) má nevýhodu v tom, že pro přibližné vyčíslení Rc musí být na zájmovém území k dispozici měřená radiační data, a proto je pro

lokality bez jakýchkoliv aktinometrických záznamů (i v minulosti) tato varianta nepoužitelná. Pro účely této práce byla použita varianta s extraterestrickou radiací. Poprvé byla tato modifikace popsána Prescottem (1940). Rovnice odhadu Rcelk má tvar: n  Rcelk = Re ⋅  a + b ⋅  N 

(2.3)

Při takto definovaném vztahu a použití extraterestrické radiace jako hypotetického maxima je možno interpretovat fyzikální význam empirických koeficientů a, b takto:

11




a je minimální podíl extraterestrického záření, která proniká na Zemi v podmínkách, kdy je obloha zcela pokryta mraky,




součet (a+b) je podíl extraterestrického záření, pronikající na Zemi v podmínkách maximální propustnosti atmosféry, tedy za zcela bezoblačné oblohy. Při aplikaci vztahu (2.3) vyvstává problém určení hodnot koeficientů a, b. Ze

stručného přehledu prací, které se touto problematikou zabývaly (viz Martínez Lozano et al. 1994) vyplývá, že mezi faktory, které hodnoty koeficientů ovlivňují patří poloha lokality, data k odhadu použitá (hodinové, denní či měsíční záznamy) a také roční období. Nejvýraznější vliv ze jmenovaných faktorů má jednoznačně poloha lokality, resp. její zeměpisná šířka. Bylo učiněno několik pokusů o nalezení funkční závislosti pro hodnoty koeficientů a, b na zeměpisné šířce, např. Glover et McCulloch (1985), ovšem s problematickými výsledky. Obdobná situace je z hlediska analýzy vlivu ročního období. Stanovení koeficientů pro konkrétní oblast je tedy nutno provádět analýzou lokálních meteorologických dat. Existují také některé modifikace Angstromovy - Prescottovy metody, které vycházejí např. z faktu, že změny v poměru n N , které nastávají v poledních hodinách, mají na celkovou denní hodnotu radiace mnohem větší vliv než obdobné výkyvy v ranních či večerních hodinách, což vede k myšlence váženého průměrování. Při tomto přístupu je ovšem nutno disponovat hodinovými záznamy slunečního svitu. Přehled příslušných modifikací s odkazy na literaturu nabízí Martínez - Lozano et al. (1994). Pro účely diplomové práce byly odhady prováděny podle vztahu (2.3), přičemž hodnoty koeficientů pro odhady denních sum globální radiace byly interpolovány z map pro střední Evropu vytvořených v rámci projektu MARS (Monitoring Agriculture with Remote Sensing) organizovaném Joint Research Centre of EU. Mapy jsou umístěny na internetové adrese http://supit.net. Využití doby trvání slunečního svitu pro odhady globální radiace s sebou nese několik nevýhod. Doba trvání např. nevypovídá o intenzitě záření, dále díky zpětnému odrazu od oblaků může nastat vyšší hodnota radiace při částečně zatažené obloze, což je proti smyslu výpočetního schématu. Ovšem z meteorologických prvků 12

má doba svitu s globální radiací nejužší vztah a metoda většinou poskytuje velmi kvalitní odhady (Trnka et al. 1999).

2.1.3 Metoda Klabzubova Také tato metoda vychází z doby trvání slunečního svitu, nicméně oproti metodě Angstromově, u které lze podat přibližnou fyzikální interpretaci empirických koeficientů, nemá její výpočetní schéma žádný fyzikální základ. Bylo odvozeno jako regresní matematický vztah mezi relativní dobou slunečního svitu pro určitý den v roce a příslušnou sumou globální radiace. K tomu byla použita data ze solární observatoře v Hradci Králové. Tvar výpočtové rovnice je následující:

Rcelk = 7.19 + 0.2508 ⋅ S r − 9.28 ⋅ 10 −6 ⋅ (S r + 22.9) ⋅ ( J − 147.7 )




2

(2.4)

R celk – celková suma globální radiace n Sr = – relativní trvání slunečního svitu (viz výše) N J – pořadové číslo dne v roce (1 - 366)

Fakt, že celý vztah je založen na analýze dat z jedné stanice, naznačuje jeho omezenou využitelnost pro odlišné klimatické podmínky. Rovněž tak je problematické jeho využití v chladné části roku, neboť v tomto období jsou i při velkých hodnotách relativního svitu dosahovány malé radiační úhrny. Model byl totiž primárně koncipován pro aplikaci v růstových modelech plodin ve vegetačním období a při použití v zimních měsících může algoritmus vést i k záporným hodnotám radiace (Klabzuba et al. 1999). Nicméně pro teplejší části roku na území střední Evropy je model dobře využitelný, ve srovnání s Angstromovým modelem navíc odpadá nutnost zjišťovat hodnoty empirických koeficientů pro danou lokalitu (Trnka et al. 1999).

2.1.4 Metoda Hargreavesova V případě, že na zájmové lokalitě jsou k dispozici pouze záznamy denních teplotních extrémů, je hojně využívaným nástrojem odhadu globální radiace Hargreavesova metoda. Vychází z práce Hargreavese et al. (1985), kdy primárním cílem autorů bylo vyvinout metodiku závlahového režimu pro aridní oblasti západní Afriky (Senegal). Za tím účelem byly testovány různé vztahy pro výpočet referenční evapotranspirace, která je kritickým údajem, nutným k úvahám o zavlažovacím režimu. Hodnota evapotranspirace je v úzkém vztahu se sumou globální radiace, bez

13

které je její výpočet nemožný. V podmínkách rozvojových zemí je ale nutno uskrovnit požadavky na meteorologická data, tudíž proto byl hledán vztah umožňující odhady radiace z co možná nejdostupnějších hodnot (což jsou právě hodnoty teplotního maxima a minima). Původní vztah pro odhady radiace má tento tvar: Rcelk = Re ⋅ 0.16 ⋅ Tmax − Tmin


(2.5)

Tmax ,Tmin – denní teplotní maximum a minimum

Kalibrace vztahu probíhala pouze na lokalitách s velmi specifickým klimatem západní Afriky. Aby bylo umožněno použití i v odlišných klimatických podmínkách, byl vztah (2.5) upraven na tvar: Rcelk = Re ⋅ a ⋅ Tmax − Tmin + b


(2.6)

a, b – empirické koeficienty.

Přestože je metoda původně určena pro subtropické oblasti a její použitelnost v jiných podmínkách je i přes uvedenou úpravu omezená, je v praxi využívána velmi

často, a to zejména pro svoji jednoduchost a nenáročnost na vstupní data (Trnka et al. 1999). Hodnoty empirických koeficientů jsou pro podmínky střední Evropy relativně stálé a pro území České Republiky mohou být považovány za konstanty. Lze je interpolovat z map a tabulek umístěných na http://supit.net nebo jsou k dispozici v Trnka et al.(1999).

2.1.5 Metoda Supitova Metoda vznikla pro potřeby hodnocení vlivu počasí na produkci zemědělských plodin na území Evropy, prováděném systémem CGMS (Crop Growing Monitoring System). Tato aktivita byla jednou ze součástí projektu MARS, který byl již zmiňován výše, podrobnější informace viz např. Supit (1997). V rámci CGMS byl s pomocí denních hodnot meteorologických veličin, půdních charakteristik a údajů o způsobu zemědělského managementu simulován potenciální růst plodin. CGMS operuje v prostorovém rozlišení 50 x 50 km nad celým evropským územím a globální radiace je jedním z požadovaných vstupů do systému. Supit et van Kappel (1998) hledali optimální metodu, která by poskytovala kvalitní odhady

14

globální radiace pro každý čtverec. Byly uvažovány možnosti prostorové interpolace skutečných radiačních měření, ale síť meteorologických stanic ve kterých jsou měření prováděna, byla pro tyto účely příliš řídká, neboť např. podle Hubbarda (1994) je pro zachycení prostorové variability globální radiace nutno dodržet vzdálenost odhadované a měřené lokality na maximální hranici 30 km. Rovněž nebyly využity odhady pomocí satelitních dat, které poskytují sice prostorově rozložené, ale kvalitativně diskutabilní výsledky. Z výše uvedených důvodů byla zkoumána možnost využít metod odhadu z náhradních meteorogických prvků. Po vyloučení několika metod, z různých důvodů nevhodných pro specifický úkol plošných odhadů nad celou Evropou (Hargreaves, Angstrom-Prescott), byla vyvinuta metoda nová, která vyžaduje jako vstupní údaje hodnoty denních teplotních extrémů a hodnotu celkové oblačnosti. Teplotní extrémy jsou běžně dostupné a použití oblačnosti skýtá výhodu v tom, že ačkoliv není příliš hustě měřeným meteorologickým prvkem, může být odečítána ze satelitních snímků. Dříve publikované analýzy globální radiace a oblačnosti naznačují nelineární vztah mezi těmito veličinami (Worner 1967). Výsledná rovnice odhadu podle Supit et van Kappel (1998) je kombinací rovnic Wornera (1967) a Hargreavese et al. (1985) a má tento tvar:   O  Rcelk = Re ⋅ a ⋅ Tmax − Tmin + b ⋅ 1 −   + c 8    




O – celková denní oblačnost [1/8]




a, b, c – empirické koeficienty

(2.7)

Hodnoty empirických koeficientů byly kalibrovány a validovány s využitím dat z množství meteorologických stanic po celé Evropě. Jejich hodnoty je opět možno získat na http://supit.net, popřípadě jsou uvedeny v práci Supit et van Kappel (1998). V případě, že jsou k dispozici záznamy oblačnosti, je tato metoda kvalitním nástrojem k odhadu globální radiace, a přestože byla primárně koncipována pro evropské podmínky, je v současnosti využívána v systémech predikce růstu plodin v mnoha částech světa (Trnka et al. 1999).

15

2.1.6 Metoda Winslowova Tato metoda je založena na detailním rozpracování transmitance atmosféry. Motivací při jejím návrhu byla práce Bristow et Campbell (1984), která se při odhadu radiace rovněž soustředí na hodnotu transmitance. Její výpočet však vyžaduje znalost tří místně specifických empirických koeficientů, které nemají jasnou fyzikální interpretaci a které je nutno odvodit z dlouhodobých klimatologických pozorování na dané lokalitě. To činí postup nepraktickým pro běžné použití. Winslow et al. (2001) se zaměřili na odstranění tohoto nedostatku s cílem vytvořit model dobře použitelný při různých klimatických podmínkách včetně tropických regionů. Základem je empirický inverzní vztah mezi transmitancí a relativní vlhkostí vzduchu. Relativní vlhkost vzduchu je poměr mezi atmosférickým tlakem vodních par a tlakem nasycených vodních par a může být odhadnuta z denní maximální a minimální teploty. Příchozí radiace během dne způsobuje zvýšení teploty vzduchu, což má za následek pokles relativní vlhkosti. Během nárůstu teploty z denního minima na maximum je množství radiace, které tento nárůst způsobilo, úměrné odpovídajícímu poklesu relativní vlhkosti:

R T min − T max ≈ rh T min − rh T max


R T min −T max – příchozí radiace během nárůstu teploty z minima na

maximum


rh T – relativní vlhkost vzduchu při teplotě T.

Zpravidla je možno klást rh T min = 1 , uvedený vztah může být tedy zjednodušen na: R T min −T max ≈ 1 − rh T max Hodnota R T min −T max reprezentuje větší část denního globálního záření, nikoli celkovou sumu. Pokud je zavedena funkce D předpisem Rcelk / RT min −T max , potom pro celkovou globální radiaci dostáváme vztah:

(

Rcelk ≈ D ⋅ 1 − rh T max

)

Uvedené teoretické úvahy jsou výchozími momenty, na kterých je založeno odvození konečného výpočetního vztahu Winslowovy metody. Metoda počítá též s

16

dalšími faktory zajišťujícími univerzálnější použitelnost pro různé klimatické podmínky (původní úvahy byly prováděny nad klimatickými daty z velmi teplých humidních oblastí). Konečný algoritmus používá jako vstupní údaje tyto veličiny:
denní maximální a minimální teplota vzduchu
denní úhrn srážek
průměrná roční teplota dané lokality
zeměpisná šířka
nadmořská výška. Výpočetní rovnice má pak tento tvar:  e (T min)  Rcelk = Re ⋅ τ cf ⋅ 1 − β ⋅ s ⋅D e s (T max)  

(2.8)

Význam a výpočet jednotlivých členů je následující: o τ cf = (To ⋅ Ta ⋅ Tv ) Po – celková transmitance [-] P




(

)

To = 0.947 − 1.033 ⋅ 10 −5 ⋅ Φ -

2.22

– transmitance čistého suchého vzduchu [-]

Φ – zeměpisná šířka [°]




Ta – transmitance způsobená adsorpcí záření aerosoly (dosazuje se 1)




Tv = 0.9636 − 9.092 ⋅ 10 −5 ⋅ (Tmean + 30)

1.8232

– transmitance ovlivněná

přítomností vodní páry ve vzduchu [-] -

Tmean – průměrná roční teplota lokality [°C] (pokud denní srážky překročí 1.0 mm, redukuje se Tv o 0.13)




P = (1 − 2.2569 ⋅ 10 −5 ⋅ z ) 5.2553 – korekce tlaku na nadmořskou výšku [-] Po

-

z – nadmořská výška [m.n.m.]

 m ⋅T  o e s (T ) = 0.611 ⋅ exp  – tlak nasycených vodních par při teplotě T [kPa] n +T 




m, n – empirické koeficienty [-]

2   π   H −   4  o D = 1 −   2H 2     

−1

– korekce chyby způsobené různou délkou dne

17




H = arccos(− tgΦ ⋅ tgδ ) – „half day length“ [-], viz Sellers (1965) -

δ – solární deklinace [rad], výpočet viz Allen (1998)

  ∆T  o β = MAX 1.041;23.753 ⋅ + 273 . 16 T mean  


β – empirický koeficient [-]




∆T – průměrný roční rozsah teplot [°C]

(2.9)

o Re – denní suma extraterestrické radiace, výpočet viz Allen (1998). Pokud vypočtená hodnota globální radiace klesne pod hodnotu 0.1 ⋅ Re , doporučuje autor metody výsledek upravit nastavením na tuto hodnotu. Podle zkušeností je i při silně nepropustné atmosféře vždy přítomen dostatek difúzní radiace potřebný k tomu, aby byla dosažena alespoň hodnota jedné desetiny extraterestrické radiace. Koeficient β je ve většině případů roven 1.041, korekce podle vztahu (2.9) se provádí pouze pro horské lokality s velkým denním rozsahem teplot (Winslow et al. 2001). Z výše uvedeného popisu je patrné, že se jedná o poměrně komplikovaný algoritmus, ovšem použití Winslowovy metody má několik výrazných pozitiv. Jednak jsou vstupními hodnotami denní teplotní extrémy a srážkový úhrn, což jsou veličiny běžně dostupné, dále není potřeba zjišťovat žádné empirické koeficienty specifické pro danou lokalitu. Navíc podle Trnka et al. (2005) má metoda univerzální použitelnost pro všechny zeměpisné šířky.

Druhý díl metodické části je věnován stručnému popisu teoretických principů a výpočetních schémat interpolační metody kriging.

18

2.2 Interpolační metoda kriging Cílem interpolace je odhad neznámé hodnoty náhodného prostorového procesu v určitém bodě na základě pozorovaných hodnot v okolních bodech. Procesem může být např. sluneční záření nebo déšť, hodnotou pak suma globální radiace či srážkový úhrn v nějakém místě. Metoda kriging byla pojmenována po jihoafrickém důlním inženýrovi D. G. Krigeovi, který vyvinul empirické metody pro odhady ložisek rud založené na prostorové analýze vzorkovacích vrtů. Původní matematická formulace pochází z práce Matherona (1963), který také metodu po Kriegovi pojmenoval. Ovšem jak podotýká Cressie (1993), myšlenky velmi podobné Matheronovým lze ve stejném období nalézt u autorů ze Sovětského Svazu (Gandin 1963), kde byly téměř identické postupy vyvinuty v oblasti interpolací meteorologických dat. Podrobnosti o původní práci Kriega i o následném vývoji metody podává Cressie (1990). Kriging vychází ze základního principu prostorové statistiky, kterým je tvrzení, že jevy, které nastávají blíže sobě, mají tendenci být si podobnější než jevy, které nastávají dále od sebe. Mezi jevy v prostoru tedy můžeme pozorovat jistou závislost, která vyplývá ze vzájemné polohy míst, kde k těmto jevům došlo. Při predikci neznámé hodnoty z pozorovaných hodnot v okolních bodech je tato závislost reflektována a pozorování, která se nacházejí blíže odhadovanému bodu mají na odhad větší vliv než pozorování vzdálenější (Kraus 2007). K tomu, aby mohlo být této závislosti využito při predikci, je nutno ji nějakým způsobem popsat a kvantifikovat. Tomuto tématu se věnuje následující oddíl.

2.2.1 Vyjádření prostorové závislosti Před vlastním popisem vyjádření závislosti je nutno definovat studovaný prostorový proces. V celém dalším textu je uvažován dvourozměrný eukleidovský prostor R 2 , kde s ∈ R 2 je konkrétní bod v tomto prostoru (definovaný souřadnicemi x a y). Symbolem Z (s ) je míněna hodnota prostorového procesu Z v bodě s. Vzdáleností mezi dvěma body s1 a s 2 prostoru je míněna eukleidovská vzdálenost: s1 − s 2 =

(s1x − s 2 x )2 + (s1 y − s 2 y )2

19

Hodnotu Z (s ) v libovolném bodě s je možno rozepsat jako součet dvou složek:

Z (s ) = µ (s ) + δ (s )


µ (s ) – trendová složka, střední hodnota procesu




δ (s ) – náhodná složka

(2.10)

Trendová složka je buď konstantní v celém prostoru, nebo obecněji funkcí polohy bodu. Náhodná složka je v každém případě funkcí polohy v prostoru a má nulovou střední hodnotu. Právě v náhodné složce je obsažena prostorová závislost, kterou je nutno popsat. Pro vyjádření prostorové závislosti je možno použít několika statistických veličin, obvyklé je použití rozptylu nebo kovariance. Cressie (1993) doporučuje zvolit první možnost, protože volba rozptylu jako veličiny popisující prostorovou závislost nevyžaduje tak silné předpoklady o studovaném prostorovém procesu jako v případě kovariance, jak bude upřesněno níže. Proto bylo také v této diplomové práci za tímto účelem využito výhradně rozptylu. 2.2.1.1 Variogram Při volbě rozptylu (variance) se veličina popisující prostorovou závislost sledovaného procesu nazývá variogram, značí se symbolem 2γ (⋅) a je definována následujícím způsobem:

2γ (h ) = D[Z (s + h ) − Z (s )]


D – rozptyl




Z (s ) – hodnota realizace prostorového procesu v bodě s




Z (s + h ) – hodnota v bodě (s+h)




2γ (h ) – hodnota variogramu pro posun mezi body délky h.

(2.11)

Variogram tedy popisuje, jak se mění rozptyl rozdílu hodnot mezi dvěma body prostoru v závislosti na měnící se vzdálenosti těchto bodů. Poloviční hodnota variogramu, tj. γ (⋅) je označována jako semivariogram. Při použití kovariance je veličina popisující prostorovou závislost označována jako kovariogram C (⋅) :

C (s1 − s 2 ) = cov[Z (s1 ), Z (s 2 )]

(2.12)

20




cov – kovariance

Aby mohla být prostorová závislost popsána pomocí kovariogramu na základě pozorovaných hodnot v bodech prostoru, je nutno u prostorového procesu předpokládat tzv. stacionaritu druhého řádu. Pokud uvedený předpoklad neplatí, mohou nastat případy, kdy zatímco variogram je definován, kovariogram buďto vůbec neexistuje, nebo je jeho odhad zatížen neúměrně velikou chybou. Podrobnosti viz Cressie (1993). Odhady variogramu z pozorovaných dat v prostoru také vyžadují jisté předpoklady o sledovaném procesu, ty budou přiblíženy v následující kapitole. 2.2.1.2 Odhady variogramu Uvažujme, že byla provedena sada měření hodnot nějakého prostorového procesu Z na n lokalitách, k dispozici je tedy soubor hodnot Z (s1 ), Z (s 2 ),........., Z (s n ) ,

přičemž je předpokládáno jejich normální rozdělení. Konečným cílem je odhadnout neznámou hodnotu Z (s o ) v nějakém bodě so . Jak již bylo uvedeno, pro výpočet pomocí krigingu je nutné kvantifikovat prostorovou závislost obsaženou ve sledovaném procesu, neboli vyjádřit konkrétně jeho variogram. Prvním krokem, který je nutno provést, jsou bodové odhady variogramu na základě hodnot vstupního datového souboru. V literatuře je popsáno několik způsobů odhadů variogramu. Do softwaru vzniklého v rámci této práce byly zahrnuty dva. Prvním je klasický odhad podle Matherona (1963):

1 Nh 2γˆ (h ) = ⋅ ∑ (Z (s i + h ) − Z (s i )) N h i =1


2

(2.13)

N h – počet různých párů bodů, vzdálených od sebe o délku (popř. směr) vektoru h.

Druhým použitým typem odhadu je tzv. robustní odhad podle Cressie et Hawkins (1980), který je odolnější vůči extrémním hodnotám ve vstupním datovém souboru, tzv. „outliers“. Tyto hodnoty mohou vznikat chybou měření, chybným zápisem apod.

21

Robustní odhad má tento tvar:  1 Nh   ⋅ ∑ Z (s i + h ) − Z (s i )  N  2γ (h ) =  h i =1

4

 0.494   0.457 +  N h  

(2.14)

Aby mohl být variogram odhadován ze vstupního datového souboru, je nutno předpokládat tzv. vnitřní stacionaritu procesu. Proces je stacionární, pokud pozorovaná závislost mezi dvěma body prostoru pramení pouze z jejich vzájemné relativní polohy a nikoliv z konkrétního místa, kde se body nacházejí. Matematická formulace vnitřní stacionarity je:

E [Z (s )] = µ

(2.15)

E [Z (s + h ) − Z (s )] = 2γ (h ) 2

(2.16)

Předpoklad tedy spočívá jednak v konstantní střední hodnotě procesu, jednak v tom, že střední hodnota kvadrátů rozdílů závisí pouze na relativní vzájemné poloze bodů a je rovna hodnotě variogramu pro vektor posunu h mezi body. Pokud je navíc hodnota variogramu dána pouze délkou vektoru h a nezáleží na jeho směru, proces je označován jako izotropní (Cressie 1993). Vzhledem k těmto požadavkům je nutno vstupní data podrobit analýze, o které pojednává kapitola 2.2.2. Výsledkem bodových odhadů je soubor hodnot variogramu pro různé vzdálenosti mezi body (Obr. 1).

Obr. 1 Bodové odhady variogramu

22

Pro sestavení rovnic krigingu nestačí znalost hodnot variogramu pro vzdálenosti mezi body, které poskytl vstupní datový soubor, ale je potřeba znát hodnotu pro obecně libovolnou vzdálenost. Z tohoto důvodu je dalším krokem při kvantifikaci prostorové závislosti proložení získaných odhadů teoretickou křivkou. 2.2.1.3 Modely variogramu Proces proložení experimentálních bodů křivkou je předmětem regresní analýzy. Obecný postup při regresních úlohách je obvykle následující: a) formulace tvaru prokládané funkce (tedy sestavení teoretického modelu, který bude známý až na libovolný konečný počet neznámých parametrů), b) formulace účelové funkce (stanovení funkčního předpisu, pomocí něhož bude minimalizována suma reziduálů, tedy rozdílů mezi hodnotami teoretického modelu a hodnotami experimentálních bodů), c) výběr a aplikace metody, kterou je nalezeno minimum účelové funkce vzhledem k hledanému vektoru parametrů modelu (Javůrek et Taufer 2006a). V případě prokládání hodnot variogramu není nutné formulovat vlastní tvar prokládané funkce, v literatuře je popsáno množství teoretických modelů (semi)variogramu. První bod regrese pak spočívá ve výběru z těchto nabízených modelů. Do softwarové podoby bylo v rámci diplomové práce převedeno těchto šest modelů, které popisuje Cressie (1993): •

lineární model

γ (h, λ ) = co + bl ⋅ h

γ (h, λ ) = 0

pro h > 0 pro h = 0

(2.17)

λ = (co , bl ) – vektor parametrů modelu, co ≥ 0 , bl ≥ 0

•

•

exponenciální model   − h   pro h > 0 γ (h, λ ) = co + ce ⋅ 1 − exp a e     γ (h, λ ) = 0 pro h = 0 λ = (co , ce , a e ) co ≥ 0 , ce ≥ 0 , a e ≥ 0

(2.18)

sférický model

3 h 1  h γ (h, λ ) = co + c s ⋅  ⋅ − ⋅   2 a s 2  a s γ (h, λ ) = co + c s

  

3

  

pro 0 < h ≤ a s

(2.19)

pro h > a s

23

γ (h, λ ) = 0 λ = (co , c s , a s ) •

pro h = 0 co ≥ 0 , c s ≥ 0 , a s ≥ 0

kvadratický (“rational quadratic”) model cr ⋅ h 2 γ (h, λ ) = co + pro h > 0 1 + h 2 / ar γ (h, λ ) = 0 pro h = 0 λ = (co , c r , a r ) co ≥ 0 , c r ≥ 0 , a r ≥ 0

(

•

mocninný model γ (h, λ ) = co + c p ⋅ h ap

γ (h, λ ) = 0 λ = (co , c p , a p )

•

)

pro h > 0

(2.20)

(2.21)

pro h = 0 co ≥ 0 , c p ≥ 0 , a p ∈< 0,2)

vlnový (“wave”) model  a ⋅ sin (h / a w )  γ (h, λ ) = co + c w ⋅ 1 − w pro h > 0  h   γ (h, λ ) = 0 pro h = 0 λ = (co , c w , a w ) cw ≥ 0 , cw ≥ 0 , aw ≥ 0

(2.22)

U většiny teoretických modelů je možné na křivce pozorovat jisté význačné body, determinující prostorovou závislost sledovaného procesu. Tyto body bývají označovány jako parametry variogramu a jsou znázorněny na obr. 2.

Obr. 2 Parametry variogramu

•

bod A je označován jako rozsah („range“) a definuje vzdálenost, za kterou je již hodnota variogramu konstantní a body, které mezi sebou mají vzdálenost rovnou rozsahu a vyšší jsou vzájemně nekorelované,

24

•

bod B je nazýván „nugget effect“ (nepřekládá se) a je definován jako limita variogramu pro h → 0 . Příčinou nugget effectu může být buď chyba měření ve vstupním souboru, nebo skutečnost, že proces obsahuje prostorovou variabilitu na menším měřítku, než jaké jsme schopni postihnout pozorovanými daty,

•

bod C je označován jako práh („sill“) a odpovídá hodnotě variogramu pro vzdálenost rovnou rozsahu.

Skutečnost, že body se vzájemnou vzdáleností vyšší, než je hodnota rozsahu, se již ve svých hodnotách neovlivňují, svádí k tomu, aby pro odhad hodnoty v bodě s o byly použity jen ty body vstupního souboru, které jsou vzhledem k so „uvnitř“ rozsahu. Cressie (1993) proti tomu namítá, že tvar křivky variogramu je funkcí všech pozorovaných bodů, a proto by také měly být tyto body do odhadu zahrnuty. Kromě použití některého z definovaných modelů je samozřejmě možné vytvořit původní model variogramu. Při vytváření modelu 2γ (⋅) je však podmínkou, aby funkce exp[− a ⋅ γ (⋅)] byla pozitivně definitní pro jakékoliv kladné číslo a . Pro hlubší analýzu problematiky tvorby modelu variogramu viz Cressie (1993). V rámci diplomové práce nebyla tato možnost uvažována, ve výsledném softwaru je k dispozici výše zmiňovaných šest modelů. Následně po výběru konkrétního teoretického modelu je zapotřebí formulovat účelovou funkci. Tato funkce nějakým způsobem vyjadřuje rozdíl mezi experimentálními a modelovými hodnotami. Obvyklá definice účelové funkce (UF) je:

UF = ∑ ( y exp,i − y mod,i (λ )) n

2

(2.23)

i =1




y exp – experimentální hodnota




y mod (λ ) – modelová hodnota, funkcí vektoru parametrů modelu λ .

Řešení spočívá ve vyjádření parciálních derivací účelové funkce přes všechny hledané parametry a jejich položení rovno nule, čímž vzniká soustava tzv. normálních rovnic, jejímž řešením je vektor parametrů λ . Tento postup je znám pod názvem metoda nejmenších čtverců („ordinary least squares“). Možných definic

25

účelové funkce je více, k dispozici je např. metoda Čebyševova (minimalizace maximální chyby), metoda absolutních chyb atd. (Ralston 1978). Nalezení parametrů modelu některou z výše uvedených metod je možné, ale redukuje problematiku na čistě numerický proces, nereflektující její prostorově statistický kontext. Z tohoto důvodu zformuloval Cressie (1985) tvar účelové funkce, který obsahuje informaci o počtu párů bodů ze vstupního datového souboru, které se podílely na odhadu variogramu pro každou vzdálenost. Takto definovaná funkce automaticky „znevažuje“ ty odhady, které byly získány s pomocí malého počtu bodů, a dává větší váhu odhadům v okolí počátku. Tento postup byl pojmenován jako metoda vážených nejmenších čtverců („weighted least squares“) a jeho účelová funkce má tvar:  y exp,i  UF = ∑ N hi ⋅ − 1 i =1  y mod,i (λ )  n




2

(2.24)

N h i – počet párů bodů podílejících se na i-tém bodovém odhadu variogramu.

Do výsledného softwaru byla zapracována jednak metoda weighted least squares (wls) a jednak pro porovnání také obligátní ordinary least squares (ols). Pro nalezení minima účelové funkce je možné využít několika přístupů. Obecně jde o metody nederivační, metody derivační (gradientní) nebo metody Monte Carlo (náhodné). Pro účely diplomové práce byly uvažovány pouze gradientní metody. Tyto metody využívají vektoru parciálních derivací účelové funkce jako nástroje k postupnému přibližování se k optimu, tj. k její minimální hodnotě (Javůrek et Taufer 2006b). Po úvaze byla zvolena metoda největšího spádu. Popis její implementace do softwarové podoby je obsahem kapitoly 3.3.2. Pokud jsou provedeny odhady variogramu a jejich proložení křivkou, je tímto popsána prostorová závislost studovaného procesu a je možno přistoupit k sestavení rovnic krigingu. Jak ale bylo uvedeno v kapitole 2.2.1.2, odhadování variogramu z experimentálních dat předpokládá vnitřní stacionaritu procesu, kterou je nutno ověřit ještě před vlastním odhadem variogramu. Spolu s tímto ověřením je potřeba získat o datech další informace, které později vedou k výběru vhodné varianty krigingu. Analýze vstupních dat se věnuje následující kapitola.

26

2.2.2 Analýza vstupního datového souboru Teoretické principy vyjadřování prostorové závislosti a následných odhadů metodou kriging vycházejí z implicitního modelu prostorového procesu, který obsahuje předpoklad normality rozložení dat a vnitřní stacionaritu. Je tudíž předpokládáno, že pozorovaná data jsou výběrem z normálního rozdělení a jejich autokorelační struktura závisí pouze na relativních prostorových vztazích mezi body. Výrazné odchylky od tohoto modelu by pak měly být identifikovány a příslušná data vyloučena z dalších úvah a výpočtů (Cressie 1984). Pro potřeby kvalitního popisu prostorové závislosti a následného odhadu metodou kriging je nutné prověřit tyto tři otázky: 1. zda vstupní data obsahují hodnoty, které by porušovaly předpoklad stacionarity 2. zda je možné považovat prostorový proces za izotropní 3. zda je možné uvažovat konstantní střední hodnotu procesu. Za těmito účely vyvinul Tukey (1977) celou sadu nástrojů, které pojmenoval „Exploratory Data Analysis“ (EDA). Vylepšení některých metod pak navrhnul Cressie (1984). Většina těchto technik spočívá v různých způsobech sumarizace a grafického vynášení dat. Pro tyto účely Cressie (1984) požaduje, aby data ležela v pravidelné nebo téměř pravidelné mřížce. Pokud tomu tak není, je nutné data mřížkou proložit, tedy určit pro každý bod jemu příslušející řádek a sloupec. Principy několika z těchto postupů budou popsány v následujících odstavcích. 2.2.2.1 Detekce nestacionarity V této fázi je cílem najít takové body, jejichž hodnoty jsou nějakým způsobem extrémní vůči hodnotám bodů v jejich okolí. Tyto body musejí být pro účely dalšího zpracování dat vyloučeny, protože by nepřípustně ovlivnily výsledné odhady v okolních bodech (až do vzdálenosti dané rozsahem variogramu). Jednoduchým nástrojem, kterým je možno identifikovat přítomnost atypických hodnot je tzv. číslicový histogram („stem and leaf plot“). Histogram může vypadat např. takto:

27

3 |1 4 | 155 5 | 2333568 6 | 0011244578 7 | 022256689 8 | 002689 9 | 4558 10 | 11 | 12 | 9 Každá číslice v pravé („listové“ - leaf) části grafu odpovídá jednomu bodu, přičemž hodnota v bodě vznikne složením pravé části a levé („stonkové“ - stem) části. Údaj 9 | 4 tedy odpovídá hodnotě bodu 94. Z obrázku je patrné, že číslicový histogram podává hrubou představu o tvaru rozdělení souboru a navíc je možno pozorovat, že hodnota 129 vykazuje výraznou odlišnost od ostatních hodnot (Tukey 1977). Jiný způsob detekce odlehlých hodnot je založen na vlastnostech průměru a mediánu. Pokud jsou data rozmístěna v pravidelné mřížce, je možné vypočítat výběrový průměr a výběrový medián pro každý sloupec a řádek a tyto potom porovnat. Medián jakožto robustní míra polohy není ovlivněn přítomností extrémních hodnot, kdežto průměr z tohoto pohledu naopak citlivě reaguje. Velký rozdíl mezi hodnotami mediánu a průměru indikuje řádek či sloupec zasluhující pozornost. K posouzení velikosti tohoto rozdílu doporučuje Cressie (1984) využít statistiku:

(

D = n ⋅ (x − ~ x ) / σ ⋅ 0.5708




)

σ = (IQR )(2 ⋅ 0.6745) ,

kde IQR značí mezikvartilové rozpětí, přičemž hodnoty D kolem hodnoty tři a výše značí možný problém na příslušném úseku dat.

2.2.2.2 Izotropie Proces nazveme izotropním, pokud vykazuje v každém směru stejnou prostorovou závislost (Cressie 1993). V takovém případě je možné pro odhady variogramu využít interakce všech možných dvojic bodů vstupního datového souboru. Data však mohou

28

vykazovat v určitém směru anizotropii, tedy odlišnou autokorelaci od ostatních směrů. Potom je nutno zvolit pro odhady variogramu směr, který takové zatížení nevykazuje. Anizotropie obsažená v datech může být prokázána, ale její příčina je apriorně pokládána za neznámou, a je tedy modelována jako náhodná chyba (Kraus 2007). Výpočet výběrových mediánů pro každý řádek a sloupec a jejich následné vynesení do grafu je vhodným způsobem, jak detekovat případné směrové zatížení. Obrázek č. 3 ilustruje použití tohoto postupu.

Obr. 3 Detekce anizotropie

Z obrázku je patrný silný lineární trend při postupu ve směru osy x, ve směru osy y není směrové zatížení viditelné (Cressie 1984). 2.2.2.3 Trend Hodnotu Z (s ) v bodě s je možno podle (2.10) rozepsat jako součet trendové složky (střední hodnoty), která je vlastností procesu, a náhodné složky, která přísluší danému bodu v prostoru. Pro potřeby výběru vhodné varianty krigingu je nutné zjistit, zda je možno považovat trendovou složku za konstantní, nebo jestli je funkcí polohy bodu.

29

Technika zvaná mediánové zhlazování („median polish“) umožňuje tyto dvě složky oddělit a podrobit analýze každou zvlášť (Tukey 1977). Metoda vychází z faktu, že data umístěná v mřížce lze rozložit podle následujícího schématu: hodnota = příspěvek celého pole + příspěvek řádku + příspěvek sloupce + reziduum Mediánové zhlazování je iterační proces, kdy od každé hodnoty je postupně odečítán aktuální medián příslušného sloupce, potom medián příslušného řádku a tento postup je opakován až do chvíle, kdy nově vypočtené mediány začnou konvergovat k nule a hodnoty v poli se přestanou měnit. Hodnoty, které zbyly na místě původních bodů, představují hodnotu náhodné složky pro daný bod. Hodnotu trendové složky získáme jejím odečtením od původní hodnoty. Cressie (1984) doporučuje odhadnout variogram pro trendovou i náhodnou složku a tyto variogramy porovnat, čímž lze jasně prokázat nebo vyvrátit přítomnost celkového trendu v datovém souboru. Uvedené techniky, které byly vybrány z ilustračních důvodů, nepředstavují v žádném případě vyčerpávající výčet možných způsobů analýzy vstupních dat. Pro informaci o dalších způsobech (např. „variogram cloud“, „pocket plot“ atd.) viz Tukey (1977), Cressie (1984). V dalších kapitolách jsou stručně představeny základní varianty krigingu.

2.2.3 Výpočetní schéma krigingu Je dána sada n pozorování prostorového procesu Z : Ζ = {Z (s1 ), Z (s 2 ),..., Z (s n )}

(2.25)

Cílem je odhad (predikce) Z (s o ) , tj. neznámé hodnoty procesu Z v bodě s o . Funkční předpis, pomocí kterého je odhad prováděn, lze formálně zapsat: p(Ζ, s o ) = Z~(s o )

(2.26)




p (Ζ, s o ) – funkce (predikátor), jejíž nezávisle proměnné jsou Ζ – vektor pozorovaných hodnot a s o – známá lokace




Z~ (s o ) – odhad Z (s o ) .

30

V případě metody kriging se jedná o lineární predikci, neznámá hodnota Z (s o ) je odhadována jako lineární kombinace pozorovaných vstupních hodnot. Odhadovou funkci lze vyjádřit jako: n

p(Ζ, s o ) = ∑ α i ⋅ Z (s i )

(2.27)

i =1




α i – koeficienty lineární kombinace, neboli váhy, které jsou přiřazeny jednotlivým bodům vstupního souboru.

Kriging tedy spočívá ve váženém průměrování vstupních hodnot, kdy pro vyčíslení vah jsou zásadní hodnoty variogramu mezi bodem odhadu a vstupními body. Prvním požadavkem je nalezení optimální odhadové funkce, tedy takové, při jejímž použití, bude rozdíl skutečné hodnoty Z (s ) a vypočteného odhadu Z~ (s ) co o

o

možná nejmenší. Za tímto účelem je nutné definovat předpis, který tento rozdíl kvantifikuje. Jeho formální definici lze zapsat takto: L[Z (s o ), p (Ζ, s o )]




(2.28)

L – ztráta („loss“), která nastane, pokud je Z (s o ) odhadována pomocí

p (Ζ, s o ) .

Ztrátových funkcí je možno definovat mnoho, nejvýhodnější vlastnosti pro následné kalkulace má funkce definovaná jako kvadrát rozdílů:

L[Z (s o ), p(Ζ, s o )] = [Z (s o ) − p(Ζ, s o )]

2

(2.29)

Optimální predikátor je potom taková funkce, která minimalizuje

{

E{L} = E [Z (s o ) − p(Ζ, s o )]


2

}

(2.30)

E – střední hodnota, vzhledem k společnému rozdělení Z (s o ) a Ζ .

Druhým požadavkem při odhadu je jeho nevychýlenost, kterou lze formálně zapsat: E {p (Ζ, s o )} = E {Z (s o )}

(2.31)

Zpracováno podle Cressie (1993). Jednotlivé varianty krigingu se liší ve výchozích předpokladech o prostorovém procesu, s čímž souvisí rozdílné uplatnění požadavků optimality a nevychýlenosti

31

v algebraické podobě konečných rovnic. V následujícím textu je popsáno schéma dvou základních variant krigingu, které byly zahrnuty do vyvíjeného softwaru. 2.2.3.1 Ordinary kriging Tato varianta vychází z předpokladu konstantní střední hodnoty prostorového procesu:

Z (s ) = µ + δ (s )

(2.32)




µ – neznámá konstantní hodnota




δ (s ) – náhodná složka hodnoty Z (s ) , přičemž δ (⋅) je náhodný vnitřně stacionární proces s nulovou střední hodnotou a variogramem 2γ (⋅) .

Již bylo uvedeno, že kriging vyjadřuje odhad jako lineární kombinaci pozorovaných hodnot, viz (2.27). Odhad hodnoty Z (s o ) tedy spočívá v nalezení vektoru vah (α 1 , α 2 ,...α n ) ´. Nevychýlenost predikce ve smyslu (2.31) je potom vyjádřena podmínkou n

∑α i =1

i

(2.33)

=1

Hledání optimálního predikátoru tedy znamená minimalizaci ztrátové funkce (2.30) přes hledaný vektor vah při dodržení podmínky (2.33), což vede k formulaci minimalizačního kriteria ve tvaru: 2

n    n  E  Z (s o ) − ∑ α i ⋅ Z (s i ) − 2m ∑ α i − 1 i =1    i =1 




(2.34)

m – Lagrangeův multiplikátor, neboť vzhledem k (2.33) se jedná o tzv. vázaný extrém (např. Rektorys 1963).

S pomocí podmínky (2.33) může být řadou algebraických úprav odvozeno, že 2

2

[

2

]

n n n n   ( ) ( ) ( ) ( ) Z s − α Z s = − α α Z s − Z s / 2 + 2 α i [Z (s o ) − Z (si )] / 2 ∑∑ ∑ i  i j i j  o ∑ i i =1 i =1 j =1 i =1  

Kriterium (2.34) tudíž může být převedeno do tvaru:

∑∑ α iα j γ (si − s j ) + 2∑ α i γ (so − si ) − 2m ∑ α i − 1 n

n

i =1 j =1




n

i =1



n

 i =1

 

(2.35)

γ (⋅) – semivariogram 32

Po vyjádření parciálních derivací z (2.35) pro všechny váhy α i a pro m a jejich následném položení rovno nule je obdržena soustava lineárních rovnic, kterou je možno zapsat v maticovém tvaru: r r Μ⋅a = b




(2.36)

M – čtvercová matice o délce strany (n + 1), kde n je počet bodů vstupního souboru, matice je ve tvaru

M=




γ (s1 − s1 ) γ (s1 − s 2 ) γ (s 2 − s1 ) .

...

. γ (s n − s1 ) 1

. 1

...

γ (s1 − s n )

1 1

γ (s n − s n )

1 0

1

V souladu s definicemi teoretických modelů semivariogramu v kapitole 2.2.1.3 má matice na hlavní diagonále všechny prvky rovny nule a dále je




zřejmé, že se jedná o matici symetrickou. r a – sloupcový vektor hledaných vah

r a =

α1 α2 .

m – hodnota Lagrangeova multiplikátoru

αn m




r b – sloupcový vektor hodnot semivariogramu mezi bodem s o a vstupními

body

r b =

γ (s o − s1 ) γ (s o − s 2 ) . γ (s o − s n ) 1

Po vyřešení soustavy rovnic stačí dosadit vypočtené váhy α i do vztahu (2.27), čímž je obdržena výsledná hodnota odhadu: Z~ (s o ) = ∑ α i ⋅ Z (s i ) n

i =1




Z~ (s o ) - odhad hodnoty Z (s o ) .

33

Minimalizovanou střední chybu odhadu (2.30), jinak též nazývanou rozptyl odhadu je možno vyjádřit jako: n

σ 2 (s o ) = ∑ α i γ (s o − s i ) + m

(2.37)

i =1

Protože je předpokládáno normální rozdělení hodnot procesu Z (⋅) , je dále možné určit interval predikce, který s určitou pravděpodobností 1 − α pokrývá odhadovanou hodnotu. Jako α se obvykle volí 5 %. Šířka intervalu je pak dána vztahem:

I = [Z~ (s o ) − 1.96 ⋅ σ (s o ), Z~ (s o ) + 1.96 ⋅ σ (s o )]

(2.38)

Interval I pokrývá hodnotu Z (s o ) s pravděpodobností 95 %. pozn. Při dosazení

s o = s x do soustavy (2.36), kdy s x je jedním z bodů vstupního souboru, je

řešením soustavy


αx =1




α i = 0 , pro všechna i ≠ x

Odhadovaná hodnota v jednom ze vstupních bodů je tedy rovna hodnotě tohoto bodu, což je při převádění do softwarové podoby užitečným kontrolním nástrojem.

Zpracováno podle Cressie (1993).

2.2.3.2 Universal kriging Tato modifikace vychází z předpokladu, že trendová složka µ (s ) prostorového procesu podle (2.10) není konstantní, ale mění se s polohou bodu způsobem, který můžeme odhadnout. Hodnotu v bodě s je možno v takovém případě zapsat takto: p

Z (s ) = ∑ β i ⋅ f i (s ) + δ (s )

(2.39)

i =0




δ (s ) – náhodná složka hodnoty Z (s ) , přičemž δ (⋅) je náhodný vnitřně stacionární proces s nulovou střední hodnotou a variogramem 2γ (⋅)




f i (s ) – funkce, jejichž nezávisle proměnnými jsou souřadnice polohy bodu s




β i – konstanty.

Trendová složka procesu je tedy vyjádřena jako neznámá lineární kombinace (p+1) známých funkcí f 0 ,..., f p , přičemž tyto funkce mohou být konstantami.

V praxi bývá trend vyjadřován jako lineární kombinace polynomů zvoleného stupně,

34

obvykle stupně jedna (lineární trend) nebo dva (kvadratický trend), což jsou dva případy zahrnuté do vyvíjeného softwaru. V případě lineárního trendu je trendová složka vyjádřena ve tvaru:

µ (s ) = β 0 + β 1 x + β 2 y ,

(2.40)

v případě kvadratického trendu pak ve tvaru:

µ (s ) = β 0 + β1 x + β 2 y + β 3 xy + β 4 x 2 + β 5 y 2

(2.41)




x, y – souřadnice bodu s




funkce f 0 podle definice (2.39) je konstanta, f 0 = 1 .

Odhadovaná hodnota je opět stanovena jako lineární kombinace hodnot vstupního bodového souboru: Z~ (s o ) = ∑ α i ⋅ Z (s i ) , n

i =1

přičemž podmínku nevychýlenosti odhadu (2.31) je v tomto případě možno zapsat jako:

r

r

α ⋅F = f

(2.42)




α = (α 1 ,...α n ) – vektor vah




F – matice n x (p+1), přičemž n je počet bodů vstupního souboru, p je

r

stupeň polynomu zvoleného k popisu trendové složky

F=




f 0 (s1 )

...

f p (s1 )

.

.

.

f 0 (s n )

...

f p (s n )

r f = ( f 0 (s o ),..., f p (s o )) n

pozn. Pro stupeň polynomu p = 0 přechází podmínka (2.42) na vztah

∑α i =1

i

= 1 , tedy na podmínku

nevychýlenosti odhadu ordinary krigingu.

Odvození rovnic univerzálního krigingu je zcela analogické jako v případě ordinary krigingu, celý proces je popsán v Cressie (1993), str. 151 – 157. Opět jde o hledání optimálního vektoru vah, které spočívá v minimalizaci chyby odhadu (2.30), přičemž (p+1) funkcí lineární kombinace trendu způsobuje přítomnost (p+1)

35

Lagrangeových multiplikátorů. Proces minimalizace opět vede k soustavě lineárních rovnic, která může být v maticové formě zapsána jako: r r M ⋅a = b




(2.43)

M – čtvercová matice o délce strany (n + p + 1), kde n je počet vstupních bodů a p je stupeň polynomu, matice je ve tvaru:

γ (s1 − s1 )

... f 0 (s1 ) ... γ (s1 − s n ) . . . . . γ (s n − s1 ) ... γ (s n − s n ) f 0 (s n ) ... M = f 0 (s1 ) ... f 0 (s n ) 0 ... . . . . 0 f p (s1 ) ... f p (s n ) 0 .... Matice je symetrická s nulami na hlavní diagonále. r
a – hledaný vektor vah a Lagrangeových multiplikátorů

f p (s1 ) . f p (s n ) 0 . 0

α1 .

αn r a= m0 . mp


r b – sloupcový vektor ve tvaru:

γ (s o − s1 ) . γ (s o − s n ) r b= f 0 (s o ) . f p (s o ) Vyřešením soustavy a dosazením hodnot vah do (2.27) je vypočítána odhadovaná hodnota: Z~ (s o ) = ∑ α i ⋅ Z (s i ) n

i =1

Rozptyl odhadu je dán vztahem: n

p

i =1

j =0

σ 2 (so ) = ∑ α i γ (s o − s i ) + ∑ m j f j (s o )

(2.44)

Predikční interval:

36

I = [Z~ (s o ) − 1.96 ⋅ σ (s o ), Z~(s o ) + 1.96 ⋅ σ (s o )]

(2.45)

pokrývá hodnotu Z (s o ) s pravděpodobností 95 %. Výpočetní schéma univerzálního krigingu je zřejmé, problém je v tomto případě s vyjádřením variogramu. V situaci, kdy trendová složka procesu není konstantní, nelze k odhadu variogramu použít klasický (2.13) nebo robustní (2.14) způsob odhadu, protože oba tyto vztahy vychází z předpokladu, že hodnoty trendové složky se při výpočtu vynulují, což v případě procesu s nekonstantní střední hodnotou neplatí a celý odhad je zatížen chybou. Řešením tohoto problému je odečtení hodnoty trendu od původní hodnoty bodu a provedení odhadu variogramu pouze pro náhodnou složku. Hodnota trendu je ale neznámá, proto musí být odhadnuta vzhledem ke zvolenému polynomickému vyjádření. Zpracováno podle Cressie (1993). Postupů v tomto případě může být několik, do softwarové podoby byla implementována klasická lineární regrese podle principu:

β = (X T ⋅ X ) ⋅ X T ⋅ Ζ r

(2.46)

r




β – vektor odhadů koeficientů lineární kombinace polynomů




Ζ – vektor hodnot vstupního bodového pole




X – matice ve tvaru




X =

f 0 (s1 )

...

f p (s1 )

.

.

.

f 0 (s n )

...

f p (s n )

viz např. Rektorys (1963). V následující části diplomové práce bude popsáno převedení metod odhadů globální radiace a krigingu do softwarové podoby.

37

3. POPIS AUTORSKÉHO SOFTWARE Aplikace byla vytvořena v programovacím jazyce C++ ve vývojovém prostředí Dev-C++. Výsledným produktem je .exe soubor pracující pod operačním systémem Microsoft Windows. pozn. Vývojové prostředí Dev-C++ je produktem firmy Bloodshed software. Je určeno k vytváření aplikací pro operační systém Microsoft Windows, i když v současné době existují i jeho linuxové verze. Integrovaným kompilátorem v Dev-C++ je MinGW, pro více informací o tomto kompilátoru viz http://www.mingw.org. Celé prostředí Dev-C++ je volně dostupné na internetové adrese http://sourceforge.net/projects/dev-cpp/ a podle zkušeností autora této práce je vynikajícím nástrojem k vývoji aplikací v jazycích C a C++.

Protože jazyk C++ neobsahuje ve své standardní podobě žádné prvky grafického uživatelského rozhraní, je při vývoji aplikací nutné využít nějaké externí knihovny. Použitým nástrojem pro tvorbu grafického uživatelského rozhraní bylo Windows API (WinAPI) od společnosti Microsoft. Jazyk C++ podporuje tzv. objektově orientované programování (OOP), čehož bylo při vývoji aplikace využíváno. Není v možnostech této diplomové práce popsat filozofii OOP ani konkrétní implementaci jejích principů do jazyka C++. Na tomto místě je nutné alespoň stručně definovat základní pojem OOP, často používaný v dalším textu, kterým je třída. Třída („class“) je entita, zpodobňující nějaký objekt. Tímto objektem může být jakýkoliv objekt reálného světa i čistě abstraktní pojem. Definice třídy spočívá ve vymezení datových členů, kterými třída disponuje, a ve vymezení úkonů, které s datovými členy může provádět. Konkrétní instance dané třídy se nazývá objekt. Pro názornost může být takovou třídou např. třída matic, která byla do softwaru skutečně implementována. Datovými členy potom mohou být např. počet řádků a sloupců, hodnoty jednotlivých políček matice, hodnota normy matice atd. Mezi úkony, které je dovoleno s datovými členy třídy provádět, lze zařadit např. inverzi matice, násobení s jinou maticí, transpozice atd. Objektem třídy matic je pak nějaká matice A s konkrétními hodnotami svých datových členů. Při vývoji aplikace bylo využíváno dalších typických vlastností C++, především velmi účinného způsobu manipulace s pamětí počítače prostřednictvím odkazů a ukazatelů. Případný zájemce nalezne podrobný výklad této problematiky např. v Prata (2009).

38

V mnoha částech aplikace byl pro uložení a manipulaci s daty použit datový kontejner „vector“, který je součástí standardní knihovny C++. Tato struktura umožňuje dynamicky alokovat paměť a přidávat či ubírat prvky za běhu programu bez nutnosti opětovného definování a kopírování hodnot pole. Pro více informací o standardní knihovně (STD) a standardní knihovně šablon (STL) C++ viz Josuttis (2005). Vytvořený software byl pojmenován Astacus. Celá aplikace je z hlediska funkčnosti rozdělena do dvou víceméně nezávislých modulů. První modul byl vytvořen k výpočtu odhadů globální radiace, ve druhém je implementována interpolační

metoda

kriging.

Celková

struktura

aplikace

je

znázorněna

na následujícím diagramu:

Obr. 4 Celková struktura aplikace

Z obrázku je patrné, že software může být využit k samostatným odhadům radiace. Výsledky odhadů je buď možné exportovat ve formě textového souboru, nebo uložit do struktury bodového pole, která slouží jako úložiště dat pro následný vstup do modulu krigingu. Vstupní data pro kriging ovšem nemusí pocházet

39

z odhadů radiace, ale mohou být importována do bodového pole z vnějšího prostředí. Aplikaci je tedy možno využít pro interpolaci jakýchkoliv prostorových dat. V následujících kapitolách je stručně popsáno softwarové ztvárnění výše uvedených programových komponent. Tyto kapitoly se soustředí na slovní popis vnitřní struktury programu, návod k použití spolu s obrázky uživatelského rozhraní je obsažen v příloze č.1. Pro podrobnější informace je možno nahlédnout přímo do zdrojového kódu, který je spolu se zkompilovanou aplikací a diplomovou prací uložen na přiloženém CD. Pro tyto účely bylo snahou autora vkládat do kódu komentáře, stručně naznačující účel prováděných kroků, což se z větší části podařilo.

3.1 Odhady radiace V rámci tohoto modulu byla vytvořena struktura, která pro určitou lokalitu dovoluje provádět výpočty odhadů za delší období najednou a vytváří tak časové řady těchto odhadů. Podnětem k tomuto řešení je skutečnost, že časové řady meteorologických veličin jsou častým vstupem do hydrologických modelů. Pokud má aplikace najít uplatnění i v této oblasti, je uživatelsky nepřijatelné zadávat ručně hodnoty vstupních meteorologických prvků pro každý den, obzvláště pokud by cílem mělo být vytvoření např. desetileté časové řady. Jako vstupní formát byl zvolen formát textového souboru, do kterého je snadné exportovat data z MS Excel, což je nejrozšířenější nástroj pro úschovu a zpracování dat. Vstupy do výpočtu radiace zahrnují jednak potřebné informace o lokalitě (zeměpisné souřadnice, nadmořská výška atd., závisí na zvolené metodě), dále textový soubor, kde každý řádek obsahuje hodnoty meteorologických prvků z jednoho dne (viz Obr. 5). Informace o přesném tvaru textových souborů pro jednotlivé metody jsou k dispozici v hlavním menu softwaru v záložce Nápověda / Vstupní soubory.

40

Obr. 5 Schéma výpočtu radiace

Pro pohodlnou manipulaci s textovými soubory, které byly jako vstupní a výstupní formát použity na více místech aplikace, byly vytvořeny funkčně propojené třídy C++ s názvy „txt“ a „slovo“. Funkčně nadřazená třída „txt“ odpovídá celému textovému souboru, třída „slovo“ kompaktní skupině znaků, neobsahující mezeru, tabulátor ani znak nového řádku. Implementace třídy „txt“ je nejlépe patrná z výčtu jejích základních datových členů a funkcí. Datové členy: •

pole slov – obsah souboru převedený na jednotlivé objekty třídy „slovo“

•

počet slov, počet řádků

•

počet slov na každém řádku – nutné pro kontrolu vyžadovaného tvaru vstupních textových souborů.

Funkce: •

převedení obsahu souboru do jednotlivých slov

•

vracení konkrétního slova jako pole znaků

•

funkce podávající informace o tvaru souboru (počet slov, řádků, slov na jednotlivých řádcích atd.).

Pro kontrolu vstupů uživatele byly vytvořeny funkce, které zjišťují, zda zadaná hodnota (nebo položka textového souboru) je skutečně číslem. Všechny uživatelské vstupy jsou primárně uloženy jako znakové pole a následně jsou předány k vyhodnocení. Kontrolní funkce vrací číselnou hodnotu odpovídající zadanému

41

znakovému poli, pokud je toto pole v pořádku. Pokud není, vrátí funkce hodnotu -1 a informaci o chybném vstupu uloží do zvláštní proměnné. Struktura procesu odhadů radiace je shodná pro všechny zahrnuté metody. Uživatel zadá požadované údaje o lokalitě a nastaví cestu k textovému souboru obsahujícímu meteorologická data. Software nejdříve zkontroluje platnost zadaných hodnot a analyzuje textový soubor (počet údajů na řádcích, platnost číselných hodnot atd.). Pokud je vše v pořádku, provede výpočty odhadů podle příslušného algoritmu zvolené metody. Výsledkem je časová řada odhadů radiace, která je následně uložena opět ve formě textového souboru. Pokud jsou zjištěny závady ve vstupních datech, nebo pokud během procesu některá operace selže, je podáno příslušné chybové hlášení, což platí pro všechny softwarové struktury aplikace. V dialogovém okně každé metody je přítomna nabídka uložení výsledků do bodového pole, které slouží jako úložiště vstupních dat pro kriging. Pokud uživatel tuto nabídku zaškrtne, po skončení výpočtu a uložení textového souboru s výsledky se otevře dialog pro uložení dat do bodového pole. Do pole je možné uložit jednotlivé hodnoty odhadu, průměr nebo sumu celé časové řady. Zdrojový kód k algoritmům jednotlivých metod je k dispozici v souborech umístěných v adresáři

ZDROJOVY_KOD/RADIACE

na přiloženém CD. Kód ke třídám

„txt“, „slovo“ a ke kontrolním funkcím je v adresáři ZDROJOVY_KOD/TRIDY.

3.2 Bodové pole Pro ukládání prostorových dat byla vytvořena třída „bod“. Struktura této třídy je logická vzhledem k charakteru dat, která je v tomto případě nutné uchovávat. Její základní funkcionalitu přibližuje následující výčet některých datových členů a funkcí: Datové členy: •

identifikátor (název) – znakové pole, jednoznačně identifikující daný bod

•

hodnota v bodě

•

souřadnice x a y.

Funkce: •

nastavení a vracení hodnot datových členů

•

operátor porovnání názvů – pro účely vyhledávání v bodovém poli.

•

operátor přiřazení bodu. 42

Hodnoty lze do pole bodů vložit několika způsoby. Jednak jako výsledky odhadů radiace, tato možnost byla zmíněna výše. Dalším způsobem je přímé vkládání jednotlivých bodů, kdy jsou po otevření příslušného dialogového okna ručně zapsány hodnoty bodu. Třetí možností je načtení celého pole najednou z textového souboru. Je tedy možné si celé pole připravit v jiném prostředí (např. Excel), převést do textového souboru a ten importovat do aplikace. Jeden řádek textového souboru odpovídá hodnotám jednoho bodu, přičemž musí být dodrženo následující pořadí údajů: identifikátor – hodnota – souřadnice x – souřadnice y, které je při importu kontrolováno. Údaje o identifikátoru a hodnotě jsou povinné, souřadnice při importu dat v souboru být nemusí (v takovém případě zůstane datový člen příslušné souřadnice nedefinovaný, o čemž si objekt třídy „bod“ vede záznam ve speciální proměnné). Analogicky lze vytvořené pole exportovat z aplikace ven do textového souboru. Aktuální stav bodového pole lze za běhu programu kdykoliv zobrazit a lze ho editovat. Výčet editačních možností zahrnuje: •

přidání bodu

•

editace hodnot vybraného bodu

•

smazání bodu

•

import nového pole

•

export stávajícího pole

•

smazání celého obsahu pole najednou

•

vyhledávání bodů podle hodnoty identifikátoru

Zdrojový kód třídy „bod“ je uložen v adresáři ZDROJOVY_KOD/TRIDY. Po otevření modulu krigingu je možné načíst a zobrazit mapový podklad z externího souboru a na jeho základě umisťovat jednotlivé body do mapy pomocí myši (tedy přiřazovat jim souřadnice).

3.3 Kriging Popis programátorského ztvárnění krigingu je rozdělen do tří kapitol. První se týká odhadů variogramu a s tím spojené analýzy vstupních dat. Druhá kapitola je věnována prokládání odhadů variogramu a ve třetí kapitole jsou popsány procedury spojené se sestavením a vyřešením rovnic krigingu. 43

3.3.1 Odhady variogramu a analýza dat Do softwaru Astacus jsou implementovány oba způsoby odhadu variogramu popsané v kapitole 2.2.1.2, tedy klasický odhad podle vztahu (2.13) a robustní odhad podle (2.14). Byla uvažována varianta, že datový soubor může vykazovat anizotropii a že tedy může být nutné odhady variogramu provádět v určitém směru.Možný způsob detekce anizotropie je popsán v kapitole 2.2.2.2. Jeho podstatou je proložení datového souboru pravidelnou mřížkou a výpočet výběrových mediánů pro každý řádek a sloupec této mřížky. Z hodnot mediánů je možné usuzovat, zda data v příslušném směru vykazují či nevykazují směrové zatížení. Tento postup je v aplikaci implementován do dvou funkcí. První funkce přebírá jako vstupní hodnoty: •

bodové pole

•

zvolený počet řádků či sloupců výsledné mřížky.

Uživatel zadává jednak počet elementů a jednak směr (x nebo y). Na základě rozsahu hodnot souřadnic bodů ve zvoleném směru funkce vypočítá šířku elementu a proloží celé pole čtvercovou mřížkou. Následně vytvoří kopii původního bodového pole, ve které jsou souřadnice bodů vycentrovány do příslušných čtverců mřížky. Funkce zároveň podává informace o počtech bodů v jednotlivých řádcích a sloupcích, což spolu s grafickým výstupem dovoluje zhodnotit aktuálně zvolené proložení. Uživatel má možnost proceduru opakovaně provádět, dokud výsledek nesplňuje požadavky a body nejsou mřížkou proloženy optimálně. Druhá funkce přebírá již vycentrované bodové pole, počítá mediány pro každý řádek i sloupec a podává tyto informace uživateli, který se na jejich základě rozhoduje o směru odhadu variogramu. Uživatel má rovněž možnost nechat vypočítat pro každý řádek a sloupec absolutní hodnotu rozdílu mediánu a průměru a na základě těchto informací vyhledávat body vnášející do souboru nestacionaritu. Z pohledu směru odhadu jsou k dispozici tři možnost: •

odhad ve směru osy x

•

odhad ve směru osy y

•

všesměrný (izotopický) odhad.

44

Při odhadu v určitém směru algoritmus prochází vycentrované bodové pole, načítá hodnoty rozdílů bodů a počet interagujících dvojic pro všechny přítomné délky posunu mezi body a nakonec vypočítává odhady variogramu podle zvoleného způsobu (klasický či robustní odhad). Výsledkem jsou tři vektory hodnot: •

vektor odhadů variogramu

•

vektor délek posunu (vzdáleností)

•

vektor počtu interagujících dvojic.

Pro jednotlivé způsoby odhadu (klasický a robustní) jsou odhady pro oba směry prováděny jedinou funkcí, provádějící primárně odhady ve směru osy x. Při volbě směru osy y dojde k dočasné záměně souřadnic bodů. Při izotopickém odhadu variogramu jsou zahrnuty do výpočtu vzájemné interakce všech bodů datového souboru. Po výpočtu rozdílu hodnot pro každou dvojici bodů je v tomto případě nutné setřídit hodnoty do tříd četnosti podle vzdálenosti mezi body, přičemž každá třída četnosti je zastoupena jedinou hodnotou vzdálenosti. Odhady jsou poté počítány podle zvoleného způsobu. K určení počtu tříd bylo použito Sturgesova pravidla: k = 1 + 3.3 ⋅ log n




k – počet tříd




n – celkový počet bodů , viz Anděl (2007).

Při jakékoliv volbě parametrů odhadu (tedy směru a způsobu) jsou výsledkem procesu výše uvedené tři vektory, celý proces je možno schématicky znázornit:

Obr. 6 Schéma odhadů radiace

Tyto vektory jsou vstupními daty pro následující proces prokládání odhadů variogramu křivkou. Zdrojový kód k funkcím provádějícím analýzu dat je umístěn v adresáři

ZDROJOVY_KOD/ANALYZA

DAT,

k funkcím pro odhady variogramu

v adresáři ZDROJOVY_KOD/KRIGING.

45

3.3.2 Prokládání odhadů variogramu Jak bylo uvedeno v kapitole 2.2.1.3, prokládání variogramu spočívá v nalezení takových parametrů teoretického modelu, při kterých je minimalizován rozdíl mezi bodovými odhady a hodnotami vypočítanými podle modelu. Tento rozdíl je vyjádřen pomocí tzv. účelové funkce. Do softwarové podoby byly zahrnuty dva tvary účelové funkce, jednak funkce nejmenších čtverců (viz rovnice 2.23), a jednak funkce vážených nejmenších čtverců (2.25). Pro minimalizaci byla zvolena metoda největšího spádu. Při této metodě se hodnoty parametrů mění ve směru největšího poklesu účelové funkce, kterým je směrový vektor opačný k vektoru parciálních derivací podle jednotlivých parametrů. Algoritmus této metody lze popsat následující posloupností operací: 1. určení počátečního bodu P0, ze kterého algoritmus vychází 2. vypočtení gradientu G, tj. vektoru parciálních derivací účelové funkce v bodě P0 3. testování normy gradientu ||G|| vůči předem zvolenému konvergenčnímu kriteriu K a. pokud je ||G|| < K, výpočet je ukončen, b. pokud není ||G|| < K, následuje posun hodnot parametrů ve směru – G, bod P0 je nahrazen bodem P1 a následuje návrat ke kroku 2. Tato metoda přináší dobré výsledky, ale v okolí minima, které není příliš ostré, může algoritmus oscilovat (Javůrek et Taufer 2006b). Toto základní schéma má mnoho variant. Posun při kroku 3b je možné provádět buď tak dlouho, dokud dochází ke snižování účelové funkce (tuto vzdálenost je možné v jednoduchých případech spočítat i analyticky), nebo je zvolena maximální délka posunu, při jejímž dosažení je znovu počítán gradient. Test na ukončení je možné provádět nejen vzhledem k normě gradientu, ale také k velikosti snížení hodnoty účelové funkce mezi dvěma posuny apod. V případě prokládání hodnot variogramu je situace komplikovanější z toho důvodu, že parametry jednotlivých modelů nemohou nabývat jakýchkoliv hodnot, ale jsou definovány jako nezáporné (viz rovnice 2.17 – 2.22, kapitola 2.2.1.3). Parametr

46

a p mocninného modelu je ve svých hodnotách omezen oboustranně, neboť je

definován na intervalu a p ∈< 0,2) . Po konzultaci s doc. Mayerem z katedry numerické matematiky MFF UK byla provedena modifikace metody největšího spádu. Tato modifikace spočívá v tom, že pokud se algoritmus během výpočtu dostane do oblasti nepovolených hodnot jakéhokoliv parametru, je zjištěn průsečík s hranicí přípustné oblasti a z tohoto průsečíku jsou do prostoru vyslány směrové derivace. Tento postup (pro názornost pouze ve 2D prostoru) je znázorněn na obrázku číslo 7.

Obr. 7 Směrové derivace na hranici přípustné oblasti

Porovnáním směrových derivací je následně zhodnoceno, který ze zvolených směrů vede k největšímu snížení účelové funkce a další změny parametrů jsou prováděny tímto směrem. Specifická situace nastává při hledání parametrů sférického modelu. Po vyjádření účelové funkce (ve tvaru vážených nejmenších čtverců) je obdržen tento funkční předpis: 2

    2   [a s ] n  γ (h )  γ (h )  + UF = Nh ⋅  − 1 N ⋅ − 1  h   co + cs   3 h 1  h 3   h = [a s ]+1  h =1  c + c ⋅  ⋅ − ⋅     o s  2 as 2  as       

∑

∑

(3.1)

47

přičemž hledaným vektorem parametrů je λ = (co , c s , a s ) . Symbol [a s ] , přítomný v sumačních hranicích obou sum, vyjadřuje celočíselnou část hodnoty parametru a s , z čehož je patrné, že účelová funkce je vzhledem k parametru a s nespojitá a v bodech a s = l , l ∈ Z nemá derivaci. Tato situace byla vyřešena tak, že při minimalizaci je hledáno lokální minimum pro každý podprostor ohraničený hodnotami parametrů: •

co ≥ 0

•

cs ≥ 0

•

a s ∈ (l , l + 1) ,

přičemž l ∈ 0, n − 1 , kde n je počet bodových odhadů variogramu. Takto definované podprostory jsou znázorněny na obr. 8.

Obr. 8 Podprostory pro hledání lokálních minim

Porovnáním hodnot účelových funkcí pro všechna nalezená lokální minima je nakonec vybráno celkové minimum a proces je ukončen. Pozornost byla věnována otázce prvotního přiblížení, tedy určení bodu, ze kterého bude minimalizační proces startovat. Toto určení spočívá v „ovzorkování“ prostoru a nalezení bodu s nejmenší hodnotou účelové funkce ze všech testovaných. Tato operace je z výpočetního hlediska velmi náročná, obzvláště v trojrozměrném prostoru, proto bylo nutno determinovat rozsah možných hodnot parametrů. Všechny modely variogramu (kromě lineárního) mají tři parametry, v definicích 2.17 až 2.22 označeny jako co , c x , a x . Parametr co určuje vertikální posun ve směru osy variogramu, jeho hodnota může být interpretována jako nugget effect (viz kapitola 2.2.1.3). Parametry c x a a x se ve svém významu navzájem ovlivňují, ale obecně je možné říct, že c x určuje rozpětí hodnot variogramu, kdežto a x ovlivňuje 48

tvar výsledné křivky. Na mnoha velmi odlišných modelových datech byl zkoumán potřebný rozsah a hustota vzorkování pro jednotlivé modely variogramu. Výsledné nastavení funkcí je kompromisem mezi snahou o dosažení přijatelného výpočetního času a vysoké kvality výsledného proložení. Pozitivním faktem na analýzách modelových dat bylo zjištění, že účelové funkce vykazují plynulý průběh bez výrazných zlomů, neboť i při velmi rozdílném rozsahu a hustotě vzorkování prováděném na stejných datových souborech byly konečné výsledky minimalizace podobné, i když celý proces startoval ze zcela jiného bodu. Podobným

problémem

jako

počáteční

přiblížení

se

ukázalo

nastavení

konvergenčních kriterií. Všechny funkce provádějící minimalizaci obsahují tři konvergenční kriteria: •

kriterium normy gradientu – pokud je norma gradientu v aktuálním bodě menší než tato hodnota, výpočet zde skončí,

•

kriterium snížení účelové funkce – výpočet skončí, pokud je snížení účelové funkce mezi dvěma posuny menší než tato hodnota,

•

oscilační kriterium – výpočet skončí, pokud se po určitém počtu iterací proces neposune o tuto vzdálenost nebo vyšší.

Opět bylo provedeno velké množství testů, ve snaze najít optimální poměr mezi výpočetním časem a kvalitou výsledků. Pokud jsou kriteria nastavena na příliš velké hodnoty, dojde k ukončení algoritmu daleko od minima. Pokud jsou nastavena příliš malá, je výpočetní čas neúměrně dlouhý a algoritmus může začít oscilovat bez možnosti ukončení. Nastavené hodnoty pro jednotlivé modely variogramu by měly být zárukou přijatelného kompromisu, poskytujícího velmi dobré výsledky za krátkou výpočetní dobu. Do algoritmů funkcí byla přesto přidána další pojistka proti oscilaci. Ta spočívá v tom, že vždy po určitém (velkém) počtu iterací dojde ke zvýšení všech konvergenčních kritérií. Pokud tedy algoritmus přeci jen začne v nějaké oblasti oscilovat, po čase dojde vlivem navyšování kritérií k jeho ukončení.

49

Výsledný algoritmus funkcí provádějících minimalizaci je možné přiblížit pomocí následujícího schématu: 1. nastavení konvergenčních kriterií 2. vzorkování a nalezení počátečního bodu 3. výpočet gradientu G (vektoru parciálních derivací účelové funkce) 4. test normy G proti příslušnému kriteriu – norma G je menší než kriterium → konec – norma G je větší → bod 5. 5. test na zastavení případné oscilace – algoritmus osciluje → konec – neosciluje → bod 6. 6. spočtení hodnoty účelové funkce před posunem 7. posun ve směru - G o délku maximálního kroku (posunů je ve skutečnosti deset, každý o desetinu kroku, přičemž se sleduje hodnota účelové funkce a jako cíl posunu je vybrán bod s nejmenší hodnotou) 8. test snížení účelové funkce proti kriteriu – snížení je menší než kriterium → konec – snížení je větší → bod 9. 9. vyhodnocení nepřípustné oblasti – algoritmus s nachází v přípustné oblasti → návrat k bodu 3. – zjištěna nepřípustná oblast → bod 10. 10. spočtení průsečíku s hranicí oblasti a posun do tohoto bodu 11. výpočet směrových derivací a jejich zhodnocení – žádný směr nevede ke snížení účelové funkce → konec – alespoň jeden směr vede ke snížení účelové funkce → porovnání, derivace příslušná největšímu snížení je nastavena jako - G, následuje návrat k bodu 4. V případě sférického modelu je celý algoritmus opakován pro každý interval a s ∈ (l , l + 1) , viz výše.

Funkce přebírají jako svá vstupní data tři vektory, vytvořené během procesu odhadu variogramu (vektor variogramů, vektor vzdáleností, vektor počtu dvojic). Čtvrtou vstupní položkou je prázdný vektor, který funkce naplní hledanými hodnotami parametrů modelu:

50

Obr. 9 Schéma prokládání variogramu

3.3.3 Kriging Pro potřeby krigingu byla vytvořena třída „matice“, umožňující provádět operace maticového počtu, nutné k sestavení a výpočtu rovnic krigingu. Opět pro ilustraci následuje výčet základních datových členů a funkcí třídy: Datové členy:


počet řádků, počet sloupců




pole – dvourozměrné pole prvků matice.

Funkce:


inverze matice (s výběrem hlavního prvku)




transpozice matice




operátory sčítání a násobení matic




řešení soustavy lineárních rovnic




výpis matice do textového souboru, načtení matice z textového souboru.

3.3.3.1 Ordinary kriging Proces sestavení a řešení rovnic ordinary krigingu je koncipován jako jedna funkce, stejně tak je tomu v případě universal krigingu. Výsledkem interpolace může být buď odhad v jednom bodě (tedy číslo), nebo rastr prostorového rozložení odhadů. Pro obě tyto varianty byla použita stejná funkce, jejímž výsledkem je matice hodnot, v případě bodového odhadu má tato matice velikost 1x1. Vstupními hodnotami pro finální proceduru jsou: •

bodové pole

•

číselná hodnota identifikující zvolený model variogramu

•

vektor parametrů tohoto modelu, získaný v procesu prokládání odhadů variogramu

51

•

souřadnice levého horního a pravého spodního rohu pokrývané oblasti (v případě bodového odhadu jsou souřadnice totožné)

•

prázdná matice pro hodnoty odhadů

•

prázdná matice pro hodnoty rozptylů odhadu. Funkce naplní obě prázdné matice výslednými hodnotami.

Průběh funkce spočívá v sestavení a vyřešení soustavy lineárních rovnic ve tvaru r r M ⋅ a = b , viz rovnice (2.36), kapitola 2.2.3.1. Prvním krokem je sestavení matice M , tedy matice hodnot semivariogramu mezi r

body vstupního bodového pole. Následuje sestavení vektoru b , jehož prvky jsou hodnoty semivariogramu mezi bodem odhadu a body vstupního pole. Poté je vyřešena soustava (2.36). Z několika možných způsobů řešení byla vybrána Gaussova eliminační metoda s výběrem hlavního prvku, neboť tento způsob maximálně omezuje vliv zaokrouhlovacích chyb, ke kterým nutně dochází při řešení soustavy rovnic na počítači (Ralston 1978). Výsledkem je vektor ar , jehož prvky jsou hodnoty vah příslušející jednotlivým bodů vstupního pole a hodnota Lagrangeova multiplikátoru. Prvky vektoru ar jsou dosazeny do rovnic (2.27) a (2.37), ze kterých jsou vypočteny hodnoty odhadu a rozptylu odhadu. Takto popsaný postup ilustruje řešení při výpočtu jednoho bodového odhadu. V situaci, kdy je odhady pokrývána souvislá oblast, je uvedený postup modifikován. Při volbě oblasti o m řádcích a n sloupcích musí funkce provést m ⋅ n odhadů podle výše popsaného schématu. Čísla m a n vyjadřují velikost vytvářeného rastru, celkový počet odhadů může být roven několika milionům. Při všech těchto odhadech je r

matice M totožná (protože vstupní bodové pole je stejné) a mění se jen vektor b , který je závislý na poloze aktuálně odhadovaného bodu. Provedená modifikace spočívá v řešení všech soustav pro všechny odhadované body najednou viz schéma na obrázku 10.

Obr. 10 Řešení soustavy lineárních rovnic

52

Ze schématu je patrné, že matice M je sestavena pouze jednou na počátku r procesu. V dalším průběhu jsou vyčíslovány vektory b pro jednotlivé body oblasti, ze kterých je sestavována matice B. Gaussova eliminace je pak provedena pro všechny sloupce matice B najednou, výsledkem (po zpětné substituci) je matice r řešení A, ve které každý sloupec odpovídá vektoru řešení a pro příslušný bod oblasti. r Po vyřešení soustavy jsou dosazením příslušných vektorů a do rovnic (2.27) a

(2.37) pro každý bod oblasti naplněny obě vstupní matice, tj. matice odhadů a rozptylů odhadů. 3.3.3.2 Universal kriging Struktura funkce universal krigingu je zcela shodná jako v případě ordinary krigingu, pouze místo soustavy rovnic (2.36) je řešena soustava (2.43), čemuž samozřejmě odpovídá příslušný tvar sestavovaných matic. Rozdíl mezi předešlou variantou krigingu spočívá v odhadech variogramu, prováděných před vstupem do finální procedury. Jak bylo uvedeno v kapitole 2.2.3.2, v situaci s nekonstantní trendovou složkou je nutné před odhadem variogramu odečíst hodnoty trendu od hodnot vstupních bodů a provádět odhady pro zbylá rezidua. Tento problém byl vyřešen vytvořením dvou funkcí, které jsou přidány do běhu programu před odhady variogramu. První funkce na základě uživatelem zvoleného polynomického vyjádření trendu provede odhady koeficientů lineární kombinace polynomů podle rovnice (2.46). Druhá funkce potom vypočítá rezidua všech bodů vstupního pole a vloží je jako hodnoty bodů do kopie bodového pole. Tato kopie následně vstupuje do odhadů variogramu. Všechny další procesy jsou již stejné jako v případě ordinary krigingu (přičemž do vlastní funkce universal krigingu již samozřejmě vstupuje původní bodové pole).

53

Celkové schéma procedur krigingu je patrné z obrázku číslo 11.

Obr. 11 Schéma funkcí krigingu

Ke všem výše popsaným komponentám softwaru bylo vytvořeno grafické uživatelské rozhraní, které umožňuje intuitivní ovládání aplikace pomocí tlačítek, dialogových oken, rozbalovacích seznamů a podobných standardních prvků. Dále byla ke všem prvkům uživatelského rozhraní vytvořena mutace v anglickém jazyce, uživatel má možnost jednoduchým způsobem přepnout celou aplikaci do angličtiny. Ukázky uživatelského prostředí, které zároveň mohou sloužit jako návod k použití softwaru jsou obsahem přílohy č. 1.

54

4. VÝSLEDKY – INTERPOLACE RADIAČNÍCH ODHADŮ V této části diplomové práce jsou představeny výsledky praktického použití softwaru Astacus. Cílem bylo analyzovat, zda je možné interpolační metodou kriging získat kvalitní odhady radiace v situaci, kdy samotnými vstupními daty pro interpolaci jsou odhady této veličiny. Pro tyto účely bylo z meteorologických stanic Českého hydrometeorologického ústavu (ČHMÚ) vybráno 91 stanic, které víceméně rovnoměrně pokrývají území České republiky. Přehled vybraných stanic je obsažen v příloze 2. Pro období jednoho roku (1991) byly pro každý den provedeny z pozorovaných dat na těchto stanicích odhady radiace. Tyto odhady byly následně interpolovány s různým nastavením parametrů krigingu. Takto získané hodnoty byly porovnány s měřenými hodnotami radiace. Z důvodů omezených finanční prostředků na nákup dat byly nejprve porovnány metody odhadu radiace mezi sebou a byla vybrána nejlepší z nich. Následně byla zakoupena data (z výše zmíněných 91 stanic) pouze těch meteorologických veličin, které vstupují do vybrané metody.

4.1 Porovnání odhadů radiace Pro porovnání byla použita měřená data ze stanice České Budějovice z let 1961 1998, kterými disponuje pracoviště autora této práce (Ústav pro hydrodynamiku AVČR). K dispozici byla v tomto případě data potřebná pro všechny metody odhadu i měřené hodnoty radiace. Pro porovnání výsledků bylo použito tří statistik. První byla suma rozdílů (SR): SR = ∑ x měě − xodh •

x měě – měřená hodnota [MJ . m-2 . den-1]

•

xodh – odhadovaná hodnota [MJ . m-2 . den-1]

(4.1)

Druhou byl průměrný denní rozdíl (PDR): PDR = •

∑x

měě

− xodh

n

(4.2)

n – počet záznamů

Poslední byla průměrná denní odchylka (PDO):

55

PDO =

x měě − xodh 1 ∑ n x měě

(4.3)

Pro představu dodejme, že pokud je radiace vyjádřena v [MJ . m-2 . den-1], její hodnoty na území České republiky se pohybují zhruba v rozmezí od minima 1 - 5 v zimním období do maxima 30 - 35 v letním období. Výsledné hodnoty uvedených statistik pro posuzované metody jsou uvedeny v tabulce č. 1. Tabulka č. 1. Porovnání metod odhadu radiace na základě porovnání s měřenými daty Winslow Klabzuba Hargreaves Angstrom Supit -2 -1 3 3 3 3 3 SR [MJ.m .den ] 31.3 · 10 82.1 · 10 35.9 · 10 30.1 · 10 33.6 · 10 -2 -1 PDR [MJ.m .den ] 2.25 5.91 2.58 2.22 2.42 PDO [-] 0.33 0.65 0.40 0.25 0.36 Pozn. SR – suma rozdílů, PDR – průměrný denní rozdíl, PDO – průměrná denní odchylka.

Z výsledků vyplývá, že jednoznačně nejlepší odhady poskytuje metoda Angstromova. Metoda Klabzubova dává dobré výsledky v letním období, což je v souladu s jejím původním účelem (použitím v růstových modelech plodin), v chladných částech roku se kvalita odhadů prudce zhoršuje, což značně ovlivnilo výsledné statistiky. Zbylé tři metody se pak výsledkově pohybují mezi Angstromem a Klabzubou. Představu o rozdílné přesnosti metod Angstroma a Winslowa, které byla podle statistik hodnoceny jako dvě nejlepší, je možno získat vizuálním porovnáním dvou grafů (Obr. 12 a Obr. 13), ve kterých jsou proti sobě vyneseny měřené hodnoty a odhady za rok 1991 ve stanici České Budějovice.

56

40 35

-1

odhady [MJ.m .den ]

30

-2

25 20 15 10 5 0 0

5

10

15

20

25 -2

30

35

-1

měřené hodnoty [MJ.m .den ]

Obr. 12 Výsledky odhadů radiace Angstromovy metody versus měřené hodnoty ve stanici České Budějovice 35

25

-2

-1

odhady [MJ.m .den ]

30

20

15

10

5

0 0

5

10

15

20

25 -2

30

35

-1

měřené hodnoty [MJ.m .den ]

Obr. 13 Výsledky odhadů radiace Winslowovy metody versus měřené hodnoty ve stanici České Budějovice

57

Celkově lze poznamenat, že všechny porovnávané metody jsou schopny poskytovat kvalitní odhady radiace, každá má své opodstatnění při určité dostupnosti vstupních údajů. Pro interpolace v této práci byla vybrána metoda Angstromova.

4.2 Porovnání variant krigingu Pro vytvoření sítě vstupních bodů byly použity záznamy doby trvání slunečního svitu pro výše zmíněných 91 stanic za rok 1991. Z těchto záznamů byly spočítány metodou Angstrom - Prescott odhady radiace pro každý den. Vzniklo tedy 365 vstupních souborů, každý obsahující 91 bodů v prostoru. Všechny tyto soubory byly analyzovány za účelem detekce případné nestacionarity a anizotropie. Na výsledných souborech byly následně testovány různé varianty nastavení krigingu a interpolované hodnoty byly opět porovnány s měřenými hodnotami radiace ve stanici České Budějovice za rok 1991. Prvním učiněným závěrem je fakt, že robustní odhad variogramu dává podstatně lepší výsledky než odhad klasický. Rovněž tak prokládání odhadů variogramu metodou vážených nejmenších čtverců znamená mírné zlepšení odhadů oproti standardní metodě nejmenších čtverců, i když rozdíl není tak výrazný jako v případě odhadů variogramu. Všechny následující výsledky byly proto dosaženy s nastavením robustního odhadu a vážených nejmenších čtverců. Pozornost byla dále soustředěna na výběr varianty krigingu (ordinary – universal lineární trend – universal kvadratický trend) a výběru modelu variogramu (výše uvedených šest modelů). Pro každou kombinaci byly provedeny interpolace každého ze 365 vstupních souborů a výsledné hodnoty byly porovnány s měřenými daty pomocí tří statistik, použitých pro porovnání metod odhadu radiace (formule 4.1 - 4.3). Přehled výsledků porovnání je uveden v tabulce č. 2.

58

Tabulka č. 2 Porovnání bodových odhadů radiace s měřenými daty při různém nastavení parametrů krigingu ordinary kriging variogram SR PDR PDO lineární 1113.29 3.05 0.39 exponenciální 1120.69 3.07 0.39 sférický 1076.01 2.95 0.39 kvadratický 1064.97 2.92 0.39 wave 1117.12 3.06 0.40 mocninný 1129.45 3.09 0.40 universal kriging - lineární trend lineární 1180.20 3.23 0.41 exponenciální 1160.98 3.18 0.40 sférický 1230.70 3.37 0.44 kvadratický 1227.03 3.36 0.43 wave 1221.15 3.35 0.43 mocninný 1196.21 3.28 0.41 universal. kriging - kvadratický trend lineární 1159.83 3.18 0.41 exponenciální 1163.94 3.19 0.41 sférický 1230.70 3.37 0.44 kvadratický 1229.08 3.37 0.43 wave 1236.44 3.39 0.40 mocninný 1208.17 3.31 0.43 Pozn. SR – suma rozdílů, PDR – průměrný denní rozdíl, PDO – průměrná denní odchylka.

Při pohledu na statistiky lze konstatovat, že všechny porovnávané varianty nastavení krigingu dávají velmi dobré a také velmi podobné výsledky. Obecně mírně lepší hodnoty dává ordinary kriging, obzvláště kvadratický variogram, ovšem rozdíly mezi jednotlivými nastaveními byly většinou pod úrovní chyby odhadu. Společným rysem testovaných variant bylo také to, že všechny dávají mírně nadhodnocené odhady proti skutečnosti, zvláště v období letních maximálních hodnot. Tento fakt ilustruje obrázek č. 14, který zobrazuje v jednom grafu měřené hodnoty radiace a odhady ordinary krigingem s kvadratickým variogramem.

59

měřené hodnoty

kvadratický

40

radiace [MJ.m -2.den -1]

35 30 25 20 15 10 5 0 0

50

100

150

200

250

300

350

den Obr. 14 Porovnání měřených hodnot a odhadů radiace ordinary krigingem s kvadratickým variogramem

Analogické grafy pro ostatní varianty krigingu vypadají velmi podobně. Při interpretaci výsledků je však potřeba zdržet se všeobecných závěrů, neboť byly získány na základě porovnání s daty z jediné meteorologické stanice, s čímž může souviset i zmíněné nadhodnocování, které může reflektovat specifickou polohu lokality. Výsledky přesto opravňují k závěru, že interpolace radiačních odhadů je schopna poskytnout věrohodnou představu o radiaci na zkoumané lokalitě.

60

5. DISKUZE Předložená diplomová práce měla dva hlavní cíle. Prvním cílem bylo naprogramovat software, provádějící jednak odhady radiace z meteorologických prvků a jednak interpolaci metodou kriging. Výsledný software se skládá ze dvou modulů, první v sobě zahrnuje pět metod odhadu radiace a umožňuje vytvářet časové řady těchto odhadů. Z tohoto pohledu tedy může nalézt uplatnění při tvorbě vstupních dat pro hydrologické modely. Druhý modul provádí interpolační metodu kriging. Obsahuje obě základní varianty krigingu, široký výběr modelů variogramu a umožňuje různé způsoby odhadu variogramu a prokládání těchto odhadů. V porovnání s nejrozšířenějším komerčním softwarem na poli geografických informačních systémů, kterým je ArcGIS – ArcView od společnosti ESRI, obsahuje software Astacus širší nabídku modelů variogramu a dovoluje výběr modelu variogramu i v případě universal krigingu. Nejvýraznějším rozdílem je ale skutečnost, že vyvinutý software obsahuje nástroje umožňující analyzovat vstupní datový soubor. Tato analýza je zaměřena na odhalení případné nestacionarity a anizotropie a software také dovoluje provádět směrové odhady variogramu i pro data neležící v pravidelné mřížce. Na druhé straně je ale nutné zmínit podstatnou nevýhodu softwaru, kterou je absence kartografických souřadných systémů. Astacus pohlíží na prostor jako na matici, ve které jsou prvky tvořeny jednotlivými pixely vytvářeného rastru. V důsledku zakřivení Země je ale zobrazení zemského povrchu na rovinnou plochu spojeno s chybou. Pro oblast velikosti České republiky je ale při přesném umisťování bodů chyba malá a rozhodně není limitujícím faktorem kvality výsledných odhadů. Z důvodu absence souřadných systémů je také žádoucí, aby použité mapové podklady importované do softwaru co možná nejméně zkreslovaly prostorové vztahy zájmové oblasti, které jsou pro kriging určující. Aby k takovému zkreslení nedošlo vinou samotného softwaru, byla funkce zobrazující bitmapy koncipována tak, že pokud je načítaná bitmapa v jakémkoliv směru větší než zobrazovací okno, přizpůsobí se tomuto oknu vždy při zachování poměru stran, relativní prostorové poměry mapového podkladu tak zůstanou zachovány. Druhým cílem diplomové práce bylo pomocí vyvinutého softwaru posoudit, zda je možné získat kvalitní data globální radiace interpolací radiačních odhadů. Nejprve

61

bylo provedeno porovnání samotných metod odhadu radiace, které jednoznačně upřednostnilo metodu Angstromovu, což je v souladu s Trnka et al. (2005). Závěrem následných interpolací a porovnání výsledků s naměřenými hodnotami je zjištění, že kriging může poskytovat kvalitní odhady radiace i ve chvíli, kdy nejsou k dispozici žádná radiační data a vstupními daty pro interpolaci jsou pouze odhady radiace z jiných meteorologických veličin. Je však třeba mít na zřeteli, že výsledky byly získány porovnáním interpolovaných hodnot měřenými daty pouze v jedné stanici, mohou tedy být ovlivněny její specifickou polohou. Jak již bylo řečeno v úvodu, kriging zahrnuje širokou škálu interpolačních technik, přičemž obě implementované varianty tvoří základ, ze kterého ostatní modifikace vycházejí. V případě interpolací radiačních odhadů přichází do úvahy obzvláště tzv. cokriging, který vychází z předpokladu vzájemné závislosti několika zkoumaných veličin. Odhady radiace by zřejmě bylo možné zpřesnit např. pomocí údajů o nadmořské výšce, neboť podle Martínez - Lozano et al. (1994) je hodnota radiace s hodnotou nadmořské výšky korelována.

62

6. ZÁVĚR Předložená diplomová práce splnila všechny vytyčené cíle. Byla provedena rešerše problematiky odhadů radiace a rešerše týkající se interpolační metody kriging. Na základě získaných poznatků byl naprogramován software, umožňující vytvářet časové řady odhadů radiace pomocí pěti vybraných metod a dále provádějící interpolace metodou kriging. Software zahrnuje obě kmenové varianty krigingu a také nástroje pro analýzu dat. Praktickým použitím vytvořeného softwaru byla zkoumána možnost získat údaje globální radiace pomocí interpolací odhadů této veličiny. Výsledky ukazují, že tato možnost je reálná a kriging může poskytovat kvalitní data globální radiace v situaci, kdy vstupními hodnotami jsou pouze radiační odhady. Autor práce bude nadále software rozvíjet, neboť kriging zahrnuje řadu dalších technik, které vychází z již naprogramovaných variant. V plánu

je také

implementovat do softwaru i jiné interpolační metody (IDW, spline, atp.) a především zavedení kartografických souřadných systémů.

63

Seznam literatury Allen, R.G., Pereira, L.S., Raes, D., Smith, M., (1998): Crop evapotranspiration: guidelines for computing crop requirements. Irrigation and Drainage Paper No. 56, FAO, Rome, Italy, 300 str. Anděl, J. (2007): Satistické metody. Matfyzpress: 299 str. Angstrom, A. (1924): Solar and terrestrial radiation. Q.J.R.Meteorol. Soc., 50: 121125. Angstrom, A. (1956): On the computation of global radiation from records of sunshine. Ark. Geofys., 2: 471 - 479. Bristow, K.L., Campbell, G.S. (1984): On the relationship between incoming solar radiation and daily maximum and minimum temperature. Agric. For. Meteorol., 31: 159-166. Cressie, N.A.C. (1984): Towards resistant geostatistics. In: G. Verley et al. (eds.): Geostatistics for natural resources characterization, Dordrecht, Reidel: 21:44. Cressie, N.A.C. (1985): Fitting variogram models by weighted lest squares. Mathematical Geology, 17,5: 563 – 586. Cressie, N.A.C. (1990):The origins of kriging. Mathematical geology, 22: 239 - 252. Cressie, N.A.C. (1993) : Statistics for spatial data. John Wiley and Sons, New York: 900 str. Cressie, N.A.C., Hawkins, D.M. (1980): Robust estimation of the variogram. Journal Inter. Assoc. Math. Geol., 12: 115 - 125. Gandin, L.F. (1963): Objective analysis of meteorological fields. GIMIZ, Leningrad. Glover, J., McCulloch, J.S.F. (1985): The empirical relation between solar radiation and hours of sunshine. Q.J.R. Meteorol. Soc., 84: 172 - 175. Hargreaves, G.L., Hargreaves, G.H., Riley, J. P. (1985): Irrigation water requirement for the Senegal River Basin. J. Irrig. Drain. Eng. ASCE, 111: 265–275. Hubbard, K.G. (1994): Spatial variability of daily weather variables in the high plains of the USA. Agric. For. Meteorol., 68: 29 - 41 Chalupa, R. (2003): 1001 tipů a triků pro Visual C++. Computer press, Brno: 434 str. Javůrek, M., Taufer, I. (2006a): Nebojme se nelineární regrese (1). Chemagazín 16,2: 27 - 29. Javůrek, M., Taufer, I. (2006b): Nebojme se nelineární regrese (2). Chemagazín 16,4: 13 - 14. Josuttis, M.N. (2005): C++ Standardní knihovna a STL. CP Books, Brno: 741 str. Klabzuba, J.(2001): Aplikovaná meteorologie a klimatologie, IV. Díl - Záření Slunce, Země a atmosféry. Česká zemědělská univerzita v Praze: 43 str. Klabzuba, J., Bureš, R., Kožnarová, V. (1999) In: Proceedings of the ‘‘Bioklimatologické pracovné dni 1999 Zvolen’’, Model výpočtu denních sum globálního záření pro použití v růstových modelech, 121–122.

64

Klabzuba, J., Kožnarová, V. (1991): Zářivá energie jako faktor mikroklimatu porostu. Editpress, Vysoká škola zemědělská, Praha: 118 str. Kraus, J. (2007): Geostatistika jako prostorové modelování statistických jevů. Statistika, 44, 6: 490-502. Liu, D.L., Scott, B.J. (2001): Estimation of solar radiation in Australia from rainfall and temperature observations. Agric. For.Meteorol., 106: 41–59. Martínez-Lozano, J.A., Tena, F., Onrubia, J.E., De la Rubia, J. (1994): The historical evolution of the Angström formula and its modifications: review and bibliography. Agric. For. Meteorol., 33: 109-128. Matheron, G. (1963): Principles of geostatistics. Economic Geology, 58: 1246--1266 Prata, S. (2009): Mistrovství v C++. Computer press, Brno: 1119 str. Prescott, J.A. (1940): Evaporation from a water surface in relation to solar radiation. Trans. R. Soc. Sci. Aust., 64: 114 - 125. Ralston, A. (1978) Základy numerické matematiky. Academia, Praha: 635 str. Rektorys, K. (1963): Přehled užité matematiky. Státní nakladatelství technické literatury: 1140 str. Sellers, W.D. (1965): Physical Climatology. The University of Chicago Press, Chicago: 272 str. Sobíšek, B. (1993): Meteorologický slovník výkladový a terminologický. MŽP ČR: 594 str. Soltani, A., Meinke, H., de Voil, P. (2003): Assesing linear interpolation to generate daily radiation and temperature data for use in crop simulations. Eur. J. Agron., 21: 133–158. Soukupová, J. (2009): Atmosférické procesy (základy meteorologie a klimatologie). Česká zemědělské univerzita v Praze: 201 str. Supit, I. (1997): Predicting national wheat production volumes using growth simulation results plus a trend function. Agric. For. Meteorol., 88: 199 - 214. Supit, I., van Kappel, R.R. (1998): A simple method to estimate global radiation. Solar Energy 63: 147–160. Thornton, P.E., Running, S.W., White, M.A. (1997): Generating surfaces of daily meteorological variables over large regions of complex terrain. J. Hydrol. 190: 214–251. Thorton, P.E., Running, S.W. (1999): An improved algorithm for estimating incident daily solar radiation from measurements of temperature, humidity and precipitation. Agric. For. Meteorol., 93: 211–228. Trnka, M., Žalud, Z., Eitzinger, J., Dubrovský, M. (2005): Global solar radiation in Central European lowlands estimated by various empirical formulae. Agricultural and Forest Meteorology, 131: 54–76. Tukey, J. (1977): Exploratory data analysis. Addison - Wesley: 688 str. Vaníček, K. (1994): Popis pole globálního záření na území České republiky v období 1984 - 1993. Český hydrometeorologický ústav, Praha: 85 str.

65

Winslow, J.C., Hunt, E.R., Piper, S.C. (2001): A globally applicable model of daily solar irradiance estimated from air temperature and precipitation data. Ecol. Modelling, 143: 227-243.

66

Seznam obrázků Obr. 1 Bodové odhady variogramu Obr. 2 Parametry variogramu Obr. 3 Detekce anizotropie Obr. 4 Celková struktura aplikace Obr. 5 Schéma výpočtu radiace Obr. 6 Schéma odhadů radiace Obr. 7 Směrové derivace na hranici přípustné oblasti Obr. 8 Podprostory pro hledání lokálních minim Obr. 9 Schéma prokládání variogramu Obr. 10 Řešení soustavy lineárních rovnic Obr. 11 Schéma funkcí krigingu Obr. 12 Výsledky odhadů radiace Angstromovy metody versus měřené hodnoty ve stanici České Budějovice Obr. 13 Výsledky odhadů radiace Winslowovy metody versus měřené hodnoty ve stanici České Budějovice Obr. 14 Porovnání měřených hodnot a odhadů radiace ordinary krigingem s kvadratickým variogramem

67