Univerza v Ljubljani Fakulteta za elektrotehniko Boštjan Grašič Napovedovanje povišanih koncentracij ozona z uporabo umetnih nevronskih mrež, Gaussovih procesov in mehke logike Magistrsko delo Mentor: prof. dr. Marko Munih Somentor: prof. dr. Drago Matko V Ljubljani, marec 2005 Zahvala Iskreno se zahvaljujem prof. dr. Marku Munihu za vso pomoč v času študija in pri izdelavi magistrske naloge ter prof. dr. Dragu Matku za vso pomoč pri izdelavi seminarske in magistrske naloge. Še posebej se zahvaljujem tudi sodelavcema dr. Primožu Mlakarju in dr. Mariji Božnar. Skupaj smo predebatirali marsikateri problem s področja okoljskega modeliranja, kar je bilo podlaga za raziskave zajete v tej magistrski nalogi. Prav tako se zahvaljujem tudi vsem sodelavcem podjetja AMES d.o.o, ki so mi omogočili delo na področju okoljskega monitoringa in mi nudili ustrezne materialne pogoje za delo. i ii Povzetek Namen magistrskega dela je izgradnja modela za vnaprejšnje opozarjanje ljudi na pričakovane presežene opozorilne in alarmne vrednosti koncentracije ozona. Splošno problematiko na tem področju prikazuje uvodno poglavje. V nasprotju z višjimi plastmi atmosfere, kjer je ozona premalo, je v spodnji plasti ozračja ozona preveč, zaradi česar je škodljiv. Ozon spada v skupino sekundarnih onesnaževalcev zraka in se v urbanem smogu pojavlja kot ena od glavnih sestavin. S problemom škodljivega ozona se srečujemo tudi v Sloveniji. Običajno je ukrepanje ob nastopu prekoračitve opozorilne vrednosti prepozno, saj se maksimalne koncentracije pojavljajo v relativno kratkem časovnem intervalu v času okrog poldneva, ko je sončno sevanje najmočnejše. Obveščanje ljudi o trenutno preseženih koncentracijah ozona ne pomaga bistveno, saj je bila škoda, ki je bila storjena za njihovo zdravje, že v glavnem storjena in je nepopravljiva. Za modeliranje je bila izbrana lokacija Nova Gorica, ki predstavlja tipično urbano okolje z visokimi koncentracijami ozona v poletnih mesecih. Za izgradnjo modela so bili uporabljeni in preizkušeni trije pristopi modeliranja z umetno nevronsko mrežo, Gaussovimi procesi in mehko logiko. V magistrskem delu so v drugem poglavju najprej predstavljene kemične lastnosti in fotokemični proces nastajanja in razgradnje ozona v atomosferi. Razumevanje kompleksnosti pri nastajaju in razgradnji ozona tudi razloži, zakaj sem se odločil za uporabo modeliranja s pomočjo modelov kot so nevronske mreže, Gaussovi procesi in mehka logika, namesto kakšnega drugega fizikalnega ali kemijskega modela. Sledi splošen pregled glavnih skupin povzročiteljev ozona kot so industrija, električne naprave, izpuhi motornih vozil, hlapi goriva in kemična topila, ki so glavni viri dušikovih oksidov (NOX) in hlapnih organskih spojin (VOC). Na kratko je opisan tudi vpliv povišanih koncentracij ozona na ljudi, kot ga opisuje svetovna zdravstvena organizacija. Predstavitvi problematike ozona sledi pregled dela v tretjem poglavju, ki je bilo narejeno na področju modeliranja onesnaženja ozračja. Pregled se osredotoča predvsem na modeliranje kompleksnih pojavov v ozračju, ki jih ni mogoče opisati s preprostimi fizikalnimi in kemijskimi enačbami. Umetne nevronske mreže in mehka logika so bile že uporabljene v modelih za kratkoročno napovedovanje visokih koncentracij škodljivih snovi v ozračju. V literaturi pa ni bilo zaslediti uporabe Gaussovih procesov za podobno nalogo modeliranja. Četrto poglavje zajema metode, ki sem jih uporabil za sintezo modela z umetnimi nevronskimi mrežami, Gaussovimi procesi in mehko logiko. Posebna pozornost je namenjena Gaussovim procesom, ki predstavljajo eno izmed novejših metod za modeliranje. Poleg tega pa so v tem poglavju predstavljene tudi metode za vrednotenje modelov. Pred začetkom izgradnje modela je v petem poglavju najprej podana analiza izmerjenih vrednosti okoljskih parametrov na avtomatski okoljski merilni postaji v iii Novi Gorici, med drugim tudi glede na priporočila iz literature. Opredelil sem tudi geografski položaj Nove Gorice in glavne vire onesnaževanja ozračja. Analizo sem zaključil s predstavitvijo vhodnih spremenljivk modela. Za napovedovanje povišanih koncentracij ozona sem v naslednjem poglavju glede na izbrane vhodne in izhodne spremenljivke zgradil tri različne modele z uporabo treh različnih metod modeliranja opisanih predhodno: z umetno nevronsko mrežo, Gaussovimi procesi in mehko logiko. Za izgradnjo vseh treh modelov je bila uporabljena ista učna množica. Vse tri modele sem po sintezi tudi ovrednotil z uporabo različnih cenilk, kar je podano v sedmem poglavju. Napovedi povišanih koncentracij sem ocenil z uporabo štirih različnih cenilk: korelacijskega koeficienta (r), kvadratnega korena povprečne kvadratne napake (RMSE), cenilko p6 in parametri sposobnosti (SP, SR in SI). Rezultati primerjave vrednotenja modelov kažejo, da je mogoče z različnimi pristopi modeliranja zgraditi povsem enakovredne modele. Prvo prednost modela z Gaussovimi procesi v primerjavi modeloma z umetno nevronsko mrežo in mehko logiko vidim v tem, da je predikcija izhoda modela na osnovi Gaussovih procesov normalna verjetnostna porazdelitev s srednjo vrednostjo in varianco. Drugo prednost pa v optimizaciji precej manjšega števila hiperparametrov, kot v primeru umetne nevronske mreže in mehke logike. Slabost, ki se je izkazala pri optimizaciji modela na osnovi Gaussovih procesov, pa je nastopala v primeru neustreznih začetnih vrednosti hiperparametrov. Izkušnje pridobljene pri izdelavi naloge kažejo, da je dobra identifikacija modela odvisna predvsem od pravilne izbire vhodnih spremenljivk. Glede na pregled priporočil iz literature se je izkazalo, da imajo nekatere vhodne spremenljivke bistveno manjši vpliv na izhod, kot je bilo omenjeno, ali pa celo onemogočajo dobro identifikacijo. Nova Gorica se bo kot urbano okolje v prihodnosti verjetno še širila, kar morda pomeni povečevanje števila prekoračenih vrednosti koncentracij ozona. Ker za vpliv na predvidene povišane koncentracije ozona v urbanem okolju zaenkrat še nimamo direktnih vzvodov, se kaže uporaba modela za napovedovanje povišanih koncentracij ozona edino smiselna in upravičena za varovanje zdravja ljudi. iv Abstract Aim of the thesis is to build a predictive model for high ozone concentrations. Predicted high ozone concentrations that exceed warning and alarming limits would be used to inform public in advance. General environmental aspects of ozone are given in introductory chapter. Countrary to highest layers of athmosphere where the concentration is too low the lower layer of athmosphere contains too high levels of ozone concentration, making it harmful for our health. Ozone is declared as secondary air pollutant. It occurs as one of the main components of urban smog. Slovenia is also confronted with the harmful ozone problem. Usually, it is too late for measures to be taken after ozone concentrations have already exceed warning (limiting) levels, because high ozone concentrations occur in relative short time intervals around noon time when solar radiation is the highest. Informing the public about exceeded ozone concentrations does not help much, because damage that inflicted people health has already been made and their health cannot be restored. Nova Gorica town was used as a modelling location because is represents typical urban environment with high ozone concentrations in summer months. In order to build a model three different modelling methods were used and tested: artificial neural networks, Gaussian processes and fuzzy logic. The thesis in second chapter begins with a description of ozone chemical properties and a photochemical process of ozone composition and decomposition in atmosphere. Understanding this complex process explains my decision to build a model based on modelling techniques such as artificial neural networks, Gaussian processes and fuzzy logic, instead of using some other chemical or physical model. The thesis proceeds with a general overview of main origins of ozone, such as industrial fuel combustion, electrical appliances, motor vehicles and consumer solvents; they represent the main sources of nitrogen oxides (NOX) and volatile organic compounds (VOC). The influence of high ozone concentrations on human health is also briefly described. The thesis continues with an overview of air pollution modelling research field. It mainly focuses on modelling complex occurrences in atmosphere that cannot be described by simple physical and chemical equations. Artificial neural networks and fuzzy logic have already been used to model high concentrations of harmful compounds in atmosphere. However, I did not trace any example of using Gaussian processes for this task. The fourth chapter is covering the methods that were used for model synthesis: artificial neural networks, Gaussian processes and fuzzy logic. The representation in particular focuses attention on Gaussian processes. Before the model synthesis are in chapter five given all measured environmental data from automatic environmental measuring station located in Nova Gorica. Geographical position of Nova Gorica town and its main sources of ozone were v also defined. The model inputs analysis was conducted. According to selected input variables and output variables were in chapter five created three models for high ozone concentration prediction based on different modelling approaches previously described: artificial neural networks, Gaussian processes and fuzzy logic. The construction of all three models applied to one identical learning set of data. In the seventh chapter are models also evaluated with four different evaluation methods: correlation coefficient (r), square root of mean square error (RMSE), p6 and skill parameters (SP, SR in SI). The comparison of evaluation results shows that is possible to create very similar models by applying different modelling methods. Output of Gaussian processes is given as Gaussian conditional distribution that consists of mean and standard deviation of prediction. This is one of the advantages of Gaussian processes, namely the output of the other two modelling techniques is only a means of prediction. The second advantage is small number of hyperparameters that are used during optimization process. A disadvantage of Gaussian processes modelling technique showed numeric instability of the method occuring in a case of inappropriate initial hyperparameter values. Experiences gained in this thesis is confirming that is good model identification basically based on proper selection of input variables. According to the review of recommendations found in literature some of input variables have insignificant influence on output. In the worst case they obstruct good model identification. In future Nova Gorica town as an urban environment will probably expand furthermore, meaning that the number of high ozone concentrations will also rise. The use of high ozone concentration prediction model seems to me only feasible and justified move for human health protection because, no other alternative seems to be currently available. vi Vsebina 1.Uvod ........................................................................................................... 1 2.Ozon ........................................................................................................... 5 2.1.Kemična struktura ozona ....................................................................................... 5 2.2.Ozon kot oksidant ................................................................................................... 6 2.2.1.Električna redukcija ozona ............................................................................................. 6 2.2.2.Kemična redukcija ozona ................................................................................................ 6 2.3.Reagiranje ozona .................................................................................................... 7 2.4.Fotokemični proces nastajanja in razgradnje ozona v atmosferi ..................... 7 2.4.1.Nastajanje in razgradnja ozona v naravi ....................................................................... 8 2.4.2.Nastajanje in razgradnja ozona v urbanem okolju ......................................................... 8 2.5.Povzročitelji ozona v urbanem okolju .................................................................. 8 2.5.1.NOx ................................................................................................................................. 9 2.5.2.Hlapne organske spojine (VOC) ................................................................................... 10 2.6.Vpliv ozona na ljudi .............................................................................................. 11 2.7.Meritve ozona v Sloveniji ..................................................................................... 12 3.Modeliranje onesnaženja ozračja .......................................................... 13 3.1.Uporaba modelov za modeliranje ozračja ......................................................... 16 4.Metode za modeliranje in vrednotenje modelov .................................. 21 4.1.Umetne nevronske mreže .................................................................................... 21 4.2.Gaussovi procesi .................................................................................................. 24 4.2.1.Bayesovo modeliranje .................................................................................................. 24 4.2.2.Modeliranje z Gaussovimi procesi ............................................................................... 26 4.2.3.Kovariančna funkcija .................................................................................................... 32 4.2.4.Stacionarne kovariančne funkcije ................................................................................. 32 4.2.5.Nestacionarne kovariančne funkcije ............................................................................. 35 4.2.5.1.Linearni izraz ......................................................................................................... 35 4.2.5.2.Vhodno odvisen model šuma ................................................................................ 36 4.2.5.3.Prostorsko spremenljiva dolžinska merila ............................................................. 36 4.2.6.Določanje hiperparametrov kovariančne funkcije ....................................................... 36 4.2.6.1.Metoda največje podobnosti .................................................................................. 37 4.2.6.2.Metoda z uporabo Monte Carlo pristopa ............................................................... 37 4.2.6.3.Primerjava opisanih metod .................................................................................... 38 4.3.Mehka logika ......................................................................................................... 39 4.3.1.Mehka identifikacija sistemov ....................................................................................... 40 4.3.2.Grupiranje mehkih množic ............................................................................................ 42 4.3.3.ANFIS (Adaptive Neuro-Fuzzy Inference System) ........................................................ 43 4.4.Metode za vrednotenje modelov ......................................................................... 43 4.4.1.Korelacijski koeficient .................................................................................................. 44 4.4.2.Kvadratni koren povprečne kvadratne napake (RMSE) ............................................... 44 4.4.3.Cenilka p6 ..................................................................................................................... 44 4.4.4.Parametri sposobnosti .................................................................................................. 45 5.Analiza izmerjenih podatkov ................................................................. 47 5.1.Geografski položaj Nove Gorice ......................................................................... 47 vii 5.1.1.Vipavska dolina in Goriško polje .................................................................................. 47 5.1.2.Nova Gorica .................................................................................................................. 48 5.2.Meritve meteoroloških podatkov in koncentracij snovi v zunanjem zraku v Novi Gorici ................................................................................................................... 51 5.3.Izbira vhodnih spremenljivk modela .................................................................. 56 6.Sinteza modelov za napovedovanje povišanih koncentracij ozona. .63 6.1.Sinteza modela z umetno nevronsko mrežo ..................................................... 63 6.2.Sinteza modela z Gaussovimi procesi ............................................................... 65 6.3.Sinteza modela z mehko logiko .......................................................................... 67 7.Vrednotenje modelov ............................................................................. 69 7.1.Preizkušanje modela na osnovi umetne nevronske mreže ............................. 70 7.2.Preizkušanje modela na osnovi Gaussovih procesov ..................................... 71 7.3.Preizkušanje modela na osnovi mehke logike .................................................. 72 7.4.Primerjava rezultatov preizkušanja modelov .................................................... 73 8.Zaključek ................................................................................................. 75 9.Viri in literatura ....................................................................................... 79 Dodatek A: Osnove verjetnostnega modeliranja .................................... 85 A.1. Naključna spremenljivka .................................................................................... 85 A.2. Analiza porazdelitev naključnih spremenljivk ................................................. 86 A.3. Normalna (Gaussova) porazdelitev .................................................................. 86 A.4. Naključni vektorji ............................................................................................... 87 Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi ................................................................................... 91 B.1. Identifikacija modela .......................................................................................... 91 B.2. Prikaz vpliva posameznih hiperparametrov modela ....................................... 98 Dodatek C: Odvisnost negotovosti od kovariančne funkcije .............. 103 C.1. Uporaba stacionarne kovariančne funkcije brez upoštevanja šuma .......... 103 C.2. Uporaba stacionarne kovariančne funkcije z upoštevanjem šuma ............ 104 viii Seznam slik Slika 1: Ozon v atmosferi .................................................................................................................. 1 Slika 2: Resonačni strukturi ozona .................................................................................................... 5 Slika 3: "Opazovana" struktura ozona ............................................................................................... 5 Slika 4: Redukcija ozona ................................................................................................................... 6 Slika 5: Idealizirana reakcija etilena z ozonom ................................................................................. 7 Slika 6: Viri NOx ............................................................................................................................... 9 Slika 7: Viri VOC .............................................................................................................................. 9 Slika 8: Izmerjene vrednosti avtomatske merilne postaje v Novi Gorici dne 05.07.2004 .............. 10 Slika 9: Večnivojska preceptronska mreža ...................................................................................... 22 Slika 10: Nevron .............................................................................................................................. 22 Slika 11: Log-sigmoidna prenosna funkcija .................................................................................... 23 Slika 12: Tan-sigmoidna prenosna funkcija .................................................................................... 23 Slika 13: Linearna prenosna funkcija .............................................................................................. 23 Slika 14: Primer Gaussovega procesa .............................................................................................. 27 Slika 15: Pogojna porazdelitev ........................................................................................................ 28 Slika 16: Porazdelitev v novi točki .................................................................................................. 32 Slika 17: Diagram poteka mehkega sklepanja v mehkih identifikacijskih modelih s pomočjo praktičnega primera “Dinner for two” iz programskega orodja Matlab [54] .................................. 42 Slika 18: Izsek iz državne pregledne karte Slovenije prikazuje širšo okolico Vipavske doline in Goriškega polja ................................................................................................................................ 47 Slika 19: Relief na področju Vipavske doline in Goriškega polja .................................................. 48 Slika 20: Zemljevid Nove Gorice .................................................................................................... 49 Slika 21: Letni potek mesečnih povprečnih temperatur zraka v večjih krajih v Sloveniji za obdobje od leta 1991 do leta 2000 [3] ........................................................................................................... 50 Slika 22: Letni potek mesečnega povprečnega trajanja sončnega obsevanja v urah v večjih krajih v Sloveniji za obdobje od leta 1991 do leta 2000 [3] ......................................................................... 51 Slika 23: Histogram temperature zraka ........................................................................................... 54 Slika 24: Histogram relativne vlage ................................................................................................ 54 Slika 25: Histogram moči globalnega sončnega sevanja na enoto površine ................................... 54 Slika 26: Histogram zračnega pritiska ............................................................................................. 54 Slika 27: Roža vetrov ...................................................................................................................... 54 Slika 28: Statistika rože vetrov ........................................................................................................ 54 Slika 29: Histogram koncentracij žveplovega dioksida .................................................................. 55 Slika 30: Histogram koncentracij ogljikovega monoksida .............................................................. 55 Slika 31: Histogram koncentracij dušikovega monoksida ............................................................... 55 Slika 32: Histogram koncentracij dušikovega dioksida .................................................................. 55 Slika 33: Histogram koncentracij dušikovih oksidov ...................................................................... 55 Slika 34: Histogram koncentracij ozona .......................................................................................... 55 Slika 35: Vhod 1, Temperatura zraka .............................................................................................. 58 Slika 36: Vhod 2, Moč globalnega sončnega sevanja na enoto površine ........................................ 58 Slika 37: Vhod 3, Koncentracija dušikovega monoksida (NO) ...................................................... 58 Slika 38: Vhod 4, Koncentracija dušikovega dioksida (NO2) ........................................................ 59 Slika 39: Vhod 5, Koncentracija ozona (O3) ................................................................................... 59 Slika 40: Vhod 6, Napovedana maksimalna temperatura zraka ...................................................... 59 Slika 41: Vhod 7, Napovedan veter v smeri sever-jug .................................................................... 60 Slika 42: Vhod 8, Napoved vetra v smeri vzhod-zahod .................................................................. 60 ix Slika 43: Izhod, Maksimalne dnevne koncentracije ozona ............................................................. 60 Slika 44: Vrednosti korelacijskih koeficientov med posameznimi vhodnimi spremenljivkami v model in izhodno spremenljivko modela ......................................................................................... 61 Slika 45: Rezultati preizkusa modela z umetno nevronsko mrežo s podatki iz učne množice ....... 64 Slika 46: Histogram napak preizkusa modela z umetno nevronsko mrežo s podatki iz učne množice ............................................................................................................................................ 64 Slika 47: Rezultati preizkusa modela z Gaussovimi procesi s podatki iz učne množice ................ 66 Slika 48: Histogram napak preizkusa modela z Gaussovimi procesi s podatki iz učne množice...66 Slika 49: Rezultati preizkusa modela z mehko logiko s podatki iz učne množice .......................... 68 Slika 50: Histogram napak preizkusa modela z mehko logiko s podatki iz učne množice ............ 68 Slika 51: Rezultat preizkušanja modela na osnovi umetne nevronske mreže ................................. 70 Slika 52: Odstopanje med izmerjenimi in napovedanimi vrednostmi ............................................. 70 Slika 53: Histogram napak preizkušanja modela na osnovi umetne nevronske mreže ................... 70 Slika 54: Rezultat preizkušanja modela na osnovi Gaussovih procesov ......................................... 71 Slika 55: Odstopanje med izmerjenimi in napovedanimi vrednostmi ............................................. 71 Slika 56: Histogram napak preizkušanja modela na osnovi Gaussovih procesov ........................... 71 Slika 57: Rezultat preizkušanja modela na osnovi mehke logike .................................................. 72 Slika 58: Odstopanje med izmerjenimi in napovedanimi vrednostmi ............................................. 72 Slika 59: Histogram napak preizkušanja modela na osnovi mehke logike ..................................... 72 Slika 60: Primerjava korelacijskega faktorja ................................................................................... 73 Slika 61: Primerjava RMSE ............................................................................................................ 73 Slika 62: Primerjava cenilke p6 ....................................................................................................... 73 Slika 63: Primerjava merila pravilnih napovedi SP pri meji 180 mg/m3 ........................................ 73 Slika 64: Primerjava merila uresničenih napovedi SR pri meji 180 mg/m3 ................................... 73 Slika 65: Primerjava indeksa uspeha SI pri meji 180 mg/m3 .......................................................... 73 Slika 66: Primerjava merila pravilnih napovedi SP pri meji 140 mg/m3 ........................................ 73 Slika 67: Primerjava merila uresničenih napovedi SR pri meji 140 mg/m3 ................................... 73 Slika 68: Primerjava indeksa uspeha SI pri meji 140 mg/m3 .......................................................... 73 Slika 69: Primer dvodimenzionalne Gaussove porazdelitvene funkcije ......................................... 89 Slika 70: ?1=1.0, ?2 =1.0, r=0.99 ................................................................................................... 89 Slika 71: ?1=1.0, ?2 =1.0, r=0.00 ................................................................................................... 89 Slika 72 ?1=1.0, ?2 =2.0, r=0.99 ..................................................................................................... 90 Slika 73: ?1=2.0, ?2 =1.0, r=0.99 ................................................................................................... 90 Slika 74: ?1=1.0, ?2 =2.0,r=0.0 ...................................................................................................... 90 Slika 75: ?1=2.0, ?2 =1.0,r=0.0 ...................................................................................................... 90 Slika 76: Prvi model temperature rosišča ........................................................................................ 93 Slika 77: Napaka prvega modela temperature rosišča ..................................................................... 93 Slika 78: Drugi model temperature rosišča ..................................................................................... 94 Slika 79: Napaka drugega modela temperature rosišča ................................................................... 94 Slika 80: Tretji model temperature rosišča ...................................................................................... 95 Slika 81: Napaka tretjega modela temperature rosišča .................................................................... 95 Slika 82: Četrti model temperature rosišča ...................................................................................... 96 Slika 83: Napaka četrtega modela temperature rosišča ................................................................... 96 Slika 84: Peti model temperature rosišča ........................................................................................ 97 Slika 85: Napaka petega modela temperature rosišča ..................................................................... 97 Slika 86: Vpliv prvega hiperparametra na model .......................................................................... 100 Slika 87: Vpliv drugega hiperparametra na model ........................................................................ 101 Slika 88: Vpliv tretjega hiperparametra na model ......................................................................... 102 X Seznam tabel Tabela 4.1: Tabela deležev (contingency table) .............................................................................. 45 Tabela 6.1: Vrednosti hiperparametrov ........................................................................................... 65 xi 1 1.Uvod Vrsto let nas okoljevarstveniki opozarjajo, da v višjih plasteh atmosfere primanjkuje ozona, kar povzroča rast ozonske luknje. V spodnji plasti ozračja pa se vedno bolj srečujemo z obratno težavo. V zraku, ki ga dihamo, je ozona preveč, zaradi česar je škodljiv. Slika 1: Ozon v atmosferi Ozon ima pomemben vpliv na zdravje ljudi in okolje. Glede na lokacijo pojavljanja (slika 1) in koncentracijo je lahko koristen ali škodljiv. Zelo škodljiv vpliv imajo visoke koncentracije ozona, ki se pojavljajo v najnižji plasti atmosfere ob površini Zemlje, troposferi. Ozon spada v skupino sekundarnih onesnaževalcev zraka in škoduje živim bitjem predvsem ob vdihavanju. Škodljiv pa je tudi za poljske pridelke, drevesa in drugo vegetacijo. V urbanem smogu, ki predstavlja zmes prahu, dima in izpušnih plinov, se pojavlja kot ena od glavnih sestavin. Občasne višje koncentracije ozona v troposferi pa imajo lahko tudi zelo koristen vpliv, saj deluje kot dezinfekcijsko sredstvo bolj učinkovito kot klor [30]. Med drugim se uporablja višje koncentracije za sterilizacijo različnih predmetov. Dokazano je bilo, da ozon uniči katerokoli bakterijo, virus, gobo ali plesen pri izpostavljenosti koncentraciji 850 µg/ m3 v štirih minutah. 2 1.Uvod Troposfera sega približno do višine 10.000 metrov, kjer se začenja druga plast atmosfere, to je stratosfera. Stratosfera ali plast koristnega ozona sega od višine 10.000 metrov do 50.000 metrov. V njej se nahaja ozonska plast, ki nas ščiti pred sončnimi ultravijoličnimi žarki in je najgostejša na višini od 18 do 23 kilometrov. S problemom škodljivega ozona se srečujemo tudi v Sloveniji. Glede na publikacije Agencije RS za okolje [2,4] so bile najvišje koncentracije ozona izmerjene v Novi Gorici in na Krvavcu, kjer v poletnih mesecih, ko je na voljo dovolj sončne svetlobe, tudi prekoračijo opozorilne vrednosti ali celo alarmne vrednosti [58]. Običajno je ukrepanje ob nastopu prekoračitve opozorilne vrednosti prepozno, saj se maksimalne koncentracije pojavljajo v relativno kratkem časovnem intervalu v času okrog poldneva, ko je sončno sevanje najmočnejše. Obveščanje ljudi o trenutno preseženih koncentracijah ozona ne pomaga bistveno, saj je bila škoda, ki je bila storjena za njihovo zdravje, že v glavnem storjena in je nepopravljiva. Izgradnja modela za vnaprejšnje opozarjanje ljudi na pričakovane presežene opozorilne in alarmne vrednosti koncentracije ozona je zato postala glavna rdeča nit magistrske naloge. Za modeliranje je bila izbrana lokacija Nova Gorica, ki predstavlja tipično urbano okolje z visokimi koncentracijami ozona v poletnih mesecih. Lokacija na Krvavcu ni enako zanimiva za vnaprejšnje opozarjanje, saj je gostota prebivalstva relativno redka. Napovedovanje preseženih opozorilnih in alarmnih vrednosti koncentracij ozona omogoča lokalnim oblastem vnaprejšnje opozarjanje ljudi in izvajanje določenih ukrepov za njihovo zaščito. Ljudje lahko ob opozorilu spremenijo svoje načrte za prihodnji dan tako, da ne izvajajo težjih telesnih aktivnosti v času visokih koncentracij in se tudi po možnosti čim manj gibljejo zunaj. Lahko tudi zjutraj dobro prezračijo stanovanja in čez dan ne puščajo odprtih oken, po možnosti pa se celo odmaknejo iz območja onesnaženja. Kemične lastnosti in fotokemični proces nastajanja in razgradnje ozona v atomosferi sta podrobneje opisana v drugem poglavju. Razumevanje kompleksnosti pri nastajaju in razgradnji ozona tudi razloži, zakaj sem se odločil za uporabo modeliranja s pomočjo nevronskih mrež, mehke logike in Gaussovih procesov namesto kakšnega drugega fizikalnega ali kemijskega modela. V tem poglavju so predstavljeni tudi glavni viri ozona v urbanem okolju, vpliv ozona na ljudi ter na kratko meritve ozona 1.Uvod 3 v Sloveniji. Umetne nevronske mreže in mehka logika so bile že uporabljene v modelih za kratkoročno napovedovanje visokih koncentracij škodljivih snovi v ozračju. V literaturi pa ni bilo zaslediti uporabe Gaussovih procesov za podobno nalogo modeliranja. Pregled dela, ki je bilo narejeno na tem področju, je podano v tretjem poglavju. Sledi četrto poglavje, kjer so tudi podrobneje predstavljeni in opisani omenjeni modeli. Pri opisu omenjenih modelov sem posebno pozornost namenil predstavitvi Gaussovih procesov. V petem poglavju sem izdelal analizo izmerjenih vrednosti okoljskih parametrov na avtomatski okoljski merilni postaji v Novi Gorici, med drugim tudi glede na priporočila iz literature. V poglavju je tudi opredeljen geografski položaj Nove Gorice in glavni viri oneznaževanja ozračja. Poglavje sem zaključil s predstavitvijo vhodnih spremenljivk modela. Sledi šesto poglavje, kjer je predstavljena sinteza modela napovedovanja povišanih koncentracij iz izmerjenih podatkov s pomočjo treh različnih metod modeliranja: umetne nevronske mreže, mehke logike in Gaussovih procesov. Vrednotenje izgrajenih modelov je podana v sedmem poglavju. Napovedi povišanih koncentracij sem ocenil z uporabo štirih različnih cenilk, ki so tudi podrobneje predstavljene: korelacijski koeficient (r), kvadratni koren povprečne kvadratne napake (RMSE), cenilko p6 in parametri sposobnosti (SP, SR in SI). V zaključku je podana splošna primerjava med uporabljenimi tremi načini modeliranja, izpostavljene so glavne ugotovitve ter postavljene smernice za nadaljnje raziskovanje. 5 2.Ozon 2.1.Kemična struktura ozona Ozon se pojavlja kot naravna oblika kisika [32], [3]. Ime ozon izhaja iz grške besede “ozein”, kar prevedeno pomeni duh ali vonj. Prvič je bil odkrit v laboratorijskih eksperimentih v srednini 19. stoletja. Njegovo kemično strukturo predstavlja slika 2. Slika 2: Resonačni strukturi ozona Iz ilustracije je razvidno, da je vsaka resonančna struktura sestavljena iz ene enojne vezi in ene dvojne vezi. Enojna vez je analogna peroksidnim vezem, ki so relativno šibke in vodijo k tvorbi prostih radikalov. Dvojna vez je analogna molekularnemu kisiku (O2), ki je močno povezan in relativno nereaktiven. Pretvarjanje med ilustriranima resonančnima strukturama ozona je tako hitro, da je struktura “opazovanega” ozona mešanica obeh resonančnih struktur. Posledično temu je moč vezi med osrednjim atomom kisika in preostalima atomoma kisika enaka in znaša eno in pol vezi (1.5), kar je prikazano na sliki 3. Slika 3: "Opazovana" struktura ozona Čisti ozon v obliki plina je modrikaste barve, utekočinjen (-112 °C) je modrikasto-črn, v trdni obliki (-193 °C) pa vijolično-črne barve. 6 2.2.Ozon kot oksidant 2.2.Ozon kot oksidant Razumevanje koncepta oksidacijskih stanj nam omogoča razumevanje zakaj je ozon oksidant. Oksidacijsko stanje določa, koliko elektronov je izgubil (v primeru redukcije pa pridobil) atom v molekuli. Oksidacijsko stanje obeh vodikovih atomov v molekuli vode je +1 (ker vsak atom vodika formalno preda en elektron atomu kisika), oksidacijsko stanje atoma kisika v molekuli vode pa je -2 (ker je atom kisika formalno pridobil po en elektron od vsakega atoma vodika). Oksidacijsko stanje kisikovih atomov je ponavadi -2. V molekuli ozona in molekuli kisika pa imajo vsi kisikovi atomi oksidacijsko stanje 0. Iz tega sledi, da sta tako ozon kot kisik oksidanta, ki znata dobro pridobivati manjkajoče elektrone iz različnih virov, s čimer se poveča oksidacijsko stanje vsaj enega od kisikovih atomov v procesu. V praksi je ozon dejansko močnejši oksidant kot kisik. Delno tudi zato, ker lahko ozon direktno reagira z nekaterimi spojinami, medtem, ko kisik za oksidacijo običajno potrebuje še kakšen katalizator za začetek reakcije. 2.2.1.Električna redukcija ozona Produkt električno pospešene razgradnje ozona je nastanek molekularnega kiska (O2) in atoma kisika z oksidacijskim stanjem -2 (običajno opisan v obliki molekule vode na sliki 4). Enosmerni proces vsebuje tudi redukcijski električni potencial vrednosti 2.07V (relativno na vodikovo elektrodo). Ta vrednost je večja od redukcijskih potencialov večine drugih snovi. To pomeni, da je zmožnost ozona, da oksidira večino drugih snovi, termodinamično mogoča. 2H + + 2e-+O3 - O2+H2O E0 = 2.7V Slika 4: Redukcija ozona 2.2.2. Kemična redukcija ozona Kemično vzpodbujena razgradnja ozona ni povsem preprosta zaradi nastopanja velikega števila različnih reakcij. Eden od bolj preprostih primerov uporabe ozona kot oksidanta v kemičnem procesu predstavlja reakcija ozona z etilenom na sliki 5. V tem procesu nastane pri reakciji ozona z etilenom molekulska struktura A, ki vsebuje dve povezavi kisik-ogljik. V strukturi A izgine tudi dvojna vez ogljik-ogljik. Ogljikovi atomi postanejo oksidirani, kisikovi atomi pa pri povezavi z ogljikom spremenijo oksidacijsko stanje iz 0 na -1. Reakcija se nadaljuje, njen rezultat pa je molekula kisika O2 in molekulska struktura B. Kisikov atom in struktura B sta sedaj v oksidacijskem stanju -2. Substrat (etilen) 2.2.2.Kemična redukcija ozona 7 je postal oksidiran in je en korak naprej do popolne oksidacije (v ogljikov dioksid in vodo). Slika 5: Idealizirana reakcija etilena z ozonom 2.3.Reagiranje ozona Ozon zelo rad reagira z molekulskimi strukturami, ki vsebujejo veliko število vezi (kot so C=C, N=N, C=N, itd), preko mehanizmov, ki so podobni prikazanemu na sliki 5. Niti približno v taki meri pa ne reagira rad s snovmi z eno vezjo kot so C-C, C-O, O-H. Delni vzrok je tudi v tem, da ne obstaja kemična reakcija, ki bi omogočala oksidacijo. Ozon pa dobro reagira s preprostimi oksidizabilnimi ioni kot je S2-, pri čemer nastajajo oksidni anioni kot so SO32- and SO42-. Omenjena oksidacija je preprosta in potrebuje samo kontakt iona z ozonom. Posledično nastaja oksidacija teh ionov z ozonom zelo hitro. Glede na predstavljeno je včasih boljše gledati na ozon kot na zelo reaktivno snov sposobno reagirati z velikim številom drugih snovi, kot pa na snov z zelo visokim redukcijskim potencialom. Razlog je v tem, da so praktične akcije ozona pogosto bolj odvisne od tega kako ozon reagira z ostalimi onesnaževalnimi snovmi, kot pa od njegove zmožnosti, da preprosto pridobiva elektrone (razen v primeru preproste oksidacije ionov). Tako razmišljanje poudarja, da obstaja neka logična reakcijska pot, da lahko ozon reagira z drugimi snovmi. To pomeni, da čeprav ima termodinamična lastnost ozona, da hitro oksidira, prednost, lahko v veliko primerih kinetični faktorji določajo ali bo ozon lahko oksidiral določeno onesnaževalno snov v določenem času. 2.4.Fotokemični proces nastajanja in razgradnje ozona v atmosferi Nastajanje ozona v ozračju je zelo zapleten proces katerega osnova je fotokemična reakcija, ki potrebuje svetlobo [41]. Brez svetlobe poteka proces nastajanja ozona zelo počasi. 8 2.4.1.Nastajanje in razgradnja ozona v naravi 2.4.1.Nastajanje in razgradnja ozona v naravi V naravi nastajo ozon kot posledica obsevanja segretega kisika O2 z ultravijolično svetlobo. Reakcijo prikazuje enačba 3O2+hv^2O3 (2.1). Ob ohlajanju ozračja in prenehanju obsevanja z ultravijolično svetlobo v nočnih urah se začne razgradnja ozona, kot ga opisuje enačba 2O3^3O2 (2.2). 2.4.2. Nastajanje in razgradnja ozona v urbanem okolju Fotokemični oksidanti spadajo v skupino sekundarnih osnaževalcev ozračja. Nastajajo kot posledica sončnega obsevanja atmosfere onesnažene s hlapnimi organskimi spojinami (VOC - Volatile Organic Compounds) in dušikovimi oksidi (NOx). S kemičnega vidika je sinteza fotokemičnega smoga zelo kompleksna in je še v fazi raziskovanja. Splošno znano je, da ogljikovodiki reagirajo z atomi kisika, z ozonom in z oksidacijskimi radikali zaradi povečanja prostega ravnotežja radikalov. Ti prosti radikali reagirajo z ostalimi ogljikovodiki in drugimi oksidnimi vrstami ter dušikovimi oksidi, kar povzroča smog. To je skupek zmesi kot so dušikovi oksidi, ogljikov monoksid, ogljikovodiki, aldehidi, ozon in različne organske komponente. Dušikovi oksidi se štejejo kot ena poglavitnih komponent, ki vplivajo na nastanek fotokemičnega smoga. Pri fotolizni reakciji z valovno dolžno svetlobe, ki je manjša kot 280 nm, se sproščata kisikov radikal in dušikov oksid. Prosti radikali omogočajo nastajanje sekundarnih oksidantov v atmosferi. V takem smislu hidrooksidni in hidroperoksidni radikali povzročajo oksidacijo ogljikovodikov in s tem pomembno verigo reakcij v atomosferi. Čeprav v fotokemični oksidaciji nastopa več različnih komponent, ki povzročajo reakcijo, v splošnem ni mogoče glavnih reakcijskih procesov opisati s preprostim modelom. Iz tega razloga tudi ni mogoče visokih koncentracij ozona v urbanih predelih opisati z reakcijskim procesom, ki bi upošteval samo reakcijo med NO in NO2. 2.5.Povzročitelji ozona v urbanem okolju Glavne skupine povzročiteljev predstavljajo industrija, električne naprave, izpuhi motornih vozil, hlapi goriva in kemična topila, ki so glavni viri dušikovih oksidov (NOX) in hlapnih organskih spojin (VOC) [20]. Približen odstotek posameznega onesnaževalca prikazujeta slika 6 in slika 7. Koncentracije ozona se povečajo predvsem v poletnih 2.5.Povzročitelji ozona v urbanem okolju 9 mesecih, ko je veliko sončne svetlobe in so temperature zraka visoke. Slika 6: Viri NOx Slika 7: Viri VOC 2.5.1.NOx Izmed vseh znanih dušikovih spojin (N2O, NO, NO2, NO3, N2O2, N2O4, N2O5) se pojavljata kot pomembna onesnaževalca zraka samo dušikov monoksid (NO) in dušikov dioksid (NO2). Izmerjena koncentracija obeh plinov se pogosto tudi sešteje in označi z NOX (dušikovi oksidi). Dušikov monoksid je stabilna naravna dušikova snov, ki je zelo pomembna v kemiji stratosfere. Preostali dušikovi oksidi pa igrajo predvsem pomembno vlogo posrednika v mehanizmih fotokemičnega smoga [34]. 10 2.5.1.NOx Dušikov monoksid (NO) je brezbarvna zmes brez vonja, ni gorljiv in je slabo topen v vodi. Spada v skupino zelo toksičnih plinov. V naravi nastaja NO kot produkt bakterijskega delovanja in procesov izgorevanja, ki se pojavljajo naravno. Dušikov dioksid (NO2) je plin rdečkasto-oranžno-rjave barve z ostrim dražljivim vonjem. Spada v skupino zelo toksičnih plinov in je zelo jedek in koroziven. Molekule plina absorbirajo večino svetlobe v vidnem spektru, zato se lahko v atmosferi pojavljajo rumenkaste ali oranžne meglice dušikovega dioksida. V zraku nastaja predvsem kot posledica kemične reakcije pretvorbe dušikovega monoksida v dušikov dioksid in ozon. Pretvorba poteka počasneje v čistem zraku in zelo hitro onesnaženem zraku. v 150 140 130 120 110 100 90 80 70 60 50 40 30 20 10 0 1000 900 800 700 600 h 500 400 300 h 200 100 0 012345678 9 1011121314151617181920212223 Ura dneva Slika 8: Izmerjene vrednosti avtomatske merilne postaje v Novi Gorici dne 05.07.2004 NO in NO2 koncentracije v urbanih okoljih se pojavljajo v ponavljajočem se vzorcu. V jutranjih urah se najprej pojavijo visoke koncentracije dušikovega monoksida NO kot posledica zgodnjega jutranjega prometa. Običajen primer opisanega procesa v poletnih dneh prikazuje slika 8. Na ilustraciji so prikazane izmerjene vrednosti avtomatske merilne postaje v Novi Gorici dne 05.07.2004. 2.5.2.Hlapne organske spojine (VOC) Hlapne organske spojine (VOC) [34,62] so kemične zmesi, ki imajo visok hlapni pritisk in nizko vodno topljivost. V to skupino spada širok 2.5.2.Hlapne organske spojine (VOC) 11 spekter molekul, ki temeljijo na ogljiku, kot so na primer aldehid, ciklične zmesi in dolgo-verižni ogljik. Splošno znano je, da so drevesa naravni učinkoviti viri hlapnih organskih spojin. V vsakdanjem življenju se hlapne organske spojine pojavljajo tudi v obliki različnih topil za barve, čistilnih raztopil in kot primesi v gorivu. Pri onesnaževanju zraka nastajajo kot posledica nepopolnega izgorevanja v motornih vozilih ter kot posledica izhlapevanja pri transportu in pretakanju goriva. V primerih onesnaževanja zraka se hlapne organske spojine mnogokrat razdelijo v dve podskupini, metan (CH4) in nemetanske hlapne organske spojine (NMVOC). Metan spada v skupino izredno učinkovitih plinov, ki povzročajo globalno pregrevanje ozračja. Med nemetanskimi hlapanimi spojinami pa sta najbolj znana benzen, ki lahko povzroča levkemijo pri predolgem izpostavljanju, in 1,3-butadilen, ki je pogosto uporabljen v industrijskih procesih. Podobno kot dušikovi oksidi tudi hlapne organske spojine kažejo podoben vsakodnevni vzorec v urbanem okolju [34]. V zgodnjih jutranjih urah ob povečanem prometu koncentracije narastejo. Ob pojavu sončne svetlobe pa se začne koncentracija zmanjševati, ker se porabljajo pri nastajanju fotokemičnih oksidantov. Ob poznem popoldnevu se pojavi drugi manjši maksimum, v noči pa so koncentracije relativno nizke. 2.6.Vpliv ozona na ljudi Količina ozona, ki se nahaja v troposferi, najnižji plasti ozračja, je bistveno manjša od količine stratosferskega ozona, ki nas ščiti pred ultravijoličnimi žarki [19,63]. Na povišanje koncentracije ozona so najbolj občutljivi ljudje z boleznimi dihal ter boleznimi srca in ožilja. Ozon prodre globoko v pljuča, kjer draži sluznico in pljučno tkivo, tako da ovira dihanje. Pri višjih koncentracijah pride do bolečin pri globljem dihanju in do draženja na kašelj. Prav tako draži sluznico v grlu in povzroča pekoč občutek v očeh. Svetovna zdravstvena organizacija priporoča, naj ljudje ne bodo izpostavljeni osemurnim koncentracijam, ki presegajo 120 µg/m3. Za kratkotrajno enourno koncentracijo ozona ugotavljajo, da je izpostavljenost na več kot 180 µg/m3 že neprijetna za občutljive skupine prebivalcev, med katere spadajo predvsem bolniki in starejši ljudje. Pri 240 µg/m3 pa je koncentracija že tako visoka, da lahko nastanejo škodljive posledice. To je tudi predpisana alarmna vrednost pri kateri je potrebno sprejeti nujne ukrepe. 12 2.7.Meritve ozona v Sloveniji 2.7.Meritve ozona v Sloveniji V Sloveniji se količina ozona v zraku meri že od leta 1992, število merilnih mest pa se je povečalo v letu 1997 [5]. S pomočja sredstev iz projekta Phare se je mreža avtomatskih merilnih okoljskih postaj leta 2001 povečala za osem merilnih mest. Trenutno se v okviru merilne mreže Agencije Republike Slovenije za okolje (ARSO) meri ozon na lokacijah: Ljubljana, Maribor, Nova Gorica, Murska Sobota, Trbovlje, Hrastnik, Zagorje, Krvavec, Iskrba (v bližini Kočevja) in Celje. V začetku leta 2005 pa bo zgrajena tudi postaja v Kopru. ARSO pa sprejema podatke o koncentracijah ozona tudi iz merilne mreže TE Šoštanj, TE Trbovlje in TE Brestanica ter iz avtomatskih mobilnih postaj EIMV (Elektro inštitut Milan Vidmar). Pri meritvah ozona izstopata predvsem postaji Nova Gorica in Krvavec, kjer so bile izmerjene večje vrednosti koncentracije ozona kot drugod [2]. Iz tega razloga so bili za izgradnjo modela za napovedovanje povišane vrednosti ozona izbrani meteorološki in imisijski podatki iz avtomatske merilne postaje v Novi Gorici. 13 3.Modeliranje onesnaženja ozračja Glavni namen pregleda je ugotoviti kaj nam moderna znanost do danes ponuja na področju modeliranja onesnaženja ozračja. Pregled se osredotoča predvsem na modeliranje kompleksnih pojavov v ozračju, ki jih ni mogoče opisati s preprostimi fizikalnimi in kemijskimi enačbami. Ob pregledu sem spoznal kakšne modele so uporabljali raziskovalci, s kakšnimi problemi so se srečevali, kakšna so njihova priporočila in nasveti za obvladovanje problemov ter kakšne izboljšave bi bilo še možno uvesti. Uporaba obstoječega znanja in izkušenj mi je omogočila izgradnjo izboljšanega modela za srednjeročno napovedovanje povišanih koncentracij ozona. Za modeliranje procesov, ki se odvijajo v okolju, je na voljo mnogo različnih pristopov, ki lahko nastopajo tudi v kombinaciji. Izbira najbolj primernega pristopa običajo temelji na kompleksnosti problema, ki ga obravnavamo, in na stopnji razumevanja problema. Uporaba numeričnega modela je najbrž najbolj zaželena rešitev ob predpostavki, da je na voljo zadostno velika in kvalitetna množica podatkov, ustrezna računalniška zmogljivost in dobro teoretično poznavanja problema. V splošnem pa velja, da se s povečevanjem kompleksnosti problema zmanjšuje teoretično razumevanje zaradi slabo definiranih interakcij med sistemi, kar pomeni, da je potrebno do problema pristopiti s statističnimi metodami. Tudi modeliranje ozona v ozračju spada zaradi svoje kompleksnosti v slednjo skupino. Uporaba modernejših pristopov z umetnimi nevronskimi mrežami, mehko logiko ali Gaussovimi procesi, se zdi boljša izbira in predstavlja učinkovito alternativo tradicionalnim statističnim metodam. 14 3.Modeliranje onesnaženja ozračja Razširanje snovi oziroma primesi v ozračju označujemo s skupno besedo disperzija, ki se uporablja kot skupna oznaka za dva mehanizma [40,15]. Prvi mehanizem je transport snovi v ozračju, ki se izvaja s premikom zračnih mas (z vetrom). Drugi mehanizem pa predstavlja difuzija, ki povzroča, da se primesi v ozračju razširjajo tudi prečno na smer vetra in se tako razredčujejo. Difuzija je lahko na ravni molekul, ko primesi prehajajo iz mesta z višjo koncentracijo na mesto z nižjo koncentracijo. To dogajanje je v ozračju praktično zanemarljvo v primerjavi s poglavitno drugo obliko – prisilno difuzijo. Poglavitni mehanizem prisilne difuzije so vrtinci oziroma turbulence v ozračju, ki poskrbijo, da se dim iz vira močno razredči. Kako se bodo škodljivi plini razredčili in kam se bodo razširili določajo emisijski parametri in meteorološke razmere. Pri enakih pogojih emisije je torej predvsem od vremenskih razmer odvisno, kako visoke bodo koncentracije. Razgiban relief močno vpliva predvsem na temperaturno in vetrovno polje, ter s tem posredno tudi na koncentracije škodljivih plinov. Modeliranje širjenja snovi v ozračju v grobem lahko razdelimo na dve skupini [40,15]. V prvo skupino sodijo modeli, ki skušajo ponazoriti (rekonstruirati) fizikalno dogajanje. S pomočjo podatkov o emisiji in o meteorološkem dogajanju želijo v obravnavanem območju rekonstruirati koncentracije primesi. V to skupino sodijo Gaussovi modeli [31,27,44, 46], različni “puff” modeli [25] in zapleteni numerični modeli delcev [55]. Omenjeni fizikalni modeli temeljijo na matematičnem opisu procesov v atmosferi pri katerem se rezultati računajo kot posledica vzrokov [15]. Takšni modeli skušajo razrešiti pripadajoče kemične in fizikalne enačbe, ki določajo koncentracije onesnaževalnih snovi, zaradi česar pa potrebujejo tudi natančne podatke o emisiji in meteorologiji na opazovanem geografskem področju. V splošnem pa obstajajo resne omejitve v prostorski in časovni natančnosti podatkov. Dodatno pa nekaterih podatkov in rezultatov ni možno pridobiti od okoljskih agencij ali lokalne industrije. To pa pomeni resen problem pri uporabi fizikalnih modelov, saj so nekateri vhodi v model neznani. Iz tega razloga so v praksi uporabljajo statistični modeli. V drugo skupino pa sodijo modeli, ki za določena značilna mesta v opazovanem območju skušajo napovedati koncentracije za določeno časovno obdobje vnaprej. Pri tem pa ni nujno da je v model eksplicitno vpleten zapleten mehanizem disperzije. V to skupino sodijo predvsem različni statistični modeli [11,13,22], modeli na osnovi umetnih nevronskih mrež [17,64,16,23,24,14,1], mehke logike [21] in Gaussovi 3.Modeliranje onesnaženja ozračja 15 procesi. Glede na časovni interval napovedi se lahko napovedovalni modeli delijo v štiri skupine [7]: • dolgoročni (dva do tri dni vnaprej), • srednjeročni (do 24 ur vnaprej), • kratkoročni (od 6 do 12 ur vnaprej), • super kratkoročni (od 1 do 2 uri vnaprej). Glede na cilj, ki sem si ga zadal v tej magistrski nalogi, sem v nadaljevanju posvetil večjo pozornost drugi skupini modelov: umetnim nevronskim mrežam, mehki logiki in Gaussovim procesom. Omenjeni modeli omogočajo razvoj relativno dobrega modela tudi v primeru izredno kompleksnega procesa, kot je nastajanje ozona [6]. Mešanice NOx in hlapnih organskih snovi (VOC) v prisotnosti sončne svetlobe igrajo najpomembnejšno vlogo v procesu nastajanja ozona. To je bilo prvič dokazano leta 1940 v dolini Los Angeles [28]. Od takrat naprej je bilo veliko naporov vloženih v razumevanje troposferske fotokemije. Ta prizadevanja dodatno otežuje ogromna množica hlapnih organskih spojin v onesnaženem okolju in izredno veliko načinov, kako lahko posamezne spojine reagirajo. Ocenjeno je bilo, da bi eksplicitna obravnava nastajanja ozona vsebovala več kot 20 000 reakcij z upoštevanjem nekaj tisoč reaktantov in produktov [18]. Za omejitev števila reaktantov in produktov je bilo razvitih kar nekaj kemičnih modelov, ki v svojih mehanizmih uporabljajo povprečno od 30 do 50 spojin ter med 80 in 200 reakcij [18]. 16 3.1.Uporaba modelov za modeliranje ozračja 3.1.Uporaba modelov za modeliranje ozračja M. Z. Božnar in P. Mlakar v svojem članku [14] opisujeta dve vrsti nevronskih mrež, ki sta uporabni na področju modeliranja onesnaženja ozračja: večnivojska perceptronska nevronska mreža in Kohonenova nevronska mreža. Obe mreži sta na kratko opisani, poleg pa so podani tudi praktični primeri uporabe obeh. Prvi primer praktične uporabe predstavlja napovedovanje koncentracij žveplovega dioksida (SO2) na kompleksnem geografskem področju s pomočjo večnivojske perceptronske nevronske mreže [16]. Omenjena nevronska mreža je uporabljena tudi v dveh drugih praktičnih primerih: za napovedovanje hitrosti in smeri vetra na kompleksnem geografskem področju in za rekonstrukcijo meritev SODAR-ja1. Meritve SODAR-ja v različnih plasteh nad kotlino so bile rekonstruirane iz talnih meritev vetra v kotlini in na sosednjih hribih. Kohonenova nevronska mreža je uporabno orodje za samo-organizirajoče razvrščanje. Še posebno je uporabna za analizo gruč. V prvem praktičnem primeru uporabe je razloženo, kako je bila Kohonenova nevronska mreža uporabna za iskanje skupin vetrovnih polj na majhnem kompleksnem geografskem področju. Vetrovno polje je bilo predstavljeno s šestimi talnimi meritvami smeri in hitrosti vetra v kotlini in na sosednjih hribih. V študiji je bila za preiskovanje uporabljena velika množica meritev. Kohonenova nevronska mreža je razvrstila množico podatkov v približno trideset skupin. Za ocenitev kvalitete razvrščanja so izrisane tudi rože vetrov za vsako merilno točko po skupinah. Grafična prestavitev izbranih skupin je dokazala izredno učinkovitost razvrščanja. Drugi primer učinkovite uporabe Kohonenove nevronske mreže pa prikazuje uporabo v pripravljalni fazi izgradnje napovedovalnega modela onesnaženja ozračja. V ta namen je bila razvita 1 SODAR (SOnic Detection And Ranging) sistemi se uporabljajo za merjenje na daljavo strukture vertikalne turbulence in profila vetra v nižjih plasteh atmosfere. SODAR sistemi so podobni radarskim sistemom, (radio detection and ranging), razlikujejo se samo v tem, da so namesto radijskih valov uporabljeni za zaznavanje zvočni valovi. Za SODAR sisteme se uporabljajo tudi druga imena kot na primer “sounder”, “echosounder” ali akustični radar. Večina SODAR sistemov deluje tako, da najprej odda zvočni pulz (pisk) in zatem kratek čas posluša vračajoč se signal. V splošnem se v vračajočem se signalu analizirata tako jakost kot Dopplerjev (frekvenčni) premik. S tem je se ugotovi hitrost vetra, smer vetra in turbulenten značaj atmosfere. Profil atomosfere, ki je funkcija višine, se pridobi z analizo vračajočega signala v zaporednih časovnih intervalih, ki sledijo oddaji vsakega zvočnega impulza. Vračajoč se signal, ki je zabeležen ob določenem časovnem zamiku, vsebuje podatke o atmosferi na višini, ki jo je mogoče izračunati iz hitrosti zvoka. Maksimalni doseg običajnih SODAR sistemov se giblje med nekaj sto metrov do nekaj tisoč metrov. Maksimalni doseg se običajno doseže na lokacijah, ki imajo nizek okoliški šum in zmerno do visoko relativno vlažnost. Na puščavskih lokacijah imajo SODAR sistemi omejeno višinsko zmogljivost, ker se zvok razredčuje hitreje v suhem zraku [8]. 3.1.Uporaba modelov za modeliranje ozračja 17 tehnika razvrščanja vzorcev s pomočjo Kohonenove nevronske mreže [15]. Obsežen pregled modelov na področju atmosferskih znanosti nam ponujata avtorja M. W. Gardner in S. R. Dorling v svojem članku [23], kjer se osredotočata predvsem na modeliranje s pomočjo umetnih nevronskih mrež. V članku je prikazana uporabnost umetnih nevronskih mrež, ki postajajo vse bolj uporabna alternativa tradicionalnim statističnim metodam. Predstavljeno je kar nekaj praktičnih primerov uporabe nevronskih mrež za aproksimacijo, napovedovanje in klasifikacijo v atmosferskih znanostih. Omenjena avtorja sta v drugem članku [24] tudi predstavila uporabo večnivojske preceptronske nevronske mreže za aproksimacijo funkcije za boljše modeliranje relacij med podatki. Prikazan je bil primer modeliranja relacij med urnimi koncentracijami ozona in različnimi meteorološkimi spremenljivkami na lokaciji ob obali. Glavni cilj tega projekta je bilo oceniti pomembnost meteorologije pri ugotavljanju urne koncentracije ozona. Nazorno je prikazano, da prejšnji poskusi reševanja tega problema niso bili primerni zaradi nelinearnega značaja problema. Linearna regresija je podcenjevala pomembnost meteorologije. J. Yi in R. Prybutok sta se v svojem projektu [64] lotila napovedovanja koncentracij ozona v industrijskih predelih severne Amerike. Vhod v njun model predstavlja 9 spremenljivk: jutranja koncentracija ozona, maksimalna dnevna temperatura, koncentracije CO2, NO, NO2, NOx, ter hitrost in smer vetra. Rezultati večnivojske nevronske mreže so bili boljši od rezulatov iz regresijske analize, pri čemer so bili uporabljeni isti podatki. Avtorja omenjata tudi, da je uporabljena nevronska mreža prekosila ARIMA (autoregressive integrated moving average) pristop modeliranja, čeprav bi morale biti take primerjave različnih tehnik modeliranja narejene bolj previdno. Za pošteno primerjavo obeh modelov bi bilo potrebno zgraditi oba modela z istimi podatki. V konkretnem primeru pa je bila nevronska mreža učena z vsemi devetimi meteorološkimi spremenljivkami, medtem, ko je bil ARIMA model zgrajen zgolj s časovnim potekom ozona. Comrie je v svojem prispevku [17] naredil primerjavo napovedovanja ozona s pomočjo večnivojske preceptronske mreže in z regresijskimi modeli. Napovedovanje maksimalnih dnevnih urnih koncentracij ozona v različnih ameriških urbanih okoljih je bilo narejeno z uporabo dnevnih povprečnih meteoroloških podatkov, ki so v bližnji relaciji z ozonom: 18 3.1.Uporaba modelov za modeliranje ozračja maksimalna dnevna temperatura, povprečna dnevna temperatura, povprečna dnevna hitrost vetra in dnevna doza sončnega sevanja. Comrie je v svojem zaključku navedel, da pristop z umetnimi nevronskimi mrežami ni doprinesel dramatičnih izboljšav v primerjavi z linearnimi regresijskimi modeli za napovedovanje ozona s pomočjo meteorološkega stanja okolja. Vsi uprabljeni modeli pa so zelo slabo napovedovali epizode visokih koncentracij. Avtor zaključuje, da bi se dalo model še dodatno izboljšati z ustrezno prilagoditvijo ustrezni geografski lokaciji in z uporabo dodatnih vhodnih spremenljivk. Raziskovanja in napovedovanja troposferskega ozona s pomočjo nevronskih mrež sta se lotila tudi S. A. Abdul-Wahab in S. M. Al-Alawi [1]. Model za napovedovanje ozona je zgrajen na osnovi predpostavke, da je koncentracija ozona funkcija meteroloških pogojev in različnih parametrov ozračja. V svojem delu raziskovalca ugotavljata, da predstavljajo nevronske mreže uporabno orodje zaradi sposobnosti učenja iz historičnih podatkov in zmožnosti modeliranja nelinearnosti. Mrežo sta učila na podlagi meteoroloških podatkov in podatkov o koncentracijah2 snovi v zunanjem zraku v poletnem času, ko so bile izmerjene najvišje koncentracije ozona. Vsi podatki so bili izmerjeni v urbanem okolju. Končni produkt njunega raziskovanja so bili trije modeli nevronskih mrež. Glavni namen prvega modela je bilo ugotavljanje parametrov, ki imajo največji vpliv na koncentracije ozona v časovnem intervalu 24 ur (vključeni so bili dnevni in nočni podatki). Drugi model je bil razvit za opazovanje vpliva posamezih parametrov na koncentracije ozona v dnevnem času, ko so bile izmerjene višje koncentracije ozona. Tretji, zadnji model, pa je bil razvit za napovedovanje maksimalnih dnevnih koncentracij ozona. V svojem delu sta ugotovila, da je znašal prispevek meteroloških parametrov na variacijo ozona približno od 33% do 40%. Ugotovila sta tudi, da imajo dušikov oksid, žveplov dioksid, relativna vlažnost, nemetanski ogljikovodiki in dušikov dioksid najpomembnejši vpliv na napovedane koncentracije ozona. Dodatno je pomembno vlogo igrala tudi temperatura zraka, medtem, ko je imelo globalno sončno sevanje bistveno manjši vpliv, kot je bilo pričakovano. G. Finzi in G. Nunnari sta se v svojem delu [21] lotila raziskave različnih pristopov modeliranja z uporabo nevronskih mrež in mehkih modelov. Poleg tega sta tudi ovrednotila rezultate realnih primerov modeliranja z 2 Koncentracija snovi v zunanjem zraku se izraža v masnih enotah (miligramih, mikrogramih) na enoto volumna zraka (m3) pri temperaturi 293 K (20°C) in zračnem tlaku 101.3 kPa [59] 3.1.Uporaba modelov za modeliranje ozračja 19 omenjenimi modeli na praktičnih primerih v urbanih in industrijskih področjih v Italiji. Ocene rezultatov so predstavljene in razložene v obliki urnih in dnevnih indeksov učinkovitosti napovedovanja in statistčnih indikatorjev. Z umetnimi nevronskimi mrežami sta bila zgrajena dva modela za napovedovanje koncentracij ozona od 1 do 6 ur vnaprej v dveh mestih: Brescia (v severni Italiji) in Catania (na Siciliji). Za vhodne parametre modela so bile izbrane izmerjene urne vrednosti ozona in njegovih pozvročiteljev (CO, NO in NO2 koncentracije), medtem ko avtorja poročata, da meteoroloških podatkov ni bilo na voljo. Pri oceni uporabnosti se je izkazalo, da imata modela slabši učinek pri napovedovanju izrednih situacij v primerjavi z modelom sive škatle, ki pa obratno povzroča preveliko število lažnih alarmov. V drugem primeru pa sta bila zgrajena dva modela na osnovi mehke logike za napovedovanje preseženih vrednosti koncentracij ozona v dveh mestih: Brescia in Siracusa (na Siciliji). Za vhodne podatke modela so bile tokrat na voljo izmerjene koncentracije snovi v zunanjem zraku in meteorološki parametri. Končna izbira dnevnih povprečnih vrednosti vhodov je bila sledeča: O3, NO2, NOx, temperatura zraka, sončno sevanje, zračni pritisk in smer vetra. Končna analiza napovedi je pokazala relativno zadovoljivo učinkovitost napovedovanja izrednih dogodkov ter dobro izogibanje lažnim alarmom. Avtorja v svojem delu zaključujeta, da so pristopi modeliranja z uporabo nevronskih mrež in mehke logike zelo obetajoči za napovedovanje izrednih dogodkov, potrebne pa so dodatne raziskave. Ballester in ostali so se v svoji raziskavi [11] lotili napovedovanja urnih koncentracij za 24 ur vnaprej v treh mestih v Španiji: Paterna, Alcoy in Carcagente. Za napovedovanje so uporabili tri vrste modelov: ARMAX (Autoregressive-moving average with exogenous inputs), večnivojsko perceptronsko nevronsko mrežo in FIR nevronsko mrežo [54]. Posebnost večnivojskega perceptronskega modela je v tem, da uporablja za napovedovanje koncentracije ozona za eno uro vnaprej pretekle izmerjene vrednosti in že izračunane napovedi. Za takšno modeliranje dinamike procesa so bili za vhode v model izbrani sledeči regresorji: predhodne vrednosti ozona (O3(t-1), O3(t-2), ..., O3(t-24)), meteoroloških podatkov (hitrost vetra(t-24), smer vetra(t-24), temperatura zraka (t-24), ...) in dušikovega dioksida (NO2(t-24)). Učinkovitost opisanega modela je primerjana z ostalima dvema modeloma, za katera avtorji trdijo, da sta se izkazala slabše, še posebej FIR nevronska mreža, katere glavna omejitev je, da je bila za vhod uporabljena samo koncentracija 20 3.1.Uporaba modelov za modeliranje ozračja ozona brez ostalih meteoroloških podatkov. Pri analizi slabših napovedi so uporabili metodo gruč (cluster analysis). Izkazalo se je da, so bile slabše napovedi povezane z gručami, ki so bile manj razpršene. Za izboljšavo v takih primerih avtorji predlagajo uporabo dodatnih spremenljivk (na primer posebne vremenske elemente). 21 4.Metode za modeliranje in vrednotenje modelov Za izgradnjo modela za napovedovanje povišanih koncentracij ozona se mi zdi uporaba modelov na osnovi umetne nevronske mreže, Gaussovih procesov in mehke logike glede na opisano kompleksnost ozona in predstavljene izkušnje ostalih raziskovalcev ustrezna. V nadaljevanju sem se najprej posvetil predstavitvi metod, ki sem jih uporabil za sintezo modela: umetnim nevronskim mrežam, Gaussovim procesom in mehki logiki. Posebno pozornost sem namenil Gaussovim procesom, ki predstavljajo eno izmed novejših metod za modeliranje. Za primerjavo rezultatov modelov zgrajenih z različnimi metodami pa sem uporabil štiri različne metode določanje uspešnosti modelov na preizkusnih podatkih. 4.1.Umetne nevronske mreže Kot je bilo že omenjeno, predstavljajo umetne nevronske mreže, natančneje večnivojske perceptronske mreže, učinkovito alternativo tradicionalnim statističnim metodam. Dokazano je bilo [29], da je možno uporabiti in naučiti večnivojsko perceptronsko mrežo za aproksimacijo kakršnekoli gladke omejene funkcije. V nasprotju z ostalimi statističnimi metodami ni potrebno pri učenju večnivojske perceptronske mreže upoštevati nobenih predpostavk o verjetnostni porazdelitvi podatkov. Mreža je sposobna modelirati izredno nelinearne funkcije in natančno posploševati pri učenju z novimi podatki. Predhodno je bilo opisanih kar nekaj praktičnih primerov uporabe nevronskih mrež v atmosferski znanosti. Večnivojska perceptronska mreža je sestavljena iz sistema povezanih preprostih nevronov (vozlišč), kot je prikazano na sliki 9 [54]. Vozlišča so razdeljena v več nivojev, običajno na vhodni nivo, poljubno število skritih nivojev in izhodni nivo. V vhodnem nivoju je število vozlišč enako številu vhodov, v izhodnem nivoju pa številu izhodov. V skritih nivojih pa je število nevronov poljubno in se ga običajo določi glede na izkušnje, ki se pridobijo pri procesu učenja nevronske mreže. V posameznih nivojih vozlišča niso med seboj povezana. Vsi izhodi posameznega nivoja so povezani preko uteži z vsemi vhodi naslednega nivoja. Smer toka podatkov tudi določa ime naprej podajajoča nevronska mreža (feedforward neural network). 22 4.1.Umetne nevronske mreže (Ur, VHODNI NIVO R...ŠTEVILO VHODOV S...ŠTEVILO NEVRONOV V SKRITEM NIVOJU vi...UTEŽI IZHODNEGA NIVOJA wi,j...UTEŽI SKRITEGA NIVOJA ci...PRAGI IZHODNEGA NIVOJA bi...PRAGI SKRITEGA NIVOJA Slika 9: Večnivojska preceptronska mreža Nevron (vozlišče) predstavlja osnoven gradnik večnivojske vnaprejšne perceptronske nevronske mreže. Prikazuje ga slika 10. Vsak vhod v nevron je utežen z ustrezno utežjo wi,R. Vsota vseh zmnožkov vhodov z utežmi in vrednostjo praga b tvorijo vhod v prenosno (aktivacijsko) funkcijo f, ki je lahko poljubna odvedljiva funkcija. Uporaba mnogo preprostih nelinearnih prenosnih funkcij v večnivojski vnaprejšni mreži lahko po predpostavki omogoča aproksimacijo izredno nelinearnih funkcij. V primeru uporabe linearnih prenosnih funkcij pa je možno modeliranje samo linearnih funkcij [29]. Običajo se uprabljata za modeliranje log-sigmoidna (slika 11) in tan-sigmoidna (slika 12) prenosni funkciji, ki imata preprosto izračunljiv odvod, občasno pa tudi linearna prenosna funkcija (slika 13). R...ŠTEVILO VHODOV Slika 10: Nevron 4.1.Umetne nevronske mreže 23 Slika 11: Log- Slika 12: Tan- Slika 13: Linearna sigmoidna prenosna sigmoidna prenosna prenosna funkcija funkcija funkcija Eno od ključnih vprašanj za uspešno modeliranje je določitev števila nevronov v skritih plasteh. To vprašanje je v literaturi [35] sicer obdelano, vendar ne daje zadovoljivih odgovorov na splošne primere. Število nevronov sem določal ekperimentalno. Osnova pa je bilo število nevronov, ki jih je predlagal programski paket NeuroShell2 [43], izračunanih po naslednji formuli: m = ^+W) (4.1) Pri čemer je pomen oznak sledeč: m ... število skritih nevronov, n ... število vhodov, q ... število izhodov in q ... število učnih vzorcev. Število skritih nevronov, ki ga izračunamo po enačbi (4.1) je odvisno od števila vzorcev v učni množici s katero učimo nevronsko mrežo. Ker se zavedam, da je obširna množica učnih vzorcev zagotovo redundantna, je torej tudi tako ocenjeno število skritih nevronov precenjeno. Zato sem to število tudi zmanjševal dokler nisem dobil zadovoljivih rezultatov. Navedena formula je sicer dober pripomoček za dosego cilja, se pa zavedam, da mora biti število skritih nevronov odvisno od zapletenosti procesa in ne od števila učnih vzorcev, ki so na voljo, saj se le-ti lahko ponavljajo. V primeru prevelikega števila skritih nevronov se lahko nevronska mreža predobro nauči učne podatke s čimer se zmanjša njena zmožnost posploševanja. Obratno pa lahko pride do prevelikega posploševanja v primeru premajhnega števila skritih nevronov. Po določitvi števila skritih plasti, števila nevronov in začetnih vrednosti uteži in pragov, je nevronska mreža pripravljena za učenje. Učenje 24 4.1.Umetne nevronske mreže zahteva množico učnih podatkov, ki so sestavljeni iz dvojic {vhodni vektor-izhodnih vektor}, kar pomeni, da učenje nevronskih mrež spada v obliko nadzorovanega učenja (supervised learning). Med učenjem so nevronski mreži večkrat predstavljeni učni podatki. Ob vsaki predstavitvi (iteraciji) se uteži spreminja (prilagaja), dokler ni dosežena ustrezna vhodno-izhodna relacija. Proces učenja dejansko predstavlja iterativni postopek, pri katerem nastavimo vrednosti uteži tako, da minimiziramo določeno kriterijsko funkcijo (običajno je to funkcija povprečne kvadratne napake). Za učenje perceptronskih nevronskih mrež je na voljo mnogo algoritmov, najbolj pogosto uporabljen je algoritem vzvratnega širjenja (backpropagation) [40,15,54,50]. O kvaliteti učenja odločata predvsem dva nastavljiva parametra, to sta hitrost učenja in moment. Hitrost učenja pove v kolikšni meri se pri tekočem popravku uteži upoštevajo izračunane napake, moment pa določa, koliko se pri tekočem popravku uteži upoštevajo predhodni popravki uteži. Premajhna hitrost učenja ima za posledico veliko potrato časa, povečana pa je tudi možnost, da bo mreža ostala v lokalnem minimumu, ki ga bo našla. Nasprotno pa prevelika hitrost učenja lahko povzroči oscilacije. Nizke vrednosti momenta pomenijo nevarnost oscilacij med učenjem in počasno učenje. Do oscilacij pa lahko pride tudi v primeru previsoke vrednosti momenta. 4.2.Gaussovi procesi Teorija Gaussovih procesov (GP) izhaja iz področja verjetnosti in statistike, njihova uporaba za identifikacijo nelinearnih dinamičnih sistemov pa predstavlja zelo mlado vejo raziskovanja. Kot osnovo za identifikacijo z GP se uporablja njihova podobnost z umetnimi nevronskimi mrežami. Rezultat predikcije z modelom GP je normalno porazdeljena gostota verjetnosti modelirane izhodne vrednosti procesa s srednjo vrednostjo in varianco. To je tudi njihova prednost, saj v nasprotju z drugimi nelinearnimi modeli, dobimo tudi informacijo o zaupanju v napovedano vrednost izhoda. Osnove verjetnostnega modeliranja povzete po [36], ki so potrebne za razumevanje modeliranja z Gaussovimi procesi, so podane v dodatku A. 4.2.1.Bayesovo modeliranje Najprej definirajmo problem, ki ga želimo rešiti [26,47]. Podane imamo izmerjene podatke D={X,t}={(xn,tn),n=1..N)} (izmerjene vrednosti vsebujejo tudi šum), ki so sestavljeni iz N parov vhodnih vektorjev {xn} 4.2.1.Bayesovo modeliranje 25 (dimenzija vhodnega vektorja je L) in skalarnih izhodov {tj. Iz podanih podatkov želimo sestaviti model s pomočjo katerega bomo lahko izdelovali napovedi. Glede na novi vhodni vektor xN+1 želimo poiskati verjetnostno porazdelitev ustreznega izhoda tN+1. Splošen model za opisan problem predstavlja enačba tn=y(*J+vn (4.2) kjer označuje y(xn) funkcijo modela in v„ šum. V prostoru možnih funkcij, ki bi bile lahko uporabljene za model, definirajmo predpostavljeno verjetnostno porazdelitev p(y\A), kjer je z A označena množica hiperparametrov. Poleg tega definirajmo tudi prepostavko o verjetnostni porazdelitvi šuma p(v\B), kjer v označuje vektor šumnih vrednosti v=(vuv2,...vn) in B označuje množico hiperparametrov modela šuma. Glede na podane hiperparametre lahko z enačbo p(tN\{xn),A,B)=] J p(tN\{xn),y,v)p(y\A)p(v\B)dydv (4.3) — 30 —30 zapišemo verjetnostno porazdelitev podatkov, kjer je tN=(t1}t2,...,tN). Če sedaj definiramo še vektor t=(tut2,...,tN,tN+i), lahko zapišemo pogojno verjetnostno porazdelitev vrednosti tN+1 z enačbo (4.4) in s tem tudi uporabimo pogojno verjetnostno porazdelitev za izdelavo napovedi o vrednosti tN+1 . P(tNjD,A,B,xNJ = ^f^^^ (4.4) Na tej točki je potrebno poudariti, da je potrebno ločiti tN+1 (neznana izhodna vrednost pri vhodnem vektorju xN+1 ) od tN+1 (vektor, ki je sestavljen iz tN in tN+1 ). Z vidika Bayesovega verjetnostnega modeliranja, modeliranje statičnega procesa z nelinearno funkcijo f=f(x,(o), katere parametre co={A,B} določimo glede na podatke Z), lahko predstavimo v obliki Bayesovega teorema [36,26,47]: pHd)=eijMeM (4.5) kjer je: p(w)...“a-prior” verjetnostna porazdelitev, ki vsebuje predhodno znanje o parametrih funkcije (ang. prior), ki ponavadi predpostavljajo zveznost, frekvenčno razporeditev moči itd, p(D ^...verjetnostna porazdelitev učne množice pri danih parametrih 26 4.2.1.Bayesovo modeliranje funkcije (ang. likelihood), p(D)...verjetnostna porazdelitev učne množice (ang. evidence), v primeru Bayesovega modeliranja služi kot normalizacij ska konstanta, /?(ft>|DJ...“a-posteriori” verjetnostna porazdelitev parametrov co pri dani učni množici D (ang. posterior). Bayesov teorem združuje predhodno znanje o parametrih p(co) z znanjem, dobljenim v obliki učne množice D Z upoštevanjem zapisa D={X,t} lahko Bayesov teorem zapišemo v obliki: p^^ ^jm^ (4.6) 4.2.2.Modeliranje z Gaussovimi procesi Z neparametričnimi metodami določamo predikcijo izhoda brez eksplicitne parametrizacije funkcije/ Tak primer je tudi modeliranje z Gaussovimi procesi [26,47,36]. Gaussov proces predstavlja množica naključnih spremenljivk t=(t1(x1),t2(x2),...), ki imajo skupno verjetnostno porazdelitev opisano z enačbo 1 - 1(t-tfc-\t-v) P(tN\C,{xn}) = -e 2 (4.7) za poljubno množico vhodov {x„}, kjer Z predstavlja ustrezno normirno konstanto. V primerjavi z enačbo (A.20) iz dodatka A je kovariančna matrika v tem primeru označena z C, konstatna Z pa definira izraz Z=(V2^)"aH . (4.8) Na sliki 14 je predstavljen primer Gaussovega procesa z množico štirih naključnih spremenljivk t=(t1(x1),t2(x2),t3(x3),t4(x4)). Vsaka naključna spremenljivka t„ prestavlja eno učno točko. Porazdelitev naključne spremenljivke p(tn) v učni točki t„ določa kovariančna matrika C, katere dimenzija je v prikazanem praktičnem primeru enaka 4x4. Preprost in nazoren praktičen primer modeliranja nelinearne funkcije je podan v dodatku B. 4.2.2.Modeliranje z Gaussovimi procesi 27 Slika 14: Primer Gaussovega procesa Kovariančna matrika podatkov C je definirana s parametrično kovariančno funkcijo Cmn(xmxn;0), kjer se označuje Cmn=Cm(xmxn;0), 0 pa označuje množico hiperparametrov. Določanje vrednosti hiperparametrov in dimenzije množice hiperparametrov bo predstavljeno v nadaljevanju. Zaenkrat pa za naš praktičen primer predpostavimo, da imamo že določeno množico hiperparametrov in kovariančno funkcijo. Indeks mn zgoraj označuje, da kovariančna funkcija razlikuje med xm in xn tudi, ko velja xm=xn . To je še posebej pomembno pri obravnavi modelov šuma v okviru kovariančne funkcije. V nadaljevanju bo indeks zaradi preglednosti opuščen, kljub temu pa se nanaša na vse kovariančne funkcije v nadaljevanju. Z // je označen vektor srednjih vrednosti. Glavni razlog, da je množica hiperparametrov označena s & je v tem, da imajo parametri kovariančne funkcije podobno vlogo kot parametri nevronske mreže. Sedaj predpostavimo, da je Gaussov proces s kovariančno matriko CN in vektorjem srednjih vrednosti // ustvaril vektor podatkov tN tako kot opisuje enačba (4.7). Očitno je, da ni možno vseh podatkov modelirati s procesom, ki bi imel srednjo vrednost enako 0. V nadaljevanju bomo predpostavljali, da so bili podatki ustrezno skalirani, tako da zagotavljajo domnevno o srednji vrednosti enaki 0. Z eksplicitno definicijo verjetnostne porazdelitve podatkov smo sedaj preskočili korak izražanja 28 4.2.2.Modeliranje z Gaussovimi procesi individualnih predpostavk o šumu in modelirni funkciji. Obe predpostavki sta združeni v kovariančni matriki C. Iz predstavljenega je razvidno, da je narava kovariančne matrike C in posledično kovariančne funkcije bistvena pri modeliranju z Gaussovimi procesi. Kovariančna matrika bo podrobneje predstavljena v naslednjem podpoglavju. Slika 15: Pogojna porazdelitev Funkcija zapisana z enačbo (4.7) izraža korelacijo med različni točkami v vhodnem prostoru, kar lahko najbolje ilustriramo s preprostim primerom. Upoštevajmo sliko 15. Slika prikazuje gostoto skupne Gaussove verjetnostne porazdelitve dveh točk (naključnih spremenljivk) h in t2, ki ju določa kovariančna matrika C. Pri podani vrednosti U=t lahko najdemo Gaussovo pogojno verjetnostno porazdelitev P(t2\ti=t,C) in iz tega tudi najbolj verjetno napoved vrednosti t2 pri podani vrednosti tj=t. Skala in orientacija gostote skupne Gaussove verjetnostne porazdelitve je popolnoma določena z kovariančno matriko, ki zopet kaže na glavno vlogo C pri modeliranju z Gaussovimi procesi. V splošnem primeru napovedovanja tN+1 s podanimi podatki D=({xn),tn) lahko uporabimo enačbo (4.7) ter dejstvo, da je skupna verjetnostna porazdelitev množice podatkov Gaussova, za izračun Gaussove pogojne verjetnostne porazdelitve tN+1 . P{tN+1\D,C{xm, xn;@), xN+1) = P{tN+1\D,C{xm,xn;@),xN+1) = P(tN+1\C{xm,xn;@),xN+1,{xn\) P(tN\C(xm,xn;0),{xn}) Z Z N e 2 - 1™c»+A+1 -^c^) (4.9) (4.10) JV + 1 ZN in ZN+1 označujeta primerni normirni konstanti. Enačba 4.2.2.Modeliranje z Gaussovimi procesi___________________________________29 prikazuje relacije med CN+1 ter CN, kjer je pomen oznake K = C(xN+1,xN+1;0) (4.12) ter oznake kN+l = {C{xlxN+l;©),...,C{xN,xN+l;©)) . (4.13) Gaussovo pogojno verjetnostno porazdelitev fo+,pri xw+7 lahko izpeljemo iz enačbe (4.10) po postopku opisanem v nadaljevanju. Če določimo inverz matrike CN+1 z enačbo \mn mNJ ~k+1 /i I c-L r (4.14) in uporabimo dejstvo, da je (4.15) lahko zapišemo CnMn+kN+lrnN+l = In (4.16) C^m^+zi^^O (4.17) *J+1MAr + Km;+1 = 0r (4.18) kTN+lmN+l+Kn = 1 (4.19) 30 4.2.2.Modeliranje z Gaussovimi procesi S premultiplikacijo enačbe (4.16) z C"1 in nadomestitvijo dobljenega izraza za M* v enačbi (4.18) ter upoštevanjem C-NT=C~N1 dobimo — C~1k m N+1 =, t r-1,: • (4.20) Če dobljeni izraz vstavimo v enačbo (4.17), dobimo 1 /i=t n-1 . (4.21) K-k N+1CN k N+1 Izraz (4.20) lahko zapišemo z upoštevanjem enačbe (4.21) tudi v drugačni obliki m^^-uC^k^^ml^-ixkl^Cr^-ixkl^ . (4.22) S pomočjo izrazov (4.16) in (4.22) pa lahko določimo tudi izraz za MN=c-N1-c^kN+1ml+1-c-N1+^crN1kN+1{c-N1kN+1)T . (4.23) Če se vrnemo nazaj k enačbi (4.10) lahko zapišemo odvisnost eksponenta enačbe (4.10) z elementi matrike C"^ definirane v enačbi 4.14 Mn mN + 1 tN 7--1 C"1 t -tTC~1t =\tT t \M n m™|'"l '-'JV+1'jV + 1 N^NlN [N lN + l\ j -t C~ t = T t N N N mN+1 V 'JV + 1 = Vtl+l + tN+l(m^N+tTNmN+1)+tTNMNtN-tTNC-N1t = m2N+1 + 2tN+1mTN+1tN+tTNMNtN-tTNC-N1tN (4.24) Z uporabo izrazov (4.22) in (4.23) je nadaljevanje izpeljave (4.24) sledeče tT C1 t -tTC~1t = Ju4+1-)u2riV+1Ä;+1C-;^+^C-1^+)u^C-1^+1(C-1ÄAr+1)r^-^C-1^ -M+,-2tNAkl+1C-;tN)H{C-N1kNJTtN)T{kTN+1C-N1tN))- N + l N + lN+1N NN N+1 N N+1 N N = /i(4+i-2W1(*;+1C],1^)+(*;+1C-1^)7'(*;+1C-1U = = lA(tN+l-(kTN+1C-N1tN))2 = (tN+1-(kTN+1C-N1tN)f KK KN + 1C N tN) (4.25) 4.2.2.Modeliranje z Gaussovimi procesi 31 Če v rezultatu izpeljave (4.25) uporabimo dve novi oznaki tN+1 = kTN+1C~N1tN in (4.26) = K-kT N+1 C N 1k N+1 , (4.27) lahko enačbo za izračun Gaussove pogojne verjetnostne porazdelitve tN+1 preoblikujemo v novo obliko ( (tN+1-îN+1f P(tN+1\D,C(xm,xn;0),xN+1)= 1 e 2 2 1 , (4.28) kjer î„+1 označuje srednjo vrednost napovedi v novi točki, o\+1 pa standardno deviacijo napovedi. Pomembno je tudi poudariti, da obe enačbi vsebujeta samo CN in ne CN+1. Posledično potrebujemo za napovedovanje pri konstatni učni množici ter pri konstantnih hiperparametrih kovariančne funkcije samo en izračun inverzne matrike reda NxN. Glede na primer modeliranja s štirimi učnimi točkami, ki je bil podan na sliki 14 na začetku podpoglavja, je prikazan na sliki 16 primer napovedi v novi točki t5. Kot je nazorno ilustrirano v dodatku C, sledi iz enačbe (4.27), da je negotovost, ki jo določa kovariančna funkcija, v novi napovedani točki večja od negotovosti v ostalih učnih točkah. Negotovost se z oddaljevanjem od učnih točk povečuje, z bližanjem učnim točkam pa zmanjšuje. V nobenem primeru pa negotovost ne more biti manjša, kot je v najbližji učni točki. 32 4.2.2.Modeliranje z Gaussovimi procesi Slika 16: Porazdelitev v novi točki 4.2.3.Kovariančna funkcija Kot je bilo že predhodno omenjeno, je vloga kovariančne funkcije najbolj pomembna pri modeliranju z Gaussovimi procesi [26]. Napovedane verjetnostne porazdelitve, ki nastopajo pri podanih podatkih, so v glavnem odvisne od kovariančne funkcije in njenih hiperparametrov. V nadaljevanju sledi obravnava različnih omejitev pri obliki kovariančnih funkcij in nekateri primeri kovariančnih funkcij za različne probleme. 4.2.4.Stacionarne kovariančne funkcije Za kovariančno matriko CN obstaja samo ena omejitev. Kovariančna matrika CN, ki jo določa kovariančna funkcija C(xmx„;&) in množica vhodov {xX=1 , mora biti pozitivno definitna za katerokoli rang N in katerokoli množico vhodov. Najprej preučimo samo stacionarne kovariančne funkcije, ki so odvisne samo od relativne pozicije dveh vhodnih vektorjev in zadoščajo enačbi C{xitxi+h;0)=C{h;0) . (4.29) Taka kovariančna funkcija C(h;0) zadošča pogoju pozitivne definitnosti, če velja zanjo pozitivna omejena simetrična metrika G(.), kot je predstavljena v enačbi 4.2.4.Stacionarne kovariančne funkcije 33 C(xm-xn;®)stat=C(h;@)sta = J ... J ewThG(dw) (4.30) — 30 —30 Enačba je znana pod imenom Bochnerjev teorem [12]. Enostavno se ga da prikazati z upoštevanjem izraza (4.30) o pozitivni definitnosti za katerikoli realni vektor r. rTCr>0 (4.31) Glede na obliko stacionarne kovariančne funkcije lahko to zapišemo z enačbo TrmrnCsJxm-xn;©)>0 (432) m, n Če sedaj uporabimo v enačbi (4.32) izraz za Cstat(h) iz enačbe (4.30)dobimo enačbo HrmrnCstat(xm-xn;0)^rmrJe wT(x^G(dw) m, n m, n = SOrewTxm)Orne-wTxlG(dw) (4 33) m n = S\TrewTx12G(dw) m Iz tega sledi, da je Cstat(h;0) pozitivno definitna za katerokoli množico vhodov {xm}, če je glede na enačbo (4.30) G(.) pozitivna metrika. Dejansko enačba (4.30) ne omejuje drastično možnosti izbire kovariančnih funkcij. Primer je katerakoli Fourierjeva transformacija f(h) , kjer /(h)>0forallh zadošča zahtevi o pozitivni definitnosti. Za mnogo primerov modeliranja pa želimo, da vsebuje kovariančna funkcija dve lastnosti. 1. Za večino primerov modeliranja pričakujemo, da je korelacija med dvema točkama, ki sta si blizu v vhodnem prostoru (glede na določeno metriko), večja od korelacije med dvema bolj oddaljenima točkama. 2. Pogosto ugotavljamo, da se v podatkih nahaja tudi šum, zato želimo vpeljati v kovariančno funkcijo tudi model šuma. Ker je vsota dveh pozitivno definitnih funkcij zopet pozitivno definitna funkcija, bomo razdelili definicijo kovariančne funkcije na dva dela. Prvi del bo predstavljal obliko funkcijo s katero želimo modelirati naše podatke, drugi del pa bo povezan z modelom šuma CM(xm-xn;&) = Cf(xm-xn;&) + 5mn&3 . (4.34) Z ômn je označena Kronceker-jeva delta funkcija 34 4.2.4.Stacionarne kovariančne funkcije ö I 0; m^n mn i 1 ; m — n (4.35) Eno izmed možnih oblik Cf prikazuje enačba (4.36), kjer x(1n] označuje / -to komponento vektorja x„ . Cf{xm,xn;0) = 01e 2l1 2r' +02 (4.36) Zgornji izraz vsebuje prvo lastnost, ki jo zahtevamo od kovariančne funkcije. To pomeni, da so točke, ki so si blizu v vhodnem prostoru, močno korelirane in zato omogočajo močan odziv pri podobnih vrednostih ?-jev. Definicija bližine je izražena v obliki izrazov množice dolžinskih meril {n}. Za vsak vhod obstaja ustrezna dolžinsko merilo, ki označuje razdaljo v določeni smeri preko katere pričakujemo pomembne spremembe t. Zelo veliko dolžinsko merilo označuje, da imajo napovedi izdelane z uporabo Gaussovih procesov, zelo majhno ali celo ničelno odvisno razmerje z ustreznim vhodom. Za tak vhod lahko rečemo, da je nepomemben. Takšna interpretacija je v veliki povezavi z avtomatskim določanjem značilk [37], [42]. Hiperparameter 0j določa skupno vertikalno skalo relativno glede na srednjo vrednost Gaussovega procesa v izhodnem prostoru. Hiperparameter <92 pa določa vertikalno negotovost in s tem označuje kako daleč lahko pričakujemo, da bo srednja vrednost funkcije odstopala od srednje vrednosti Gaussovega procesa. Pomen posameznih hiperparametrov je najbolje predstavljen v preprostem primeru modeliranja nelinearne funkcije v dodatku B. Drugi izraz v kovariančni funkciji 5mn03 predstavlja model šuma. Predpostavljamo, da je šum naključen. Iz tega sledi, da ne pričakujemo nobene korelacije med šumom in določenimi izhodi ter, da drugi izraz vpliva samo na diagonalne elemente kovariančne matrike. Varianco šuma podaja 03, za katero se predpostavlja, da je v primeru stacionarne kovariančne funkcije neodvisna od vhodov (modeli šuma, ki so odvisni tudi od vhodov bodo predstavljeni v nadaljevanju). Predstavljena kovariančna funkcija omogoča dobro temeljno obliko za delo. Obstajajo pa tudi druge uporabne oblike stacionarnih kovariančnih funcij. Za primer poskušajmo izdelati model z uporabo periodične funkcije, ki ima znano valovno dolžino A, v vhodni smeri /. Kovariančna funkcija, ki vsebuje takšno periodičnost, je predstavljena v enačbi 2 - 1K | " " ) (4.37) Cf{xm,xn) = 01e 2,=1 " 4.2.4.Stacionarne kovariančne funkcije 35 Funkcije, ki jih opisuje taka kovariančna funkcija nimajo močne korelacije samo med točkami, ki so si blizu, ampak tudi med točkami, ki so med seboj oddaljene za dolžinov A, v smeri l. 4.2.5.Nestacionarne kovariančne funkcije Medtem, ko se da veliko večino množic podatkov efektivno modelirati z uporabo stacionarnih kovariančnih funkcij, obstajajo določeni primeri v katerih se mora kovariančna funkcija izražati tudi v obliki nestacionarnosti. V nadaljevanju je podan opis treh posebno uporabnih nestacionarnih oblik. 4.2.5.1.Linearni izraz Pri daljšem pomiku stran od znanih podatkov, napovedi z uporabno kovariančne funkcije, kot je opisana z enačbo (4.36), težijo h konstatni vrednosti. V nasprotju s tem pa lahko predvidevamo, da obstaja neka linearna težnja v množici podatkov. Predstavljajmo si ploskev, ki jo opisuje enačba y(*)=2al x l+c . Če imata množica {al) in c obliko Gaussove verjetnostne porazdelitve s srednjo vrednostjo enako 0 in variancama ,(*)) predstavlja množico baznih funkcij in {/>,}€0 primerne koeficiente. j (Zß/I>MJ (4. 63{xm;0)=e'-1 40) Tak model šuma ni stacionaren, je pa pozitivno definiten, ker prispeva samo k diagonalnim elementom kovariančne matrike. 4.2.5.3.Prostorsko spremenljiva dolžinska merila V običajni stacionarni funkciji, ki je bila opisana predhodno, predpostavljamo, da je v posamezni smeri vhodnega vektorja dolžinsko merilo konstantno. Z lahkoto pa si predstavljamo primer slabega modela glede na podatke. Zaradi tega pa lahko določimo, da postanejo {n} funkcije x. Ne moremo pa preprosto nadomestiti parametrizirano obliko za n(x) v enačbi (4.36), ker s tem ne dobimo splošne pozitivno definitne kovariančne funkcije. Primer takšne kovariančne funkcije je podrobneje obdelan v [26] in ga opisuje enačba mXmtXm;B)=Sj[^^r^K^2M^ (4.41) 4.2.6.Določanje hiperparametrov kovariančne funkcije Za učne točke smo zgradili model. Sedaj pa moramo najti ustrezen način, da bomo lahko dosledno določili neznane vrednosti hiperparametrov 0 zgrajenega modela. V idealnem primeru bi želeli za tvorjenje napovedi izračunati integral preko vseh hiperparametrov, oziroma najti P{tN+1\xN+1,D, C{.))=jP{tN+1\xN+1,D,C{.),0)P{0\D,C{.))d0 , (4.42) kjer C(.) predstavlja kovariančno funkcijo. Za poljubno obliko kovariančne funkcije je analitičen pristop nemogoč. Lahko pa pridobimo približek končne (“a-posteriori”) verjetnostne porazdelitve s približevanjem povprečja napovedi preko vseh možnih vrednosti hiperparametrov z uporabo najbolj verjetnih vrednosti hiperparametrov (metoda največje podobnosti) ali z izvajanjem numerične integracije preko 0 z uporabo Monte Carlo metode. 4.2.6.1.Metoda največje podobnosti 37 4.2.6.1.Metoda največje podobnosti Metoda določanja najbolj verjetnih vrednosti hiperparametrov uporablja za določitev vrednosti integrala (4.42) približek P(tN+1\xN+l,D,C())^P(tN+1\xN+l,D,C(.),@Mp) , (4.43) ki temelji na najbolj verjetni množici hiperparametrov 0Mp. Približek temelji na predpostavki, da ima končna verjetnostna porazdelitev 0, P(0\D,C(.)) obliko ostre konice okrog 0mp- V splošnem je opisani približek zelo dober in tako izdelane napovedi so pogosto zelo blizu pravim [37]. Najbolj verjetne vrednosti hiperparametrov pri določeni kovariančni funkciji določimo prek “a-posteriori” zapisa verjetnostne porazdelitve parametrov [36,47,26] Optimalne vrednosti parametrov določimo z iskanjem največje vrednosti logaritma porazdelitve, ki je logaritem Gaussovega procesa [26]: log{p{0\X,t)) = -2log{\CN\)-2tTNC-N1tN-2log{2TT) (4.45) Parcialni odvod logaritma porazdelitve po hiperparametrih se glasi: V literaturi se ta metoda imenuje metoda največje podobnosti (ang. maximum likelihood - ML). Za iskanje minimuma se lahko uporablja katerakoli optimizacij ska metoda. Ena izmed možnih je uporaba metode konjugiranih gradientov zaradi enostavnih analitičnih izračunov parcialnih odvodov (4.46). Predstavljena metoda je občutljiva na začetno izbiro hiperparametrov (padec v lokalni minimum), poleg tega pa je računsko zahtevna, saj vsak korak optimizacije potrebuje izračun inverzne kovariančne matrike dimenzije n x n, kjer je n število podatkov v učni množici. 4.2.6.2.Metoda z uporabo Monte Carlo pristopa Monte Carlo pristop uporablja metodo vzorčenja za izračun približka verjetnostne porazdelitve napovedi. Z aproksimacijo integrala v enačbi (4.42) z uporabo Markove verige (Markov chain), lahko zapišemo 38 4.2.6.2.Metoda z uporabo Monte Carlo pristopa P(tN+1\xN+l,D,C{))1YJP(tNAxN+l,D,C{.),0t) , (4.47) kjer predstavljajo 0t vzorce vzete iz končne (“a-posteriori) verjetnostne porazdelitve P(0|D,CQ) hiperparametrov 0 Vsak izraz v vsoti zgornje enačbe ima Gaussovo obliko, zaradi česar je tudi Monte Carlo približek napovedane verjetnostne porazdelitve mešanice Gaussovih funkcij. Natančnost približka narašča z večanjem uporabljenega števila vzorcev. Metoda, ki jo uporabljamo za vzorčenje funkcije končne verjetnostne porazdelitve hiperparametrov 0, ima bistven vpliv na učinkovitost pristopa. Želimo si, da uporabljeni vzorci dobro predstavljajo verjetnostno porazdelitev iz katere so bili vzeti. Če pri procesu vzorčenja izpustimo določeno območje 0 prostora, ki je močno povezano z verjetnostno porazdelitvijo, bo približek integrala enačbe (4.42) zelo slab. 4.2.6.3.Primerjava opisanih metod V dodatku B je prikazan preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi z uporabo metode največje podobnosti. Pri optimizaciji hiperparametrov včasih pride do primera, ko glede na različne začetne vrednosti hiperparametrov, dobimo različne končne optimirane vrednosti hiperparamtrov. Glede na izkušnje avtorja iz [26], ta problem ne predstavlja bistvene ovire za uporabo omenjene metode pri modeliranju z Gaussovimi procesi. Omenjeni problem lahko rešimo z uporabo lastnih predpostavk o hiperparametrih. Pridobivanje najbolj verjetnih hiperparametrov tudi predstavlja najbolj privlačno odliko metode največje podobnosti. Pri interpretaciji modela z Gaussovimi procesi to pomeni, da nam takšna množica hiperparametrov lahko veliko pove o podatkih, ki jih modeliramo. V primeru samo ene končne množice hiperparametrov to pomeni, da je za izračunavanje napovedi potrebno inverzno vrednosti kovariančne matrike izračunati in shraniti samo enkrat. Pri uporabi Monte Carlo pristopa se izračuna približek končne verjetnostne porazdelitve z uporabo izračuna povprečja iz množice vzorcev. Po tem, ko pridobimo vse informacije, ki jih potrebujemo za začetek Monte Carlo simulacije, lahko začnemo izračunavati povprečja korak za korakom. Pri tem moramo vsakič izračunati inverz kovariančne matrike, ni pa inverza potrebno shranjevati. Če pa se želimo vrniti na začetek in izračunati novo statistiko o podatkih ali dodati nekaj novih učnih podatkov, je potreben ponoven zahteven izračun vseh inverznih 4.2.6.3.Primerjava opisanih metod 39 kovariančnih matrik. Iz tega je razdvidno, da nam Monte Carlo pristop ne omogoča bistveno boljše prilagodljivosti kot metoda največje podobnosti. V Gaussovih procesih, ki uporabljajo Monte Carlo pristop, je možna tudi uporaba modelov šuma, ki nimajo Gaussove porazdelitve. Metoda uporabe Monte Carlo pristopa pa se izkaže tudi za bolj učinkovito v primeru uporabe kompleksnejših kovariančnih funkcij (konkretno v primeru uporabe modela šuma, ki je odvisen od vhodov). Uporaba metode največje podobnosti v nadaljevanju pa na noben način ne pomeni, da je metoda boljša kot Monte Carlo pristop za vse probleme. Za manjše [26] množice podatkov, kjer shranjevanje matrik ne predstavlja bistvenega problema, se je izkazalo, da daje Monte Carlo pristop boljše rezultate pri nespremenljivi velikosti CPU časa. Za večje množice podatkov pa je metoda največje podobnosti hitrejša in učinkovitejša. 4.3.Mehka logika Teorija mehke logike se je začela razvijati v šestdesetih letih, ko je Lofti A. Zadeh v članku [66] prvič obravnaval tako imenovane mehke ali zamegljene veličine, ki jih ni možno opisati s formalnim matematičnim zapisom, vključno s teorijo verjetnosti. Potrebo za njihovo uvedbo, kot tudi za razvoj ustreznih matematičnih postopkov, je utemeljil z dejstvom, da je človeško zaznavanje in opisovanje bioloških in mehanskih sistemov daleč od natančno definiranega matematičnega izrazoslovja in brez natančno poznanega zakona porazdelitve. Svoje delo je nadaljeval [65] z uvedbo pojma lingvistične spremenljivke, to je spremenljivke katere vrednost ni določena s številčnimi vrednostmi, pač pa z besedami iz naravnega jezika. Ta članek predstavlja delo, ki je bilo osnova nadaljnjemu teoretičnemu razvoju teorije mehke logike in njeni aplikaciji pri reševanju številnih problemov. Danes se je mehka logika uveljavila na mnogih področjih kot so razpoznavanje vzorcev, procesiranje signalov, elektronika, zlasti pa se je uveljavila na področju vodenja in identifikacije sistemov. Metode mehke identifikacije sta prva predlagala Takagi in Sugeno [52], podrobnejša predstavitev mehkih indentifikacijskih metod pa je na voljo v literaturi [49], [51], [54], kjer se tudi natančno podane osnovne definicije in metode, potrebne za razumevanje in izgradnjo sistema mehke logike. Predstavljen je tudi izrek, ki uvršča sisteme mehke logike med univerzalne aproksimatorje in s tem postavlja teoretično osnovo za njihovo aplikacijo pri aproksimaciji neznanih funkcij. 40 4.3.1.Mehka identifikacija sistemov 4.3.1.Mehka identifikacija sistemov Pri mehki identifikaciji sistemov lahko podobno kot pri klasičnih pristopih identifikacije problem razdelimo v dva glavna dela: določitev strukture in določitev parametrov sistema mehke logike. Določitev strukture sistema mehke logike poteka po sledečem postopku: • določitev vhodov, • določitev izhodov, • določitev tipa pravil (Mamdani, Takagi-Sugeno), • določitev načina sklepanja in mehkih operatorjev, • in določitev metode mehčanja vhodov in ostrenja izhodov. Mehke identifikacijske modele lahko razvrstimo v dve skupini glede na to, kakšne tipe pravil uporabljajo. Mamdani identifikacijski modeli so sestavljeni iz pravil oblike (4.48) kjer so vhodi in izhodi mehke vrednosti, torej je uporabljeno znanje v popolnosti lingvistično interpretirano. Adaptacijskih postopkov, ki bi vodili do dobrega modela, za takšne sisteme skoraj ni na voljo. Izjema predstavljajo genetski algoritmi, ki v aplikacijah v realnem času zaradi visoke računske zahtevnosti niso preveč uporabni. Rm:ifx1AAm1Ax2jeA2A...AxnisAmntheny = Ym (4.48) Takagi-Sugeno identifikacijski modeli imajo pravila oblike (4.49). Vsako pravilo predstavlja lineariziran lokalni model sistema, identifikacijski model pa izvaja interpolacijo med njimi. Izhod Takagi-Sugeno mehkega sistema pri uporabljeni težiščni metodi mehčanja je reduciran na uteženo povprečje. Rm:ifx1isAm1Ax2isAm2A...AxnisAmntheny = f{x1x2...,xn) (4.49) Določitev parametrov sistema mehke logike poteka po sledečem postopku: • določitev števila, oblike in položaja mehkih množic • in adaptacija vseh ali samo nekaterih parametrov. Diagram poteka mehkega sklepanja v mehkih identifikacijskih modelih je prikazan s pomočjo preprostega praktičnega primera na sliki 17 [54]. Informacije v diagramu potekajo vzporedno od treh vhodov do enega izhoda. Vzporedna obravnava pravil predstavlja tudi eno izmed pomembnejših lastnosti mehkih sistemov. Proces mehkega sklepanja je 4.3.1.Mehka identifikacija sistemov 41 sestavljen iz petih postopkov. 1. Mehčanje vhodov je postopek pri katerem se ugotovi stopnja pripadnosti posameznega vhoda posamezni mehki množici s pomočjo pripadnostnih funkcij. 2. Mehki operatorji v pogojnem delu (premisi) pravil se uporabijo zatem, ko se zaključi mehčanje vhodov, v primeru, da je pogojni del (premisa) sestavljen iz večih delov. 3. Uporaba implikacijske metode sledi zatem, ko je s pogojnim delom določen tudi posledični del (konsekvenca). V praksi se uporabljajo različne vrste implikacijskih metod [54]: Mamdani implikacija, Lukasiewicz implikacija in produkt. 4. Agregacija vseh izhodov je proces pri katerem se mehke množice, ki jih določajo konsekvence pravil z uporabo implikacijske metode, združijo v združeno mehko množico. Vhod v agregacijski proces predstavlja skupina vseh funkcij iz posledičnih delov vseh pravil, ki so bile dodatno okrnjene z uporabo implikacijske metode. Izhod agregacijskega procesa pa je po ena mehka množica za vsak izhod mehkega modela. 5. Ostrenje, ki sledi agregaciji, je proces, ki priredi vhodni mehki množici realno število. Vhod v proces ostrenja je združena mehka množica, izhod pa realno število, ki tudi predstavlja izhodno vrednost mehkega modela. 42 4.3.1.Mehka identifikacija sistemov Slika 17: Diagram poteka mehkega sklepanja v mehkih identifikacijskih modelih s pomočjo praktičnega primera “Dinner for two” iz programskega orodja Matlab [54] 4.3.2.Grupiranje mehkih množic Mehko zbiranje podatkov v gruče predstavlja učinkovit pristop k identifikaciji kompleksnih nelinearnih sistemov s pomočjo mehke logike. Mnogi razvrščevalni postopki in algoritmi za modeliranje sistemov uporabljajo kot temeljeno metodo zbiranje numeričnih podatkov v gruče (angl. clustering). Glavni smisel zbiranja podatkov v gruče je določitev števila naravnih gruč v veliki množici podatkov, s katerimi je možno jedrnato predstaviti značilnosti sistema. Informacije o gručah omogočajo izgradnjo Takagi-Sugeno mehkega sistema sklepanja z uporabo relativno majhnega števila pravil. To v bistvu pomeni, da se za vsako gručo podatkov določi svoje pravilo. Za zbiranje podatkov v gruče obstaja veliko metod, ki so podrobneje opisane v literaturi [49], [51], [54]. 4.3.3.ANFIS (Adaptive Neuro-Fuzzy Inference System) 43 4.3.3.ANFIS (Adaptive Neuro-Fuzzy Inference System) Med sistemi mehke logike in umetnimi nevronskimi mrežami obstaja podobnost, ki je najočitnejša pri strukturi. Takagi-Sugeno sistem mehke množice je na primer možno prikazati tudi v obliki tronivojske nevronske mreže, kar je podrobneje predstavljeno v [51]. Sistemi mehke logike so torej primerljivi po strukturi z umetnimi nevronskimi mrežami, možna pa je tudi uporaba enakih adaptacijskih algoritmov. Oba identifikacijska modela sta tudi univerzalna aproksimatorja za zvezne, nelinearne funkcije [51,29]. Za večnivojske nevronske mreže, ki imajo prav takšno strukturo, se običajno uporabljajo gradientno osnovani algoritmi za adaptacijo parametrov. Takagi-Sugeno mehke sisteme z gradientnim adaptacijskim postopkom imenujemo tudi ANFIS (adaptive nevro-fuzzy) sisteme. 4.4.Metode za vrednotenje modelov Za primerjavo rezultatov modelov zgrajenih z različnimi metodami potrebujemo primerno kriterijsko funkcijo s katero se lahko oceni uspešnost modela na preizkusnih podatkih. Kriterijska funkcija mora biti prilagojena področju uporabe. Oceniti mora predvsem tisto sposobnost modela, ki je poglavitni cilj razvoja [15], [40]. Pri modeliranju napovedovanja povišanih koncentracij ozona predstavlja poglavitni cilj zagotovitev čim boljšega napovedovanja visokih koncentracij. Nižje koncentracije za ocenjevanje niso zanimive, razen v primeru, če jih model napove s tako napako, da ta napoved povzroči tako imenovani lažni alarm. Za oceno uspešnosti modeliranja se na področju modeliranja onesnaženja ozračja uporablja precej cenilk. Vsaka od njih pa osvetli določen vidik uspešnosti ali neuspešnosti modela [15], [40]. V nadaljevanju bom predstavil cenilke, ki sem jih uporabil za ocenjevanje zgrajenih modelov. Pred tem pa podajam še opis pomena posameznih oznak: Cm, ... i-ta izmerjena vrednost koncentracije, Cr, ... i-ta izračunana vrednost koncentracije, C ... povprečna vrednost koncentracije, vc ... standarna deviacija koncentracije, N ... skupno število preizkusnih vzorcev, / ... skupno število napovedanih prekoračitev, 44 4.4.Metode za vrednotenje modelov m ... skupno število izmerjenih prekoračitev, a ... število pravilno napovedanih prekoračitev. 4.4. 1.Korelacijski koeficient Korelacijski koeficient predstavlja normalizirano mero za ugotavljanje linearne odvisnosti dveh spremenljivk. Za povsem linearno neodvisni spremenljivki znaša vrednost koeficienta 0.0, medtem ko za linearno odvisni znaša 1.0 ali -1.0. Izračun vrednosti korelacijskega faktor med dvema spremenljivkama X in Y predstavlja enačba (4.50). r(X,Y)= cov^X,Y) ; cov(X,Y) = E((X-E(X))(Y-E(Y))) (4 50) °X°'Y . Izračun korelacijskega koeficienta za vrednotenje modelov glede na enačbo (4.50) opisuje enačba (4.51). 1 N N^(Cm-Cm)iCr-Cr) r(Cm,Cr)=^^---------------------------- (4.51). O'Cm & Cr 4.4.2.Kvadratni koren povprečne kvadratne napake (RMSE) Povprečno kvadratno napako (RMSE..Root Mean Square Error) med izračunanimi vrednostmi koncentracij iz modela in izmerjenimi vrednostmi določa enačba enačba (4.52). RMSE=1 over N-Y, [Cr-Cm)2 (4.52). i=1 4.4.3.Cenilka p6 S stališča uporabe rezultatov napovedovanja v praksi je pomembno predvsem, da zna model dobro napovedati visoke koncentracije in da se pri nizkih koncentracijah ne zmoti toliko, da bi lahko sprožil lažni alarm. Ker nobena od doslej navedenih cenilk ne pokaže zadosti uspešnosti oziroma neuspešnosti napovedovanja prekoračenih koncentracij, sem za vrednotenje uporabil tudi cenilko p6 [40], ki sem jo prilagodil velikostnemu razredu koncentracij ozona. Cenilka p6 pomeni verjetnost uspešne napovedi visokih koncentracij. Postopek izračuna je sledeč: 1. Določitev vzorcev, ki bodo upoštevani v izračunu. Za izračun se uporabijo vsi vzorci pri katerih je izmerjena visoka koncentracija in 4.4.3.Cenilka p6 45 vsi vzorci pri katerih je nastopil lažni alarm. Pogoja opisuje enačba (4.53). Vzorcev, ki ne ustrezajo danemu pogoju, se ne upošteva. Cm^120vglm3vCr^180vglm3 (4.53) 2. Določitev uspešnosti napovedi J6 za vsak vzorec iz izbrane množice posebej. Posamezna napoved J6i je uspešna (enaka 1), če je absolutna napaka napovedi manjša od 20 µg/m3 ali če je relativna napaka manjša od 0.2 in pri tem ni nastopil lažni alarm. Opisani pogoji so predstavljeni z enačbo (4.54). 1, če velja:{{\Cm -Cr\<20[xglm)v{\{Cm -Cr^lCmi\<02)) A J6=- ((Cmi>150/ig/m3)v(Cri>180/ig/m3)) 0, v ostalih primerih (4.54) 3. Določitev cenilke verjetnosti uspešne napovedi visokih koncentracij p6 tako, da preštejemo vse uspešne napovedi (napovedi z oceno 1) in vsoto delimo s številom vseh izbranih vzorcev, kot to opisuje enačba (4.55). p6= 1^J6 (4.55) i=1 4.4.4.Parametri sposobnosti Za vrednotenje in ocenjevanje sposobnosti modelov, ki napovedujejo prekoračene koncentracije ozona, je EEA (European Environment Agency – Evropska agencija za okolje) z dokumentom [61] določila standardno tabelo deležev (contingency table), ki je prikazana v tabeli 1. Pomen oznak v tabeli je podan na začetku podpoglavja 4.4. Alarmi Izmerjeni Napovedani Da Ne Skupaj Da a f-a f Ne m-a N+a-m-f N-f Skupaj m N-m N Tabela 4.1: Tabela deležev (contingency table) Z uporabo navedenih definicij so definirani trije parametri sposobnosti. • Merilo pravilnih napovedi SP kritičnih dogodkov (verjetnost detekcije) določa enačba (4.56), ki predstavlja razmerje med številom pravilno napovedanih prekoračitev in skupnim številom izmerjenih 46 4.4.4.Parametri sposobnosti prekoračitev v odstotkih. Obseg merila je od 0% do 100%, kjer 100% označuje najboljšo vrednost. a SP=—-100% (4.56) m • Merilo uresničenih napovedi SR kritičnih dogodkov določa enačba (4.57), ki predstavlja razmerje med številom pravilno napovedanih prekoračitev in številom vseh napovedanih prekoračitev v odstotkih. Obseg merila je od 0% do 100%, kjer 100% označuje najboljšo vrednost. SR= f-100% (4.57) • Indeks uspeha, ki združuje parametra sposobnosti SP in SR, določa enačba (4.58) Obseg merila je od -100% do 100%, kjer 100% označuje najboljšo vrednost. SI=fma+N+Na:^-1)-100% (4.58) Za vrednotenje modelov sem izdelal dva izračuna parametrov sposobnosti. V prvem primeru sem nastavil kritično mejo koncentracije na 180 µg/m3, kot je to predpisano v uredbi. Parametre sposobnosti sem označil z SP180, SR180 in SI180. Ker pa je bilo število vzorcev v množici za preizkušanje, ki so presegli predpisano kritično vrednost, relativno majhno, sem se odločil v namene primerjave modelov napraviti še izračun pri kritični meji 140 µg/m3. V tem primeru sem parametre sposobnosti označil z SP140, SR140 in SI140. 47 5.Analiza izmerjenih podatkov 5.1.Geografski položaj Nove Gorice 5.1.1.Vipavska dolina in Goriško polje Območje Vipavske doline in Goriškega polja leži na meji med submediteransko Slovenijo in dinarskim delom Slovenije [33,38,39], kot je prikazano na sliki 18. Vipavska dolina in Goriško polje ležita med nizko planoto Krasa in visokimi kraškimi planotami Nanosa in Trnovskega gozda. Pri Šempasu pride Vipavska dolina v Goriško ravan, zahodno od obeh Goric pa se razprostira Furlanska nižina. Relief na področju Vipavske doline in Goriškega polja prikazuje slika 19. Slika 18: Izsek iz državne pregledne karte Slovenije prikazuje širšo okolico Vipavske doline in Goriškega polja Območje Vipavske doline in Goriškega polja je nižje flišno ozemlje, katerega približno tri četrtine zavzema gričevje. Preostali del ozemlja, to je slaba četrtina, pa je ravnina. Te je več pod Trnovskim gozdom na prodnih vršajih potokov in vzdolž reke Vipave. Največja ravnina je na stiku gričevnate regije s Soško ravnino. Soča je v ledeni dobi s prodom nasula na tej strani državne meje Solkansko in Šempetrsko polje. Zajezila je Vipavo, da je nižje od Prvačine odlagala ilovico - surovino mnogih nekdanjih in sedaj redkih opekarn. 15 48 5.1.1.Vipavska dolina in Goriško polje Med reko Vipavo in Krasom se Vipavsko gričevje dviga v hribovita Vipavska brda. Dno Vipavske doline ter gričevja na njenem jugozahodnem obrobju gradijo flišne plasti. Ta mehka kamnina, ki jo sestavljajo pole laporja, peščenjaka, apnene breče in glinastih skrilavcev, se je usedla na dnu terciarnega morja, ki je takrat še pokrivalo jugozahodni del Slovenije. V Vipavski dolini in na Goriškem polju prevladuje submediteransko podnebje. Ker loči Kras Vipavsko dolino od obale, so njene temperature nižje kot ob obali, posebno ponoči in pozimi. Zaradi goratega zaledja je količina letnih padavin relativno velika. Poletna suša je omiljena. V zahodnem delu je pod vplivom sredozemskega podnebja. Pozimi je gornji del doline pogosto izpostavljen burji, ki piha povprečno 42 dni v letu od tega 30 dni precej močno. Po Vipavski dolini teče reka Vipava ter po spodnjem robu Vipavske doline reka Branica. V ti dve reki se zlivajo številni manjši potoki. Večji potoki so na primer Hubelj, Skrivšek, Malenšček, Košivec, Konjščak. Za razliko od Krasa ima Vipavska dolina dve večji tekoči površinski reki. To ji omogočajo flišne plasti, ki so le delno prepustne za vodo. Slika 19: Relief na področju Vipavske doline in Goriškega polja. 5.1.2.Nova Gorica Sodobno urejeno mesto na robu Goriškega polja, tik ob slovensko-italijanski meji, je bilo ustanovljeno leta 1947, ko je z določitvijo meje med Italijo in Jugoslavijo Gorica ostala na italijanskem ozemlju, spodnja Soška dolina, Goriška brda in spodnja Vipavska dolina na jugoslovanski 5.1.2.Nova Gorica 49 strani meje pa so ostali brez naravnega središča [33,39,38]. Solkan, edino večje naselje spodnje Vipavske doline, je bilo tedaj kraj brez tradicionalnih središčnih funkcij. Nova Gorica je torej eno najmlajših slovenskih mest. Zaradi ugodne prometne lege na stičišču Vipavske in Soške doline, Krasa in italijanske Furlanske nižine se je hitro razvila v regionalno gospodarsko, kulturno, izobraževalno, upravno in prometno središče. v Mesto se razteza po prodnih nanosih Soče in glinastih usedlinah srednjem delu porečja hudournika Koren in kaže, da se bo zraslo z naselji Kromberk, Pristava, Rožna dolina, Šempeter pri Gorici in Vrtojba. Nadmorska višina mesta znaša 92 m, število prebivalcev pa je večje od 35.000 [39]. Zemljevid mesta prikazuje slika 20. Slika 20: Zemljevid Nove Gorice Skozi Novo Gorico teče promet s tujino. Železniška proga, tri ceste in sedem mejnih prehodov (pri Vrtojbi in v Rožni Dolini sta mednarodna) povezuje njeno širšo mestno območje s sosednjo Gorico v Italiji. Bližina meje in Gorice ji daje značaj obmejnega mesta. Podnebje je prehodno submediteransko z veliko sonca glede na večino ostalih krajev v Sloveniji, kot je to prikazano na grafu 22, ter povprečno milimi in vlažnimi zimami in vročimi poletji. Slednje nazorno prikazuje graf 21, kjer je prikazan letni potek mesečnih povprečnih temperatur 50 5.1.2.Nova Gorica zraka v večjih krajih v Sloveniji za obdobje od leta 1991 do leta 2000. Podatki na omenjenih zadnjih dveh grafih ne vsebujejo podatkov izmerjenih v Novi Gorici, ker so bili na voljo samo podatki meteorološke merilne postaje Bilje v bližini Nove Gorice. Nova Gorica je pridobila avtomatsko merilno postajo za spremljanje onesnaženja ozračja in meteorologije šele v letu 2001 v okviru projekta Phare, kot je bilo že omenjeno. 25 -, 22.5 20 17.5 O 15 2- 12.5 10 7.5 5 2.5 0 -2.5 -5 #v >L> ^ ^N .#> ^ ¦ Rateče ¦ Bilje pri Novi Gorici v Postojna a Ljubljana [> Novo Mesto < Celje m Maribor x Murska Sobota ¦ Portorož 3> Slika 21: Letni potek mesečnih povprečnih temperatur zraka v večjih krajih v Sloveniji za obdobje od leta 1991 do leta 2000 [3] 5.1.2.Nova Gorica 51 325 300 275 250 225 200 175 150 125- 100-1 75 50- 25- 0 ¦ Rateče ¦ Bilje pri Novi Gorici v Postojna a Ljubljana i> Novo Mesto •* Celje m Maribor x Murska Sobota ¦ Portorož / v* i> ^ ^ ^ & ^ ^ J oÜ>N ^> Slika 22: Letni potek mesečnega povprečnega trajanja sončnega obsevanja v urah v večjih krajih v Sloveniji za obdobje od leta 1991 do leta 2000 [3] 5.2.Meritve meteoroloških podatkov in koncentracij snovi v zunanjem zraku v Novi Gorici Za identifikacijo modela so bili izbrani meteorološki in imisijski podatki, ki se merijo in shranjujejo na avtomatski okoljski postaji v Novi Gorici. Izmerjene polurne vrednosti se prenašajo iz postaje na centralno enoto Agencije RS za okolje (ARSO) ter zatem na mestno občino Nova Gorica (MO Nova Gorica). Okoljski sistem MO Nova Gorica ima poleg sistema shranjevanja polurnih podatkov vgrajen tudi mehanizem za izračun in shranjevanje vrednosti, ki se uporabljajo za vrednotenje stanja okolja glede na uredbo o ukrepih za ohranjanje in izboljšanje kakovosti zunanjega zraka [57,58,59,60]. V uredbi se za vrednotenje posameznih okoljskih veličin uporabljajo urne, dnevne, letne, osemurne drseče in štiriindvajseturne drseče vrednosti. Za ozon je predpisana opozorilna urna vrednost koncentracije 180 µg/m3 ter alarmna vrednost 240 µg/m3. Avtomatska merilna postaja meri in shranjuje meteorološke veličine: • temperaturo zraka, • relativno vlažnost, • moč globalnega sončnega sevanja na enoto površine, • vrednost zračnega pritiska, 52 5.2.Meritve meteoroloških podatkov in koncentracij snovi v zunanjem zraku v Novi Gorici • hitrost in smer vetra. Poleg meteoroloških veličin se merijo in shranjujejo tudi podatki o koncentracija snovi v zunanjem zraku: • koncentracija delcev PM(10)3, • koncentracija žveplovega dioksida (SO2), • koncentracija dušikovega monoksida (NO), • koncentracija dušikovega dioksida (NO2), • koncentracija dušikovega oksida (NOx), • koncentracija ogljikovega monoksida (CO) • in koncentracija ozona (O3). Dodatno pa je v postaji vgrajen tudi analizator hlapnih organskih spojin (VOC), ki meri sledeče parametre: • koncentracijo benzena, • koncentracijo toluena, • koncentracija (orto)ksilena • in vsoto koncentracij hlapnih organskih spojin. Zaradi tehničnih težav z merilnikom podatki o izmerjenih vrednostih hlapnih organskih spojin (VOC) niso dostopni in jih na žalost tudi ni bilo mogoče uporabiti pri modeliranju. Tudi izmerjenih vrednosti koncentracije prahu v ozračju ni bilo mogoče uporabiti zaradi neenotne statistične obdelave izmerjenih vrednosti, ki se je v izbranem časovnem intervalu spreminjala. Za modeliranje sem uporabil izmerjene urne povprečne vrednosti vseh omenjenih parametrov, ki so bili na voljo od dne 01.01.2002 do dne 31.12.2004. Vse podane urne vrednosti sem moral še dodatno urediti tako, da sem izpustil tiste urne intervale za katere niso bile na voljo urne vrednosti vseh izbranih parametrov. Po končanem urejanju mi je ostalo na voljo od 26304 možnih urnih vrednosti 14520 uporabnih urnih vrednosti (55.2%). V nekaterih urnih intervalih so manjkale izmerjene vrednosti določenih parametrov zaradi različnih vzrokov kot so kalibracije merilnika vsakih 23 ur (ker se običajno merilniki ne 3 Delci PM(10) so delci v zraku, ki jih prepušča filter s 50% nepropustnostjo za delce z aeorodinamskim premerom 10 µm. 5.2.Meritve meteoroloških podatkov in koncentracij snovi v zunanjem zraku v Novi Gorici 53 kalibrirajo hkrati, je izpad večkraten) ter popravila senzorjev in merilnikov, ki so včasih trajala tudi po več dni. Vse uporabne urne povprečne vrednosti meteoroloških in imisijskih parametrov sem tudi dodatno analiziral. Pri analizi sem za vsak parameter izrisal histogram izmerjenih vrednosti in določil minimalno vrednost, maksimalno vrednost in povprečno izmerjeno vrednost. Na slikah od 23 do 26 so prikazani histogrami meteroloških parametrov, ki se ujemajo s klimatološkimi značilnostmi Nove Gorice. Za statistično predstavitev vetra je uporabljena roža vetrov, ki v obliki frekvenčnega diagrama prikazuje odstotek vetra iz posamezne smeri. Roža vetrov je prikazana na sliki 27, statistika, iz katere je bila roža vetrov sestavljena, pa je prikazana na sliki 28. Na slikah od 29 do 34 pa so prikazani histogrami koncentracij snovi v zunanjem zraku. Iz histogramov je razvidno, da so bile glede na veljavne uredbe [57,60,58,2] izmerjene relativno nizke koncentracije žveplovega dioksida (mejna urna vrednost po uredbi znaša 350 µg/m3) in ogljikovega monoksida (mejna 8-urna vrednost po uredbi znaša 10 mg/m3). Medtem pa so izmerjene vrednosti dušikovega monoksida, duškovih oksidov in dušikovega dioksida relativno visoke, čeprav so še pod dovoljenimi vrednostmi (mejna urna vrednost po uredbi znaša za dušikov dioksid 200 µg/m3). Pri analizi podatkov sem ugotovil, da so bile opozorilne urne koncentracije ozona kar 82 krat presežene. Pri tem je potrebno poudariti, da so visoke koncentracije nad opozorilno mejo nastopile v enem dnevu kar nekajkrat, kar pomeni, da je bilo v tem času glede na podane podatke kritičnih 19 dni. 54 5.2.Meritve meteoroloških podatkov in koncentracij snovi v zunanjem zraku v Novi Gorici 3000 2500 2000 1500 1000 500 0 Temperatura zraka [st. C] min=-7.7 ¦ O o O O in P LTI d d o *- CM o LO i i O O o P d LO d LO *~ *~ Odsek maks=37.9 1 ¦ o o O o LO O m o CN co ^- O o o o O o LO CN CN co i ni 3500 3000 2500 2000 1500 1000 500 Relativna vlaga hin=0.0 povp=70.0 maks=100.0 S 9 10 Odsek Slika 23: Histogram temperature Slika 24: Histogram relativne vlage zraka 12000 10000 > 8000 0) o o > 6000 g > 55 4000 2000 Globalno sončno sevanje [W/m2] min=0.0 povp=131.8 maks=962.0 IUlSllSli 2 3 4 Ui . tj . . tj . ^¦| I LO I I (O I 5 6 7 Odsek i co i _ o _ 9 10 8000 Slika 25: Histogram moči globalnega sončnega sevanja na enoto površine Zračni pritisk povp=1003.2| [h Pa] maks=1024.3 Slika 26: Histogram zračnega pritiska ROŽA VETROV Podatki od 01.01.2402 00:00:00 do vključno 31.12.2004 00:00:00 Postaja: Nova Gorica SMERI <0.3 0.3-1.0 1.0-2.0 2.0-3.0 3.0-5.0 >5.0 SKUPAJ -11.3- 11.3 0.52 0.33 0.18 0.10 0.04 1.18 11 3 - 33 8 0 61 0 96 0.44 0.31 0.01 2.34 N ROŽA VETROV 33.8 - 56.3 0.71 0.57 0.21 0.06 0.00 1.55 Postaja: Mova Gorica Od: 01.01.2002 00:00:00 56.3 - 78.8 1.63 0.38 78.8-101.3 4.92 0.54 0.04 0.11 0.01 0.03 0.00 2.06 0.00 5.60 Do: 31.12.2004 00:00:00 101.3-123.8 4.83 1.02 0.14 0.10 0.00 6.09 / /' Veter Smeri Suma Smeri Suma 0 1.187,, 8 2.427, 123.8-146.3 2.30 1.42 146.3-168.8 1.36 1.19 0.34 0.13 0.07 0.05 0.00 4.13 0.00 2.73 ^T^ -----I—II E *- 1 2.34% S 2.46% 2 1.55% 10 2.76% 3 2.06% 11 3.32% 4 5.60% 12 4.40% 5 6.09%, 13 5.157, 6 4.13% 14 4.55% 168.8-191.3 1.36 0.84 191.3-213.8 1.32 0.99 213.8-236.3 1.46 1.14 236.3-258.8 1.64 1.48 0.22 0.14 0.16 0.17 0.00 0.01 0.01 0.02 0.00 2.42 0.00 2.46 0.00 2.76 0.00 3.32 t^!'8XZ7 7 2.73%, 15 1.44% Brezveterje: 47.03%, 258.8-281.3 2.31 1.57 0.45 0.07 0.01 4.40 ,\U *%M / 281.3-303.8 3.86 1.04 0.22 0.03 0.00 5.15 Štev. vseh inter.= 14688 Št.sprej. int.= 14688 Št.dobrih inter.= 14588 303.8-326.3 4.06 0.44 326.3-348.8 1.07 0.28 0.05 0.07 0.00 0.02 0.00 4.55 0.00 1.44 / 8 J. M SKUPAJ 47.83 33.96 14.20 3.08 0.88 0.06 \ \ --. ^ / ¦ >5.0 m/s ¦ 3.0-5.0 m/s Štev. vseh inter.= 14688 Št.dobrih inter =14688 Št.sprej int.= 14688 ~š~ / V 2.0-3.0 m/s I 1.0-2.0 m/s J 0.3-1.0 m/s O <= 0.3 m/s Slika 27: Roža vetrov Slika 28: Statistika rože vetrov 5.2.Meritve meteoroloških podatkov in koncentracij snovi v zunanjem zraku v Novi Gorici 55 Slika 29: Histogram koncentracij žveplovega dioksida Slika 30: Histogram koncentracij ogljikovega monoksida Slika 31: Histogram koncentracij dušikovega monoksida Slika 32: Histogram koncentracij dušikovega dioksida Slika 33: Histogram koncentracij dušikovih oksidov Slika 34: Histogram koncentracij ozona 56 5.3.Izbira vhodnih spremenljivk modela 5.3.Izbira vhodnih spremenljivk modela Za izgradnjo modela za napovedovanje maksimalne urne koncentracije ozona naslednjega dne glede na zbrane meterološke podatke in podatke o koncentracijah snovi v zunanjem zraku do 19 ure v tekočem dnevu sem moral najprej določiti vhodne spremenljivke v model in izhod. Za izhod modela sem določil maksimalno urno koncentracijo ozona, ki je bila izmerjena v prihajajočem dnevu od 19:00 ure v tekočem dnevu do 19:00 ure v prihajajočem dnevu. Glede na fizikalno kemijsko naravo procesa, priporočila in rezultate iz literature, izračunane vrednosti korelacijskih koeficientov med posameznimi vhodnimi veličinami in izhodno veličino (slika 44), ter na veliko število poizkusov so bili izbrane vhodne spremenljivke modela: 1. temperatura zraka (24 urna povprečna vrednost), 2. moč globalnega sončnega sevanja na enoto površine (24 urna povprečna vrednost), 3. koncentracija dušikovega monoksida (24 urna povprečna vrednost), 4. koncentracija dušikovega dioksida (24 urna povprečna vrednost), 5. koncentracija ozona (24 urna povprečna vrednost). 24 urna povprečna vrednost pomeni izračunano povprečje iz urnih povprečnih izmerjenih vrednosti od 19:00 ure prejšnjega dne do 19:00 ure današnjega dne. Med procesom modeliranja se je izkazalo, da je možno model izboljšati z uvedbo dodatnih vhodnih spremenljivk, ki vsebujejo napovedane vrednosti parametra za dan vnaprej. Za potrebe modeliranja so bile pridobljene iz izmerjenih vrednosti prihajajočega dneva. V praktični uporabi modela pa je predvideno pridobivanje teh vrednosti iz modela ALADIN Agencije RS za okolje. Vrednosti modela so objavljene na internetnih straneh ARSO [5]. Dodatne vhodne spremenljivke so: 6. napovedana maksimalna urna vrednost temperature zraka v prihajajočem dnevu, 7. napovedana komponenta hitrosti vetra v smeri sever-jug, pozitivna vrednost pomeni komponento vetra iz severne smeri, negativna pa komponento vetra iz južne smeri (povprečna vrednost izračunana iz vektorske hitrosti vetra za čas od 19:00 ure današnjega dne do 19:00 ure naslednjega dne), 5.3.Izbira vhodnih spremenljivk modela 57 8. napovedana komponenta hitrosti vetra v smeri vzhod-zahod, pozitivna vrednost pomeni komponento vetra iz vzhodne smeri, negativna pa komponento vetra iz zahodne smeri (povprečna vrednost izračunana iz vektorske hitrosti vetra za čas od 19:00 ure današnjega dne do 19:00 ure naslednjega dne). Vrednosti korelacijskih koeficientov med posameznimi vhodi in izhodom so prikazane na sliki 44 in so predstavljene v enakem vrstnem redu kot na predhodnem spisku izbranih vhodov. Pomen korelacijskega koeficienta in postopek izračuna je predstavljen v podpoglavju 6.1.1. Množico vseh vzorcev sem za potrebe modeliranja razdelil v dve množici. • Množico za učenje modelov, ki vsebuje vzorce iz štirih časovnih intervalov: 01.01.2002-31.07.2003, 01.09.2003-31.12.2003, 01.02.2004-31.08.2004, 01.10.2004-31.12.2004. Po izločevanju neuporabnih vzorcev zaradi različnih napak, ki so bile posledice različnih vzrokov (kalibracije merilnikov, okvare merilnikov in senzorjev, neuspešne komunikacije), je bilo možno za identifikacijo uporabiti 488 dobrih vzorcev izmed 1004 največ možnih vzorcev. • Množico za vrednotenje modelov, ki vsebuje vzorce iz treh časovnih intervalov: 01.08.2003-31.08.2003, 01.01.2004-31.01.2004, 01.9.2004-30.09.2004. Po izločevanju neuporabnih vzorcev zaradi različnih napak, je bilo možno za preizkušanje uporabiti 68 dobrih vzorcev izmed 92 največ možnih vzorcev. Potek posameznih vhodnih spremenljivk je prikazan na slikah od 35 do 42. Na sliki 43 pa je prikazan potek izhodne spremenljivke. Rdeča (svetlejša) krivulja označuje množico za učenje, modra (temnejša) krivulja pa množico za preizkušanje modelov. Horizontalna os grafov Dan označuje številko dneva šteto od dne 01.01.2002. To pomeni, da je na primer z indeksom 0 označen dan 01.01.2002, z indeksom 1096 pa 31.12.2004. 58 5.3.Izbira vhodnih spremenljivk modela Slika 35: Vhod 1, Temperatura zraka Slika 36: Vhod 2, Moč globalnega sončnega sevanja na enoto površine Slika 37: Vhod 3, Koncentracija dušikovega monoksida (NO) 5.3.Izbira vhodnih spremenljivk modela 59 Slika 38: Vhod 4, Koncentracija dušikovega dioksida (NO2) Slika 39: Vhod 5, Koncentracija ozona (O3) Slika 40: Vhod 6, Napovedana maksimalna temperatura zraka 60 5.3.Izbira vhodnih spremenljivk modela Slika 41: Vhod 7, Napovedan veter v smeri sever-jug Slika 42: Vhod 8, Napoved vetra v smeri vzhod-zahod Slika 43: Izhod, Maksimalne dnevne koncentracije ozona 19 5.3.Izbira vhodnih spremenljivk modela 61 Slika 44: Vrednosti korelacijskih koeficientov med posameznimi vhodnimi spremenljivkami v model in izhodno spremenljivko modela 63 6.Sinteza modelov za napovedovanje povišanih koncentracij ozona Za napovedovanje povišanih koncentracij ozona sem glede na izbrane vhodne in izhodne spremenljivke zgradil tri različne modele z uporabo treh različnih metod modeliranja opisanih predhodno: z umetno nevronsko mrežo, Gaussovimi procesi in mehko logiko. Za izgradnjo vseh treh modelov je bila uporabljena ista učna množica, ki je bila predstavljena v prejšnjem poglavju. 6.1.Sinteza modela z umetno nevronsko mrežo Za model je bila izbrana večnivojska perceptronska (multi-layer feedforward) umetna nevronska mreža, ki je bila zgrajena s pomočjo orodja za uporabo umetnih nevronskih mrež v programskem paketu Matlab poimenovanega Neural Network Toolbox [54] (funkcija newff()). Model vsebuje dva nivoja, nevroni pa uporabljajo v obeh primerih tangens-sigmoidno prenosno funkcijo. Pri procesu modeliranja sem po velikem številu poizkusov določil, da je najbolj primerno število skritih nevronov 16. Začetne vrednosti uteži in odmikov so bile nastavljene z uporabo funkcije initnw(). Učenje nevronske mreže je potekalo s pomočjo dveh množic podatkov: množice za učenje, ki je vsebovala 437 vzorcev, ter množice za sprotno preverjanje podatkov, ki je vsebovala 51 vzorcev. Za učenje je bila izbrana metoda Levenberg-Marquardt. Vsi izmerjeni podatki v model so bili predhodno normalizirani. Po končanem učenju je bil najprej preizkušen model na podatkih iz učne množice. Odvisnost med izmerjenimi in izračunanimi vrednostmi je prikazana na sliki 45. Na sliki 46 pa je prikazan še histogram napak, ki kaže značilnosti normalne porazdelitve. Napake so izračunane kot razlika med izmerjeno in modelirano vrednostjo. Iz podatkov, ki so bili uporabljeni za indentifikacijo modela, sem izračunal tudi vrednosti cenilk, ki sem jih uporabljal za vrednotenje modela (poglavje 4.4): r=0.949, RMSE=14.8, p6=0.93, SP180=40.0, SR180=85.7, SI180=39.8, SP140=86.1, SR140=76.7, SI140=82.1. Poudarim pa naj, da prikazane vrednosti ne odražajo prave kvalitete modela (v smislu zmožnosti posploševanja), ampak so samo informativne narave in prikazujejo kako dobro se je model naučil učne podatke. 64 6.1.Sinteza modela z umetno nevronsko mrežo Slika 45: Rezultati preizkusa modela z umetno nevronsko mrežo s podatki iz učne množice Slika 46: Histogram napak preizkusa modela z umetno nevronsko mrežo s podatki iz učne množice 6.2.Sinteza modela z Gaussovimi procesi 65 6.2.Sinteza modela z Gaussovimi procesi Model z Gaussovimi procesi je bil zgrajen s pomočjo programskega orodja v okolju Matlab, ki je na voljo na domači spletni strani dr. C. E. Rasmussen-a [48]. Programsko orodje predstavljata dve Matlab funkciji. • gp01lik() funkcija, ki omogoča optimizacijo hiperparametrov modela z Gaussovimi procesi, ki uporablja stacionarno kovariančno funkcijo. Optimizacija poteka po metodi največje podobnosti. Vhod v funkcijo predstavljajo začetne vrednosti hiperparametrov in množica učnih točk, ki je sestavljena iz parov {izmerjeni vhodi, izmerjen izhod}. Izhod funkcije pa so optimirane vrednosti hiperparametrov. • gp01pred() funkcija, ki omogoča na podlagi podanih hiperparametrov in učnih točk napoved izhoda modela pri podanem vhodu. Vhod v funkcijo predstavljajo optimalne vrednosti hiperparametrov in množica učnih točk, ki je sestavljena iz parov {izmerjeni vhodi, izmerjen izhod}, ter množica vhodov, za katere želimo napovedati izhodne vrednosti. Izhod pa predstavlja množica napovedanih izhodnih vrednosti, ki jo sestavljajo pari {napovedana srednja vrednost izhoda, napovedana dvojna vrednost standardne deviacije izhoda}. Glede na število vhodov, ki je znašalo v konkretnem primeru 8, sem določil 8 hiperparametrov GP modela (hiperparametri ?1...?8 predstavljajo dolžinska merila vsake dimenzije (vhoda), ?9 določa skupen vpliv vseh učnih točk na napovedane vrednosti, ?10 določa minimalno negotovost napovedanih točk). Začetne vrednosti vseh hiperparametrov sem pred optimizacijo nastavil na konstantno vrednost 0.5. Po optimizaciji so bile vrednosti hiperparametrov sledeče: 01 02 03 04 05 06 07 08 09 01O 0,270 0,061 0,019 0,150 0,342 0,885 2,941 1,772 0,256 0,022 Tabela 6.1: Vrednosti hiperparametrov Po končanem učenju je bil najprej preizkušen model na podatkih iz učne množice. Odvisnost med izmerjenimi in izračunanimi vrednosti je prikazana na sliki 47. Na sliki 48 pa je prikazan še histogram napak, ki kaže značilnosti normalne porazdelitve. Napake so bile izračunane kot razlika med izmerjeno in modelirano vrednostjo. Iz podatkov, ki so bili uporabljeni za indentifikacijo modela sem izračunal tudi vrednosti cenilk, ki sem jih uporabljal za vrednotenje modela (poglavje 4.4): r=0.947, RMSE=15.1, p6=0.93, SP180=26.7, SR180=100.0, SI180=26.7, SP140=84.6, SR140=79.7, SI140=81.3. 66 6.2.Sinteza modela z Gaussovimi procesi Zopet poudarjam, da prikazane vrednosti ne odražajo prave kvalitete modela, ampak so samo informativne narave. Slika 47: Rezultati preizkusa modela z Gaussovimi procesi s podatki iz učne množice Slika 48: Histogram napak preizkusa modela z Gaussovimi procesi s podatki iz učne množice 6.3.Sinteza modela z mehko logiko 67 6.3.Sinteza modela z mehko logiko Mehki model je bil zgrajen s pomočjo orodja FMID (Fuzzy Modeling and Identification Toolbox) v programskem okolju Matlab, katerega avtor je prof. dr. Robert Babuška [10]. Orodje sestavlja skupina Matlab funkcij za izgradjo Takagi-Sugeno (TS) mehkih modelov iz množice izmerjenih podatkov. Model temelji na grupiranju izmerjenih podatkov z mehko kovariančno matriko s pomočjo Gustafson-Kessel algoritma [51, 54]. Algoritem predstavlja razširitev običajnega algoritma mehke metode c-mehkih srednjih vrednosti (fuzzy c-means) z uporabo prilagodljive dolžinske norme. S tem je možno v množici vseh izmerjenih podatkov detektirati podmnožnice podatkov različnih geometrijskih oblik. Za izgradjo modela je bila bila uporabljena funkcija fmclust(). Končne vrednosti vhodnih parametrov funkcije (podrobnejši opis je podan v [10]) so bile določene glede na izkušnje pridobljene med procesom modeliranja: • število grup (določa tudi število pravil): 3, • parameter mehkosti (večja vrednost določa mehkejše (bolj prekrivajoče se grupe): 2, • tip pripadnostnih funkcij pogojnega dela mehkega modela: 2 (uporaba projeciranih pripadnostnih funkcij, ki naj bi glede na priložena navodila orodja omogočala natančnejši model), • tip posledičnega dela mehkega modela: 1 (uporaba globalne metode najmanjših kvadratov). Parametri modela so bili dodatno izboljšani z uporabo funkcije anfis() [54], ki se poslužuje metode povratnega učenja. Za proces učenja je bila učna množica podatkov razdeljena na množico za učenje (437 vzorcev) in množico za sprotno preverjanje učenja (51 vzorcev). Po končanem učenju je bil najprej preizkušen model na podatkih iz učne množice. Odvisnost med izmerjenimi in izračunanimi vrednosti je prikazana na sliki 49. Na sliki 50 pa je prikazan še histogram napak, ki kaže značilnosti normalne porazdelitve. Napake so bile izračunane kot razlika med izmerjeno in modelirano vrednostjo. Iz podatkov, ki so bili uporabljeni za indentifikacijo modela sem izračunal tudi vrednosti cenilk, ki sem jih uporabljal za vrednotenje modela (poglavje 4.4): r=0.941, RMSE=16.0, p6=0.92, SP180=26.7, SR180=100.0, SI180=26.7, SP140=81.5, SR140=76.8, SI140=77.8. 68 6.3.Sinteza modela z mehko logiko Zopet poudarjam, da prikazane vrednosti ne odražajo prave kvalitete modela, ampak so samo informativne narave. Slika 49: Rezultati preizkusa modela z mehko logiko s podatki iz učne množice Slika 50: Histogram napak preizkusa modela z mehko logiko s podatki iz učne množice 69 7.Vrednotenje modelov Po končanem procesu učenja sem vse tri zgrajene modele tudi preizkusil s pomočjo preizkusne množice, ki je vsebovala podatke za mesece: avgust 2003, januar 2004 in september 2004. Kot sem že predhodno omenil, ti podatki niso bili uporabljeni v procesu učenja modelov. Za vse tri modele sem pri preizkušanju izračunal tudi vse štiri vrste cenilk, ki so predstavljene v prejšnjem poglavju. Na slikah v nadaljevanju so prikazani rezultati preizkušanja modelov zgrajenih z umetno nevronsko mrežo (slika 51), Gaussovimi procesi (slika 54) in mehko logiko (slika 57) . Z modro (temnejšo) krivuljo in krogci so prikazane dejansko izmerjene vrednosti, z rdečo (svetlejšo) krivuljo in križci pa so prikazane izračunane (napovedane) vrednosti izhoda modela. Model z Gaussovimi procesi podaja tudi dvojno vrednost standardne deviacije za vsako točko napovedi, kar je na grafu prikazano s črtkano črto. Na slikah 55, 58 in 61 je za vsak model predstavljeno tudi odstopanje med izmerjenimi in izračunanimi vrednostmi. Na slikah 53, 56 in 59 so za vsak model prikazani še histogrami napak. 70 7.1.Preizkušanje modela na osnovi umetne nevronske mreže 7.1.Preizkušanje modela na osnovi umetne nevronske mreže Rezultati simulacije Slika 51: Rezultat preizkušanja modela na osnovi umetne nevronske mreže Slika 52: Odstopanje med izmerjenimi in napovedanimi vrednostmi Slika 53: Histogram napak preizkušanja modela na osnovi umetne nevronske mreže r=0.955, RMSE=17.5, p6=1.00, SPm=75.0, SR180=60.0, SIm=71.8, SP140=100.0, SR140=88.2, SI140=96.2 7.2.Preizkušanje modela na osnovi Gaussovih procesov 71 7.2.Preizkušanje modela na osnovi Gaussovih procesov Rezultati simulacije Slika 54: Rezultat preizkušanja modela na osnovi Gaussovih procesov Slika 55: Odstopanje med izmerjenimi in napovedanimi vrednostmi Slika 56: Histogram napak preizkušanja modela na osnovi Gaussovih procesov r=0.956, RMSE=16.2, p6=0.95, SPm=25.0, SR180=33.3, SIm=21.9, SP140=93.3, SR140=93.3, SI140=91.5 72 7.3.Preizkušanje modela na osnovi mehke logike 7.3.Preizkušanje modela na osnovi mehke logike Rezultati simulacije Slika 57: Rezultat preizkušanja modela na osnovi mehke logike Slika 58: Odstopanje med Slika 59: Histogram napak izmerjenimi in napovedanimi preizkušanja modela na osnovi vrednostmi mehke logike r=0.953, RMSE=18.1, p6=1.00, SP180=50.0, SR180=50.0, SI180=46.8, SP140=100.0, SR140=88.2, SI140=96.2 7.4.Primerjava rezultatov preizkušanja modelov 73 7.4.Primerjava rezultatov preizkušanja modelov 1 0,98 0,96 0,94 0,92 0,9 U UNM ¦ GP Ü ML 20 15 10 5 0 ? UNM ¦ GP D ML RMSE 1 0,98 0,96 0,94 0,92 0,9 ¦ UNM ? GP ? ML p6 Slika 60: Primerjava Slika 61: Primerjava Slika 62: Primerjava korelacijskega faktorja RMSE cenilke p6 80 70 60 50 40 30 20 10 0 ¦ UNM ¦ GP ? ML SP180 Slika 63: Primerjava merila pravilnih napovedi SP pri meji 180 mg/m3 100 95 90 85 80 ? UNM ¦ GP ? ML SP140 Slika 66: Primerjava merila pravilnih napovedi SP pri meji 140 mg/m3 70 60 50 40 30 20 10 0 ¦ UNM ¦ GP ? ML SR180 Slika 64: Primerjava merila uresničenih napovedi SR pri meji 180 mg/m3 100 95 90 85 80 ¦ UNM ¦ GP ? ML SR140 Slika 67: Primerjava merila uresničenih napovedi SR pri meji 140 mg/m3 ? UNM ¦ gp ? ml SI180 Slika 65: Primerjava indeksa uspeha SI pri meji 180 mg/m3 100 95 90 85 80 ¦ UNM ? gp ? ml SI140 Slika 68: Primerjava indeksa uspeha SI pri meji 140 mg/m3 r 74 7.4.Primerjava rezultatov preizkušanja modelov Primerjave posameznih cenilk so v splošnem pokazale, da mi je uspelo zgraditi tri skoraj enakovredne modele. Primerjava korelacijskega faktorja (slika 60) in kvadratnega korena povprečne kvadratne napake (slika 61) kaže, da je morda model z Gaussovimi procesi za odtenek boljši od ostalih dveh. Vendar pa je glede na glavni namen te naloge, ki je napovedovanje povišanih koncentracij ozona, vrednost cenilke p6, ki določa verjetnost uspešne napovedi visokih koncentracij, izkazalo da se je model na preizkusni množici odrezal malce slabše. Primerjava parametrov sposobnosti pri kritični meji koncentracije 180 mg/m3 je sicer pokazala, da se je v primeru ocenjavanja vseh treh parametrov sposobnosti SP, SR in SI, najbolje izkazal model na osnovi umetne nevronske mreže. Grafikoni rezulatov so prikazani na slikah 63, 64 in 65. Ker pa so omenejni rezultati zaradi majhnega števila vzorcev, ki so prekoračili kritično mejo (4 vzorci) vprašljivi, sem napravil tudi oceno parametrov sposobnosti pri nižji kritični meji koncentracije 140 mg/m3. Grafikoni rezultatov so prikazani na slikah 66, 67 in 68. V tem primeru se je zopet izkazalo, da je morda model na osnovi Gaussovih procesov za odtenek slabši od ostalih dveh, v glavnem pa se je izkazalo, da imamo opravka z enakovrednimi modeli. 75 8.Zaključek Prikazani rezultati primerjave vrednotenja modelov v prejšnjem poglavju kažejo, da je mogoče z različnimi pristopi modeliranja zgraditi povsem enakovredne modele. Pri tem je seveda potrebno upoštevati, da so bili vsi modeli zgrajeni na podlagi istih učnih podatkov in ovrednoteni na podlagi istih podatkov za preizkušanje. Postopek izgradnje modela je za vse tri postopke modeliranja po mojih izkušnjah precej podoben in zahteva približno enako mero predznanja. Kar je tudi razumljivo, saj so si metode sorodne in vse tri spadajo v isto skupino metod modeliranja. Sorodnost med metodami pa mi je pomagala, da sem po pridobitvi znanja in izkušenj ob proučevanju in učenju umetnih nevronskih mrež hitreje osvojil modeliranje z Gaussovimi procesi in mehko logiko. Pri učenju modelov z istimi učnimi podatki se je tudi izkazalo, da je bil za metodi modeliranja z umetno nevronsko mrežo in mehko logiko potreben približno enak čas učenja, ki je odvisen predvsem od računske zmogljivosti računalnika. Pri metodi učenja modela z Gaussovimi procesi pa je prišla računska zmogljivost računalnika še bolj do izraza, saj le ta pri optimizaciji hiperparametrov raste s tretjo potenco velikosti učne množice podatkov [53]. Za vse tri modele je bil potreben čas za izračun nove napovedi enak. Prednost modela z Gaussovimi procesi v primerjavi z modeloma z umetno nevronsko mrežo in mehko logiko vidim v tem, da je predikcija izhoda modela na osnovi Gaussovih procesov normalna verjetnostna porazdelitev s srednjo vrednostjo in varianco. Velikost variance, ki je odvisna od območja vhodnega prostora, podaja mero zaupanja v napovedano vrednost izhoda modela. Varianca je majhna na območjih, kjer je bila gostota točk učne množice velika, raste pa s približevanjem področjem, o katerih ima GP model relativno malo informacij. Druga prednost modela na osnovi Gaussovih procesov pa je optimizacija precej manjšega števila hiperparametrov, kot v primeru umetne nevronske mreže in mehke logike. Število hiperparametrov je vedno vnaprej določeno z obliko kovariančne funkcije in je običajno enako številu vhodov. 76 8.Zaključek Slabost, ki se je izkazala pri optimizaciji modela na osnovi Gaussovih procesov, pa je nastopala v primeru neustreznih začetnih vrednosti hiperparametrov. Uporabljena metoda največje podobnosti je v primeru slabo nastavljenih začetnih vrednosti hiperparametrov zelo hitro postala numerično nestabilna ali pa končala v lokalnem minimumu. Izkušnje pridobljene pri izdelavi naloge kažejo, da se je potrdila trditev iz [40,15], da je dobra identifikacija modela odvisna predvsem od pravilne izbire vhodnih spremenljivk. Glede na pregled priporočil iz literature se je izkazalo, da imajo nekatere vhodne spremenljivke bistveno manjši vpliv na izhod, kot je bilo omenjeno, ali pa celo onemogočajo dobro identifikacijo. Vhodni spremenljivki kot sta prah in žveplov dioksid v konkretnem primeru modeliranja za Novo Gorico pa celo poslabšujeta identifikacijo modela. Precej pa se je identifikacija izboljšala z uvedbo napovedi temperature in vetra. Med izvajanjem naloge se je izkazalo, da je zelo močna tudi povezava med koncentracijo ozona in temperaturo zraka, kar je tudi smiselno. Fotokemična reakcija nastajanja in razgradnje potrebuje določeno energijo, ki je v glavni meri sorazmerna temperaturi. Meje pri vrednotenju s cenilko p6 so bile postavljene razmeroma strogo zaradi primerjanja med različnimi modeli med procesom identifikacije. Dobljeni rezultati sicer povsem ustrezajo uredbi [4], ki predpisuje dovoljeno relativno napako napovedovanja ozona 50%, kljub temu pa rezultati kažejo, da bi se dalo model še malce izboljšati z večjo učno množico predvsem prekoračenih vrednosti, kar bo možno v prihodnosti. Obenem pa bo ob večji učni množici podatkov za preizkušanje tudi priložnost za sprotno ocenjevanje modela. Nova Gorica se bo kot urbano okolje v prihodnosti verjetno še širila, kar morda pomeni povečevanje števila prekoračenih vrednosti koncentracij ozona. Ker za vpliv na predvidene povišane koncentracije ozona v urbanem okolju zaenkrat še nimamo direktnih vzvodov (zaenkrat študije kažejo, da zmanjšanje prometa v mestih samo še poveča koncentracije ozona, ker se zmanjša koncentracija dušikov oksidov, ki porabljajo ozon), se mi zdi uporaba modela za napovedovanje povišanih koncentracij ozona edino smiselna in upravičena za varovanje zdravja ljudi. Zgrajene modele bi se dalo po mojem mnenju in glede na priporočila iz literature dodatno izboljšati z uporabo dodatnih vhodnih spremenljivk. Takšno skupino vhodnih spremenljivk, ki bi lahko glede na priporočila iz literature vplivala na izboljšavo predstavljajo koncentracije hlapnih organskih spojin (VOC). Merilna postaja v Novi Gorici je sicer imela ob postavitvi vgrajen merilnik hlapnih organskih spojih, vendar pa 8.Zaključek 77 zaradi tehničnih težav z merilnikom zaenkrat ni mogoče pridobiti uporabnih meritev. Pri tem pa bi bilo tudi potrebno raziskati vpliv vsake posamezne hlapne organske spojine na nastajanje ozona in njen pomen za potrebo izboljšanja napovedovanja povišanih koncentracij ozona. 79 9.Viri in literatura [1] Abdul-Wahab S. A., Al-Alawi S.M., Assessment and Prediction of TroposphericOzone Concentration Levels using Artificial Neural Netwoks, Environmental Modelling & Software, 17(3), 219-228, 2002. [2] Agencija Republike Slovenije za okolje, Mesečni bilten, http://www.arso.gov.si/o_agenciji/knji~znica/publikacije/bilten.htm, 12.11.2004 [3] Agencija Republike Slovenije za okolje, Podatki za nekatere postaje v obdobju 1991-2000, http://www.arso.gov.si/podro~cja/vreme_in_podnebje/napovedi_in_ podatki/podneb_10_tabele.html, 21.12.2004 [4] Agencija Republike Slovenije za okolje, Poročilo o stanju okolja v Sloveniji, http://www.arso.gov.si/poro~cila/Poro~cila_o_stanju_okolja_v_Slov eniji/zrak.pdf, 12.11.2004 [5] Agencija Republike Slovenije za okolje, Vreme in podnebje – napovedi in podatki, http://www.arso.gov.si/podro~cja/vreme_in_podnebje/napovedi_in_ podatki/index.html, 10.01.2004 [6] Ainslie B. , Steyn D. G., Scaling analysis of ozone precursor relationships, Air pollution modeling and its application XVI, Edited by Borrego and Incecik, Kluwer Academic/Plenum Publishers, New York, 321-328, 2004 [7] Aleksić N. M., Milutinović M., Milutinović P., Informacioni sistemi za kontrolu i upravljanje kvalitetom vazduha u urbanim sredninama, Institut za fiziku, Beograd [8] Atmospheric Research and Technology, LLC, About SODAR, http:// www.sodar.com/about_sodar.htm, 12.12.2004 [9] Ažman K., Identifikacija dinamičnih sistemov z Gaussovimi procesi z vključenimi linearnimi modeli, Magistrsko delo, Fakulteta za elektrotehniko, 2004. [10]Babuška R., Fuzzy Identification Toolbox for MATLAB 6.5 (ver. 3.03).http://www.dcsc.tudelft.nl/~babuska/, 28.02.2004 80 9.Viri in literatura [11]Ballester E. B., Valls G. C., Carrasco-Rodriguez J. L., Olivas E. S., S. del Valle-Tascon, Effective 1-day ahead prediction of hourly surface ozone concentrations in eastern Spain using linear models and neural networks, Ecological Modelling 156, 27-41, Elsevier, 2002 [12]Bochner S., Harmonic analysis and the theory of probability, Berkeley, 1979. [13]Box George E. P., Jenkins Gwilym M., Time series analysis forcasting and control, San Francisco, Holden-Day, 1976 [14]Božnar M. Z., Mlakar P., Use of neural networks in the field of air pollution modelling, Ath, Edited by Borrego and Schayes, Kluwer Academic/Plenum Publishers, New York, 2002 [15]Božnar M., Izbira učnih vzorcev za model napovedovanja ozračja na osnovi nevronskih mrež, Doktorska disertacija, Fakulteta za elektrotehniko, Univerza v Ljubljani, 1997 [16]Božnar M., Lesjak M., Mlakar P., A neural network-based method for short-term predictions of ambient SO2 concentrations in highly polluted industrial areas of complex terrain, Atmospheric Environment Vol. 27B, Pregamon Press Ltd, Great Britain, 221-230, 1993 [17]Comrie A. C., Comparing neural networks and regression models for ozone forcasting, Journal of Air & Wast Management Association 47, 653-663, 1997 [18]Dodge M., Chemical oxidant mechanisms for air quality modeling: critical review. Athmospheric Environment 34, 2103-2130, 2000 [19]EPA – Environmental Protection Agency, http://www.epa.gov/air/oaqps/eog/ozonehealth/aqi.html, 06.11.2004 [20]EPA – Environmental Protection Agency, http://www.epa.gov/airnow/, 06.11.2004 [21]Finzi G., Nunnari G., Air quality forecast and alarm systems, Air Quality Modelling- Theories, Methodologies, Computational Techniques and Available Databases and Software, Vol. II, Ch. 16, EnviroComp Institute and Air & Waste Management Association co-publishers, 2004 [22]Finzi G., Tebaldi G., A mathematical model for air pollution forecast and alarm in an urban area, Athmospheric Environment, vol. 16, 9, 9.Viri in literatura 81 2055-2059, 1982 [23]Gardner M. W., Dorling S. R., Artificial neural networks (the multilayer perceptron)-A review of applications in the atmosferic sciences, Atmospheric Environment Vol. 32, Elsevier Science Ltd., Great Britain, 1998 [24]Gardner M. W., Dorling S. R., Neural network modelling of the influence of local meteorology on surface ozone concentrations, Proceedings 1st International Conference on GeoComputation, University of Leeds, 359-370, 1996. [25]Geai P., Methode D'interpolation et de Reconstruction Tridimesionelle d'un Champ de Vent, le Code d'Analyse Objective MINERVE, Technical Report EDF/DER/HE/34-87.03 [26]Gibbs M. N., Bayesian Gaussian Processes for Regression and Classification, Ph.D. Disertation, Cambridge University, 1997. [27]Gutfreund P., Liu C., Nicholson B., Roberts E., COMPLEX I and II Model Performance Evaluation in Nevada and New Mexico, Journal of Air Pollution Control Association, Vol. 33, 846-871, 1983 [28]Haagen-Smit A. J., Chemistry and physiology of Los Angeles Smog. Industrial and Engineering Chemistry 44, 1342-1346, 1952 [29]Hornik K., Stinchcombe M., White H., Multi-layer feedforward networks are universal approximators, Neural Networks 2, 1989. [30]How Ozone Affects Bacteria, Fungus, Molds And Viruses, http://mold-help.org/content/view/436/, 11.12.2004 [31]IAEA Training course in Planning, Preparedness & Response to Radiological Emergencies, lecture 20.5.3, Basic meteorology – plume transport, Argonne national laboratory, Argonne, Illinois, United States of America, 1982 [32]Kloss S. D., A Discussion on Ozone Chemistry, http://www.environmental-expert.com/article316.htm, 06.11.2004 [33]Krušič M., Slovenija, Turistični vodnik, Mladinska knjiga, Ljubljana, 2002 [34]Lamb B., Gaseous pollutant characteristics. Handbook of air pollution technology. John Wiley and Sons, Inc. USA, 1984 [35]Lawrence J., Introduction to Neural Networks, California Scientific Software, Grass Valley, USA, 1991 82 9.Viri in literatura [36]Likar B., Prediktivno vodenje nelinearnih sistemov na osnovi Gaussovih procesov, Magistrsko delo, Fakulteta za elektrotehniko, 2004. [37]MacKay D. J. C., Hyperparameters: Optimize, or integrate out? In Maximum Entropy and Bayesian Methods, Santa Barbara 1993 , ed. by G. Heidbreder, pp. 43--60, Dordrecht, Kluwer, 1996. [38]Mestna občina Nova Gorica, Spletna stran mestne občine Nova Gorica, http://www.nova-gorica.si, 22.12.2004 [39]Miklavčič Brezigar I., Nova Gorica, Vodnik po mestu, IKI – Institut za komunikacije in informatiko, Ljubljana, 1998 [40]Mlakar P., Izbira značilk vzorcev za model kratkoročnega napovedovanja onesnaževanja ozračja, Doktorska disertacija, Fakulteta za elektrotehniko, Univerza v Ljubljani, 1996 [41]Moreno J., Moreno-Grau S., Bayo J., M. Angosto J., Jimenez J. E., Moreno-Clavel J., Photochemical pollution precursors and ozone in the athmosphere of Cartagena from 1995 to 1997. Wessex institute of techology, Air pollution 98, 1998. [42]Neal R. M., Bayesian Learning for Neural Networks, Number 118 in Lecture Notes in Statistics, New York, Springer, 1996. [43]NeuroShell 2, Users manual, Ward Systems Group, Inc., Frederick, Maryland, 1993 [44]Paine R. J., Eagan B. A., User's guide to Rough Terrain – Diffusion Model (RTDM), EPA/SW/MT-88/041a, U.S. EPA Research Triangle Park, NC, 1987 [45]Periodic elements, Oxygen, http://pearl1.lanl.gov/periodic/elements/ 8.html, 06.11.2004 [46]Perry S. G., CTDMPLUS: A Dispersion Model for Source near Complex Topography, Part I: Technical Formulation, Journal of Applied Meteorology, Vol. 31, 633-645, 1992 [47]Rasmussen C. E., Evaluation of Gaussian Processes and other Methods for Non-Linear Regression, Ph.D. Disertation, Graduate department of Computer Science, University of Toronto, Toronto, 1996. [48]Rasmussen C. E., Home page, Gaussian Processes for Regression software code, http://www.gatsby.ucl.ac.uk/~edward/code/gp/, 9.Viri in literatura 83 10.10.2004 [49]Rojko A., Položajno vodenje nelinearnih mehanizmov z uporabo mehke logike, Doktorat, Fakulteta za elektrotehniko, računalništvo in informatiko, Univerza v Mariboru, 2002 [50]Rumelhart D. E., G. E. Hinton, R. J. Williams, Learning internal representations by error propagation, Paraller distributed processing: Explorations in the Microstucture of Cognition, Vol. 1, MIT Press, Cambridge, MA, 1986. [51]Škrjanc I., Fuzzy modelling and identification, Fakulteta za elektrotehniko, Univerza v Ljubjani, 2003 [52]Takagi T., Sugeno M., Fuzzy Identification of Systems and Its Applications to Modeling and Control, IEEE Trans. on Systems, man and Cybernetics, vol.15, no.1, str. 116-132, januar 1985. [53]The Math Forum @ Drexel, Ask Dr. Math, Questions and Answers from our archives, Complexity of Matrix Inversion, http://mathforum.org/library/drmath/view/51908.html, 06.12.2005 [54]The MathWorks, Inc. Matlab® version 6.5. Matlab documentation. [55]Tinarelli G., et all, Lagrangian Particle Simulation of Tracer Dispersion in the Lee of a Shematic Two-dimensional Hill, Journal of Applied Meteorology, Vol. 33, No. 6, 744-756, 1994 [56]Tomšič G., Matematika IV, Fakulteta za elektrotehniko,1998. [57]Uredba o benzenu in ogljikovemu monoksidu v zunanjem zraku. Ur. list RS, št. 52/2002. [58]Uredba o ozonu v zunanjem zraku. Ur. list RS, št. 8/2003. [59]Uredba o ukrepih za ohranjanje in izboljšanje kakovosti zunanjega zraka. Ur. list RS, št. 52/2002. [60]Uredba o žveplovem dioksidu, dušikovih oksidih, delcih in svincu v zunanjem zraku. Ur. list RS, št. 52/2002. [61]Van Aalst R.M., F.A.A.M. de Leeuw (urednika), National Ozone Forecasting System And International Data Exchange In Northwest Europe, European Topic Centre on Air Quality, http://reports.eea.eu.int/TEC09/en/tech09.pdf, september 1997 [62]Wikipedia, http://www.wikipedia.org, 06.11.2004 [63]WMO/UNEP Scientific Assessment of Ozone Depletion: 2002, 84 9.Viri in literatura http://www.al.noaa.gov/WWWHD/Pubdocs/assessment02.html, 12.12.2004 [64]Yi J., Prybutok R., A neural network model forcasting for prediction of daily maximum ozone concentration in an industrialized urban area. Environmental Pollution 92(3), 349-357, 1996. [65]Zadeh L. A., Fuzzy logic = Computing with worlds, IEEE Trans. Fuzzy Syst., vol. 4, no.2, str. 103-110, Maj 1996. [66]Zadeh L. A., Fuzzy sets, Information and Control, 8, str. 338-353, 1965 85 Dodatek A: Osnove verjetnostnega modeliranja A.1. Naključna spremenljivka Naključna spremenljivka je spremenljivka, katere vrednost je odvisna od naključja. O naključnem procesu (ang. stochastic process) govorimo, ko imamo več neodvisnih (časovno, prostorsko itd.) realizacij naključne spremenljivke. Porazdelitve vrednosti posameznih realizacij naključne spremenljivke so lahko poljubne, opišemo jih s porazdelitvenim zakonom in zalogo vrednosti, s čimer je naključna spremenljivka tudi popolnoma določena. Glede na zalogo vrednosti ločimo diskretne in zvezne naključne spremenljivke, glede na obliko porazdelitvenega zakona pa različne standardne (enakomerna, normalna ali Gaussova, Studentova itd) in nestandardne porazdelitve [56,36]. Splošna oblika porazdelitvenega zakona za naključno spremenljivko X je porazdelitvena funkcija ali distribucija F(x): F{x)=P{X-ad) F{-oo)= Um F{x)=0 (A4) in v primeru gotovega dogodka (X->aó) F{oo) = limF{x) = 1 (A5) x—>x> . . 86 Dodatek A: Osnove verjetnostnega modeliranja A.2. Analiza porazdelitev naključnih spremenljivk Pri analizi porazdelitev naključnih spremenljivk si pomagamo z različnimi številskimi karakteristikami. Najpomembnejše so predstvljene v nadaljevanju. Matematično upanje ali povprečna vrednost naključne spremenljivke X, E(X): E{x)=] xp{x)dx . (A.6) — 30 Varianca ali disperzija naključne spremenljivke D(X) je definirana kot povprečje kvadratov odstopanj spremenljivke X od njene povprečne vrednosti E(X): var{X) = E{{X-E{X))2) = J {x-E{X))2v{x)dx = E{X2)-{E{X))2 (A.7) — GO Standardna deviacija je definirana kot pozitivni kvadratni koren variance naključne spremenljivke X: Xn) = P(X1X2...>Xn...temperatura rosišča. V nadaljevanju so prikazani različni modeli nelinearne funkcije temperature rosišča, ki so bili modelirani z Gaussovimi procesi. Skupne lastnosti vseh modelov so: • zgrajeni so bili pri konstantni temperaturi zraka Tz=20°C, • imajo samo eno vhodno spremenljivko, to je relativna vlaga zraka Rv, • za indentifikacijo je bilo uporabljeno enako število učnih točk, • izhod modela predstavlja temperatura rosišča v odvisnosti od relativne vlage zraka pri konstatni temperaturi zraka. Izhod vsakega modela je prikazan na dveh slikah. Na prvih slikah so prikazani odzivi modela: • z znakom o (krogec) so prikazane učne točke, • z znakom . (pika) pa izhodi modela, • pikčasta črta prikazuje vrednost dvojne standarne deviacije izhoda iz modela, 92 Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi • polna črta prikazuje temperaturo rosišča glede na enačbo (B.1), Druga slika pa prikazuje: • napako, ki je izračunana kot razlika med izhodom modela in pravo vrednostjo glede na enačbo (B.1), • pikčasta črta prikazuje vrednost dvojno standarno deviacijo izhoda iz modela, Prvi model je bil identificiran s pomočjo petih učnih točk, ki so bile naključno izbrane na območju relativne vlage od 20% do 100%. Prikazan je na slikah 76 in 77. Drugi model je bil identificiran s pomočjo petih učnih točk, ki so bile naključno izbrane na območju relativne vlage od 20% do 100%. Od prvega modela se razlikuje v tem, da so je bil učnim točkam dodan še šum. Model je prikazan na slikah 78 in 79. Tretji model je bil identificiran s pomočjo petih učnih točk, ki so bile izbrane na območju relativne vlage od 20% do 100% v enakomernih presledkih. Prikazan je na slikah 80 in 81. Četrti model je bil identificiran s pomočjo petih učnih točk, ki so bile izbrane na območju relativne vlage od 20% do 100% v enakomernih presledkih. Od tretjega modela se razlikuje v tem, da so je bil učnim točkam dodan še šum. Prikazan je na slikah 82 in 83. Peti model je bil identificiran s pomočjo petnajstih učnih točk, ki so bile izbrane na območju relativne vlage od 20% do 100% v enakomernih presledkih. V vsakem presledku so bile določene tri učne točke, ki so imele dodan šum. Prikazan je na slikah 84 in 85. Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi 93 40 60 Relativna vlaga zraka [%] Slika 76: Prvi model temperature rosišča 20 40 60 80 Relativna vlaga zraka [%] Slika 77: Napaka prvega modela temperature rosišča 94 Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi 20 40 60 80 Relativna vlaga zraka [%] Slika 78: Drugi model temperature rosišča -4[ 0 20 40 60 80 Relativna vlaga zraka [%] Slika 79: Napaka drugega modela temperature rosišča Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi 95 40 60 Relativna vlaga zraka [%] Slika 80: Tretji model temperature rosišča 20 40 60 80 Relativna vlaga zraka [%] Slika 81 : Napaka tretjega modela temperature rosišča 96 Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi 20 40 60 80 Relativna vlaga zraka [%] Slika 82: Četrti model temperature rosišča 20 40 60 80 Relativna vlaga zraka [%] Slika 83: Napaka četrtega modela temperature rosišča Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi 97 40 60 80 Relativna vlaga zraka [%] Slika 84: Peti model temperature rosišča 20 40 60 80 Relativna vlaga zraka [%] Slika 85: Napaka petega modela temperature rosišča 98 Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi B.2. Prikaz vpliva posameznih hiperparametrov modela Model, ki je bil identificiran s pomočjo petih učnih točk, ki so bile izbrane na območju relativne vlage od 20% do 100% v enakomernih presledkih in je prikazan na slikah 80 in 81, sem uporabil za praktičen prikaz vpliva posameznih hiperparametrov v modelu. Za modeliranje je bila uporabljena oblika kovariančne funkcije 1 N ixu,—x[,h2 Cf{xm,xn;0) = 0N+1e2'-1 20< +ömn@N+2 , (B.2) kjer N označuje dimezijo vhodnega vektorja. V praktičnem primeru je bila dimenzija N=l. To pomeni, da je bilo potrebno poiskati s pomočjo optimizacijskega algoritma, ki uporablja metodo največje podobnosti, optimalne vrednosti treh hiperparametrov. Po indentifikaciji modela so bile določene vrednosti hiperparametrov sledeče 01-0.00029, @2=471, ©3-0.071 . (B.3) Za prikaz pomena posameznega hiperparametra sem napravil tri poizkuse. V vsakem poizkusu sem spreminjal vrednost enega hiperparametra, vrednosti preostalih dveh hiperparametrov pa se med poizkusom niso spreminjale in so bile nastavljene na optimalne vrednosti. V prvem poizkusu, ki je prikazan na sliki 86 sem spreminjal vrednost prvega hiperparametra. Prvi hiperparameter &1 predstavlja dolžinsko merilo dimenzije. Iz slike je razvidno, da hiperparameter v bistvu določa področje (območje) vpliva posamezne učne točke. V primeru majhne vrednosti je območje zelo široko, kar v skrajnem primeru zelo majhne vrednosti pomeni, da bodo vse napovedane točke enake srednji vrednosti vseh učnih točk. Medtem, ko so v drugem skrajnem primeru zelo velike vrednosti hiperparametra napovedi z nizko varianco možne samo v neposredni bližini učnih točk. V drugem poizkusu, ki je prikazan na sliki 87 sem spreminjal vrednost drugega hiperparametra. Drugi hiperparameter 02 določa skupen vpliv vseh učnih točk na napovedane vrednosti. V skrajnem primeru zelo majhne vrednosti hiperparametra imajo učne točke skoraj ničen vpliv na napovedi (vse vrednosti napovedi so blizu konstantne vrednosti 0). V nasprotnem skrajnem primeru velike vrednosti hiperparametra pa imajo učne točke zelo velik vpliv na napovedane vrednosti v taki meri, da imajo napovedi v točkah, ki so malce oddaljene od učnih točk, že zelo veliko negotovost. Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi 99 V tretjem poizkusu, ki je prikazan na sliki 88 sem spreminjal vrednost tretjega hiperparametra. Tretji hiperparameter ?3 določa minimalno negotovost napovedanih točk, kar pomeni, da nobena napovedana vrednost ne more imeti manjše negotovosti. V skrajen primeru zelo majhne vrednosti hiperparametra imajo tudi vse napovedi zelo majhno negotovost, v drugem skrajem primeru velike vrednosti hiperparametra pa imajo vse napovedi izredno velike negotovosti, vrednosti napovedi pa blizu konstante vrednosti 0. 100 Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi Slika 86: Vpliv prvega hiperparametra na model Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi 101 Slika 87: Vpliv drugega hiperparametra na model 102 Dodatek B: Preprost primer modeliranja nelinearne funkcije z Gaussovimi procesi Slika 88: Vpliv tretjega hiperparametra na model 103 Dodatek C: Odvisnost negotovosti od kovariančne funkcije Negotovost napovedi, ki je odvisna od oblike kovariančne funkcije, opisuje enačba = K-kT N+1 C N 1k N+1 , (C.1) kjer je pomen oznake K = C(xN+1,xN+1;0) (C.2) ter oznake kN+=(C(x1xN+1;0),...,C(xN,xN+1;0)) . (C.3) V nadaljevanju bom predstavil, kako se negotovost napovedi spreminja v odvisnosti od oddaljevanja iz področja učnih točk. V prvem primeru je bila za predstavitev uporabljena kovariančna funkcija brez modela šuma, v drugem primeru pa kovariančna funkcija z modelom šuma. C.1. Uporaba stacionarne kovariančne funkcije brez upoštevanja šuma Če uporabimo pri modeliranju z Gaussovimi procesi kovariančno funcijo oblike Cf{xm,xn;0) = 01e 2'=1 2r' , (C.4) ki ne upošteva šuma, se negotovost napovedi z oddaljevanjem iz področja učnih točk ekponentno povečuje od vrednosti 0 do neke končne vrednosti, ki je enaka <9,. V takšnem primeru stacionarne kovariančne funkcije je vrednost vseh diagonalnih elementov kovariančne matrike (elementi predstavljajo varianco posameznih učnih točk) enaka <97 (saj je vrednost ekponentnega dela v enačbi (C.4) enaka 1, ker je oddaljenost med točkama enaka 0). Poleg tega pa je tudi vrednost k enaka <9,. V enem skrajnem primeru, ko je vhod v model enak enemu izmed učnih vhodov (na vhod modela pripeljemo vrednost vhoda s katero smo ga pred tem učili), postane en element kovariančnega vektorja kN enak varianci učne točke, ostali elementi pa so enaki kovariancam med učnimi točkami in napovedano točko kN=[C{x1xN+1;0),...,01...,C{xN,xN+1;©)] . (C.5) 104 Dodatek C: Odvisnost negotovosti od kovariančne funkcije Vrednost omenjenega elementa je enaka 0j (saj je vrednost ekponentnega dela v enačbi (C.4) enaka 1, ker je oddaljenost med točkama enaka 0). Kovariance med učnimi točkami in napovedano točko so manjše od <97 in se ekponentno zmanjšujejo z oddaljenostjo od napovedovane točke. Zaradi enakosti vhoda, za katerega želimo napoved, z enim od učnih vhodov, postane tudi vrednost izraza kTN+lC^kN+l enaka <9,. Iz tega pa sledi, da je vrednost negotovosti napovedi, kot ga opisuje enačba (C.1), enaka 0. To je povsem razumljivo, saj model lahko v točki s katero je bil učen, s popolno gotovostjo napove izhod. Vrednost izraza kTN+lC~Nl k N+l se z oddaljevanjem od učnih točk ekponentno zmanjšuje in doseže v neki veliki oddaljenosti končno vrednost, ki znaša 0. Zaradi tega pa negotovost eksponentno narašča proti vrednosti &j. C.2. Uporaba stacionarne kova ria nč ne funkcije z upoštevanjem šuma Predpostavljajmo, da je šum naključen. Iz tega sledi, da ne pričakujemo nobene korelacije med šumom in določenimi izhodi ter, da drugi izraz vpliva samo na diagonalne elemente kovariančne matrike. Varianco šuma označimo z &2, za katero se predpostavlja, da je v primeru stacionarne kovariančne funkcije neodvisna od vhodov. Z ômn pa je označena Kronceker-jeva delta funkcija amn=\0'm*n\ . (C.6) Če uporabimo pri modeliranju z Gaussovimi procesi kovariančno funcijo oblike Cf{xm,xn;&) = &1e 2,=1 2r ' +5mn@2 , (C.7) ki upošteva šum, se negotovost napovedi z oddaljevanjem iz področja učnih točk ekponentno povečuje od vrednosti malce večje od 02 do neke končne vrednosti, ki je enaka &+&2 . V takšnem primeru stacionarne kovariančne funkcije je vrednost vseh diagonalnih elementov kovariančne matrike (elementi predstavljajo varianco posameznih učnih točk) enaka 0+02 (saj je vrednost ekponentnega dela v enačbi (C.7) enaka 1, ker je oddaljenost med točkama enaka 0, vrednost Kronecker-jeve delta funkcije pa je enaka 1). Dodatek C: Odvisnost negotovosti od kovariančne funkcije 105 Poleg tega pa je vrednost ? v tem primeru enaka ?1+?2 (drugačna kot v prejšnjem primeru). Zopet v enem skrajnem primeru, ko je vhod v model enak enemu izmed učnih vhodov (na vhod modela pripeljemo vrednost vhoda s katero smo ga pred tem učili), postane en element kovariančnega vektorja kN enak varianci učne točke, ostali elementi pa so enaki kovariancam med učnimi točkami in napovedano točko kN=[C{xlxN+l;0),...,01...,C{xN,xN+l;©)] . (9.1) Vrednost omenjenega elementa je enaka ?1 (saj je vrednost ekponentnega dela v enačbi (C.7) enaka 1, ker je oddaljenost med točkama enaka 0, vrednost Kronecker-jeve delta funkcije pa je enaka 0). Zaradi enakosti vhoda, za katerega želimo napoved, z enim od učnih vhodov, postane tudi vrednost izraza kTN+lCNkN+l blizu vrednosti ?1. To pomeni, • da je v skrajnem primeru, ko nimamo šuma, izraz kTN+lC^kN+l enak ?1, • da je v primeru majhne vrednosti šuma ?2 izraz kTN+lC~NlkN+l malce manjši od ?1, • v primeru večje vrednosti šuma ?2, pa postaja vrednost izraza kTN+lCNkN+l čedalje manjša od ?1,. Glede na enačbo (C.1) in dejstvo, da je ? v tem primeru enaka ?1+?2, sledi, da je vrednost negotovosti napovedi, kot ga opisuje enačba (C.1), za faktor, ki je funkcija parametrov ?1 in ?2 ter kovarianc med učnimi točkami, večja od ?2. Tudi to je povsem razumljivo, saj model v točki s katero je bil učen, ne more izdelati napovedi z manjšo negotovostjo, kot je negotovost zaradi šuma. Vrednost izraza kTN+lC~NlkN+l se z oddaljevanjem od učnih točk ekponentno zmanjšuje in doseže v neki veliki oddaljenosti končno vrednost, ki znaša 0. Zaradi tega pa negotovost eksponentno narašča proti vrednosti ?1+?2. 107 Izjava Izjavljam, da sem magistrsko nalogo izdelal samostojno pod vodstvom prof. dr. Marka Muniha in prof. dr. Draga Matka. Izkazano pomoč ostalih sodelavcev sem v celoti navedel v zahvali. Boštjan Grašič