7 Uporaba kopul v hidrologiji Nejc Bezak * , Mitja Brilly * , Matjaž Mikoš * in Mojca Šraj * Povzetek Funkcije kopula se v svetu vse pogosteje uporabljajo na različ nih znanstvenih področ jih, v hidrologiji pa se njihova uporaba v več ji meri pojavlja šele v zadnjem desetletju. V prispevku so prikazani rezultati nekaterih praktič nih možnosti uporabe kopul v hidrologiji. Naredili smo bivariatne ter trivariatne verjetnostne analize elementov visokovodnih valov (konice pretokov, volumni ter č asi trajanja), analizirali smo povezanost pretokov, volumnov valov ter vrednosti koncentracij suspendiranih snovi ter definirali model, s katerim lahko na podlagi znanih vrednosti pretokov in padavin ocenimo vrednosti koncentracij suspendiranih snovi. Z uporabo statistič nih testov smo pokazali, da je kopula model dal boljše rezultate kot nekateri pogosto uporabljeni regresijski modeli. S tem smo v oceno koncentracij suspendiranih snovi vpeljali dodatne informacije (padavine) in dobili ocenjene vrednosti bližje dejanskim, izmerjenim vrednostim. Pri multivariatnih verjetnostnih analizah pa smo izrač unali različ ne skupne povratne dobe. Ključ ne besede: Kopula, bivariatne analize, trivariatne anlize, pretok, volumen vala, suspendirane snovi. Keywords: Kopula, bivariate analysis, trivariate analysis, discharge, wave volume, suspended solids. Uvod V hidrologiji so se do sedaj za izvedbo verjetnostnih analiz več inoma uporabljale univariatne porazdelitvene funkcije, kjer upoštevamo le eno spremenljivko (več inoma je to konica pretoka). V zač etku 21. stoletja so se v hidrologiji pojavili prvi č lanki, kjer so raziskovalci za izvedbo različ nih analiz uporabili funkcijo kopula (Favre et al., 2004; Salvadori & De Michele, 2004). V drugih vedah, kot sta ekonomija ali biologija, so se kopule uporabljale že nekoliko prej. Teorija kopul temelji na matematič nem teoremu, ki ga je že leta 1959 predstavil Sklar (Sklar, 1959), in bo nekoliko podrobneje predstavljen v nadaljevanju. Temeljno literaturo o funkcijah kopula pa predstavlja naslednje gradivo: Joe (1997), Nelsen (2006), Salvadori et al. (2007). Kopule omogoč ajo izgradnjo multivariatnega modela, kjer hkrati upoštevamo dve ali več v naravi odvisnih spremenljivk. Tako lahko pri verjetnostnih analizah poleg konic pretokov hkrati upoštevamo tudi volumne in č ase trajanja visokovodnih valov ali druge spremenljivke. Ta koncept lahko uporabimo pri vseh multivariatnih problemih, kjer nastopa d medsebojno bolj ali manj odvisnih spremenljivk, kjer odvisnost določ a koeficient korelacije (Kendallov ali Spearmanov) oziroma parameter kopule. Tako so bile kopule uporabljene za izvedbo multivariatnih verjetnostnih analiz visokovodnih valov (Favre et al., 2004; Salvadori & De Michele, 2004; Ganguli & Reddy, 2013), verjetnostne analize padavin (Zhang & Singh, 2007), analize suše (Shiau et al., 2007; Wong et al., 2010; Ma et al., 2011), geostatistič ne interpolacije kot alternativa obič ajnemu krigingu (Bardossy & Li, 2008), preverjanje ustreznosti prelivnega objekta na jezu (De Michele et al., 2005) ter še pri mnogih drugih * Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo, Jamova 2, 1000 Ljubljana, Slovenija 8 hidroloških problemih. Na spletni strani mednarodnega hidrološkega društva STAHY (angl. International Association of Hydrological Sciences; www.stahy.org) so zbrani č lanki, kjer je bila za izvedbo hidroloških analiz uporabljena funkcija kopula. Da se kopule v hidrologiji vse pogosteje uporabljajo prič a tudi dejstvo, da se je število č lankov s to metodologijo v bazi Web of knowledge (kategorija Water resources), od leta 2004 povzpelo na več kot 180, v še ne konč anem letu 2013 pa je bilo do sedaj objavljenih največ prispevkov v zvezi s kopulami (do oktobra je bilo v letu 2013 indeksiranih 37 č lankov). Namen prispevka je prikaz praktič ne uporabe kopul v hidrologiji. Prikazan in opisan je postopek uporabe bivariatne ter trivariatne verjetnostne analize visokovodnih valov na primeru hidrološke postaje Litija na reki Savi. Funkcijo kopula smo uporabili tudi za modeliranje pretokov, volumnov valov ter koncentracij suspendiranih snovi. Za izvedbo te analize so uporabljeni podatki z vodomerne postaje Gornja Radgona na reki Muri, kjer so se več kot 25 let neprekinjeno izvajale dnevne meritve koncentracij suspendiranih snovi (ARSO, 2013). Poleg tega smo poskušali funkcijo kopula uporabiti za ocenjevanje vrednosti suspendiranih snovi na podlagi izmerjenih vrednosti pretokov in padavin. Za ta namen smo uporabili podatke z vodomerne postaje Ranca na reki Pesnici ter padavinske postaje Polič ki vrh. Na teh praktič nih primerih so opisani osnovni koraki za izvedbo analiz z uporabo funkcij kopula. Podatki in metode Za izvedbo različ nih analiz smo uporabili podatke s hidroloških postaj Litija na reki Savi (vzorec letnih maksimalnih pretokov), Gornja Radgona na reki Muri (vzorec letnih maksimalnih pretokov in suspendiranih snovi), Ranca na reki Pesnici (meseč ne vsote pretoka in suspendiranih snovi) ter podatke s padavinske postaje Polič ki vrh (meseč ne vsote padavin), ki leži blizu vodomerne postaje Ranca. Preglednica 1 prikazuje nekatere osnovne znač ilnosti obravnavanih postaj ter lastnosti obravnavanih vzorcev. Zaradi različ ne narave obravnavanih problemov smo pri različ nih postajah obravnavali različ na obdobja meritev. Pri analizah koncentracij suspendiranih snovi smo bili omejeni z relativno kratkimi nizi meritev (Bezak et al., 2013b). Več informacij o podatkih in lastnostih vzorcev so podali Bezak et al., 2013a (postaja Litija) ter Bezak et al., 2013b (postaji Gornja Radgona in Ranca). Preglednica 1: Pregled obravnavanih postaj in nekatere osnovne znač ilnosti pripadajoč ih vzorcev Postaja Litija Gornja Radgona Ranca Polič ki vrh Reka / Območ je Sava Mura Pesnica Poreč je Pesnice Tip postaje vodomerna vodomerna vodomerna padavinska Obdobje 1953-2010 1977-2005 1970-1973 1970-1973 Tip vzorca letni maksimumi letni maksimumi meseč ne vsote meseč ne vsote Opazovane spremenljivke pretok pretok, suspendirane snovi pretok, suspendirane snovi padavine Velikost vzorca (število podatkov) 58 29 45 45 Za izvedbo analiz smo torej uporabili funkcije kopula. Osnovni teorem je podal Sklar (1959): , … , = ,… , , (1) 9 kjer je C d dimenzionalna funkcija kopula. Multivariatno d dimenzionalno porazdelitev lahko zapišemo kot kombinacijo funkcije kopula C in robnih porazdelitvenih funkcij F 1 ,…, F d (angl. marginal distributions) (enač ba 1). Gre v bistvu za obič ajne univariatne porazdelitvene funkcije, kjer lahko uporabimo različ ne parametrič ne (npr. normalna, log- normalna, Pearsonova III, generalizirana porazdelitev ekstremnih vrednosti,…) ali neparametrič ne (kernelova gostota; angl. kernel density) porazdelitvene funkcije. Glavna prednost funkcij kopula pred obič ajnimi multivariatnimi porazdelitvenimi funkcijami (npr. multivariatna normalna porazdelitev) je prav loč enost robnih porazdelitev in funkcij kopula. Ocenjevanje parametrov robnih porazdelitev (obič ajen postopek) in parametrov funkcije kopula sta tako dva loč ena postopka. To dejstvo je pomembno zato, ker lahko kot robne porazdelitve posameznih spremenljivk izberemo različ ne porazdelitvene funkcije. Tako lahko npr. pri analizah konic, volumnov in č asov trajanja visokovodnih valov kot robne porazdelitve izberemo različ ne porazdelitvene funkcije, npr. Gumbelovo za konice, Pearsonovo III za volumne ter normalno za č ase trajanja. Sklarov teorem (Sklar, 1959) pa omogoč a, da nato ocenimo parametre kopule povsem loč eno od robnih porazdelitvenih funkcij. Slednje pomeni, da lahko odvisnost med izbranimi spremenljivkami preuč ujemo loč eno od robnih porazdelitev. Za modeliranje različ nih problemov se uporabljajo različ ne funkcije kopula iz različ nih družin. Tako se v hidrologiji najpogosteje uporabljajo kopule iz Arhimedove družine (angl. Arhcimedean family), kopule iz družine ekstremnih vrednosti (angl. extreme value family), nekoliko manj pogosto pa tudi kopule iz eliptič ne družine (angl. elliptical family), ki pa se npr. v ekonomiji uporablja pogosteje (Grimaldi & Serinaldi, 2006; Genest & Favre, 2007; Ganguli & Reddy, 2013). Za modeliranje suše se najpogosteje uporablja kopula Clayton iz Arhimedove družine (Ma et al., 2011), za analize visokovodnih valov se pogosto uporabijo kopule iz družine ekstremnih vrednosti (npr. kopuli Gumbel-Hougaard ali Galambos) (Genest & Favre, 2007), za geostatistič ne interpolacije se pogosto uporabljajo normalne kopule iz eliptič ne družine (Bardossy & Li, 2008). Simetrič ne kopule iz Arhimedove družine, ki se najpogosteje uporabljajo, imajo več inoma le en parameter. Ta parameter lahko podobno kot pri univariatnih porazdelitvah ocenimo z različ nimi metodami. Uporabimo lahko metodo momentov (Kendallov ali Spearmanov koeficient korelacije), metodo največ jega verjetja (angl. maximum likelihood method), pseudo metodo največ jega verjetja (angl. pseudo likelihood method), kjer naj bi bila slednja najprimernejša (Joe, 1997; Nelsen, 2006; Salvadori et al., 2007). V primeru, da obravnavamo vsaj trivariaten primer, lahko analize naredimo tudi z asimetrič nimi kopulami, kjer imajo te en parameter več kot simetrič ne različ ice kopul (Grimaldi & Serinaldi, 2006; Hofert & Mä chler, 2011). Slednje se izkažejo za posebej uporabne v primerih, ko je korelacija med enim parom spremenljivk izrazito več ja kot korelaciji med drugima dvema paroma (trivariaten primer). Podobno kot pri univariatnih porazdelitvenih funkcijah so bili tudi za primer funkcij kopula razviti nekateri statistič ni testi za primerjavo različ nih teoretič nih modelov in izbiro najustreznejše kopule. V literaturi prevladujeta dva tipa testov, najpogostejši so testi, ki temeljijo na primerjavi teoretič ne ter empirič ne funkcije kopula, nekoliko manj pogosto pa se uporabljajo testi, ki temeljijo na Kendallovi porazdelitveni funkciji (Genest et al., 2009; Kojadinovic et al., 2011). Za kontrolo ustreznosti izbrane kopule pa je priporoč ljivo uporabiti tudi grafič ne teste, s katerimi lahko potrdimo ali zavrnemo ugotovitve statistič nih testov (Genest & Favre, 2007). Opisane korake (določ itev robnih porazdelitvenih funkcij, ocenjevanje parametrov kopul, izbira najustreznejše kopule z uporabo statistič nih ter grafič nih testov) je potrebno izvesti pri vseh analizah s kopulami, nato pa se postopki razlikujejo od namena uporabe 10 modela. Pri verjetnostnih analizah moramo izrač unati različ ne povratne dobe, pri ocenjevanju vrednosti pa moramo izvesti simulacije. Oba koraka bosta malce podrobneje opisana v naslednjih dveh odstavkih. Pri multivariatnih verjetnostnih analizah z uporabo funkcij kopula lahko izrač unamo različ ne vrednosti pogojnih (angl. conditional) ali skupnih (angl. joint) povratnih dob (Salvadori et al., 2007; Graler et al., 2013). Skupni povratni dobi imenovani OR ter AND sta definirani z naslednjima izrazoma (bivariatni primer): , = , , (2) , = , , kjer sta u in v izbrani robni porazdelitveni funkciji, pa predstavlja teoretič no funkcijo kopula (v tem primeru bivariatno). Ker pa gre za analize v več dimenzijah, te definicije povratnih dob seveda niso primerljive s tistimi, ki jih več inoma uporabljamo pri univariatnih analizah. To lastnost kopul bi morda lahko ocenili kot največ jo težavo, saj je težko (kar naenkrat) v praktič no uporabo vpeljati nov koncept povratnih dob. Tako so se v literaturi že pojavili poskusi določ itve t.i. kritič ne plasti (angl. critical layer), kjer gre za matematič no definicijo prej omenjenega pojma, ki je (bolj) analogna z obič ajnimi povratnimi dobami (Salvadori et al., 2011). Npr. pri trivariatnem primeru verjetnostnih analiz visokovodnih valov izhajamo iz izbrane vrednosti verjetnosti in na podlagi te vrednosti s pomoč jo simulacij določ imo ploskev v kopula prostoru, kjer vse toč ke, ki ležijo nad to ploskvijo, označ ujejo kritič ne dogodke (Salvadori et al., 2011). To ploskev lahko nato z uporabo inverznih porazdelitvenih funkcij transformiramo v realni prostor. Povratna doba OR predstavlja primer, ko se zgodi vsaj ena izmed obravnavanih spremenljivk, povratna doba AND pa primer, ko se hkrati zgodijo vse obravnavane spremenljivke. Iz definicije sledi, da je verjetnost pojava povratne dobe AND precej manjša od verjetnosti pojava povratne dobe OR, kar seveda pomeni, da je povratna doba za AND primer precej več ja (obravnavamo multivariaten primer). Funkcije kopula lahko uporabimo tudi za ocenjevanje vrednosti različ nih hidroloških spremenljivk na podlagi podatkov, ki jih imamo na razpolago. Prikazana bo praktič na uporaba v primeru ocenjevanja vrednosti koncentracij suspendiranih snovi. Pri modelu, kjer bi vrednosti koncentracij suspendiranih snovi ocenili le na podlagi pretokov, lahko izhajamo iz naslednje enač be: = { ≤ |! = "} = $ ,%%& $’ , (3) kjer moramo z uporabo statistič nih ter grafič nih testov ponovno izbrati ustrezno funkcijo kopula in kjer funkcija f vedno zavzame vrednost med 0 in 1. To pomeni, da lahko to funkcijo zamenjamo z naključ no generiranimi števili, ki pripadajo enakomerni porazdelitveni funkciji. Tako pridemo do enač be, kjer imamo le eno neznanko, to je vrednost porazdelitvene funkcije za spremenljivko SSC. V več ini primerov, zaradi kompleksnosti parcialnih odvodov ne moremo poiskati analitič ne rešitve enač be 3, lahko pa z uporabo numerič nih metod (npr. Newton-ova metoda) najdemo nič le, ki kot že reč eno ležijo med 0 ter 1. Te nič le pa so v bistvu vrednosti porazdelitvenih funkcij, ki jih želimo oceniti. Ker smo generirali veliko število podatkov (npr. 10000), dobimo raztros možnih vrednosti. Nato lahko določ imo najverjetnejšo vrednost koncentracije (glede na maksimum gostote verjetnosti izbrane neparametrič ne funkcije) in oblikujemo tudi empirič ne intervale zaupanja (npr. 10, 20 ali 50 % intervali). Te mejne vrednosti nato še transformiramo v realni prostor z uporabo inverzne oblike porazdelitvene funkcije. Tudi v tem primeru za ocenjevanje vrednosti koncentracij uporabimo le vrednosti pretokov, to pomeni, da ne 11 moremo prič akovati veliko boljših rezultatov kot pri preprostejših regresijskih modelih (npr. linearna ali eksponentna funkcija), ki se več inoma uporabljajo za ocenjevanje koncentracij suspendiranih snovi. Zato smo se odloč ili, da bomo poskušali oceno izboljšati tako, da bomo v analizah upoštevali še vrednosti padavin. Model, kjer bi za ocenjevanje vrednosti koncentracij suspendiranih snovi uporabili tudi vrednosti padavin, ima naslednjo obliko: ,( = { ≤ | ! = " | = )} = $ * ( ,(,%%&) $’$+ / $ * (,() $’$+ . (4) Koncept določ itve najverjetnejše vrednosti koncentracije suspendiranih snovi je enak kot pri uporabi bivariatnih kopul (enač ba 3). Rezultati uporabe kopul na praktič nih primerih V naslednjih podpoglavjih je prikazana uporaba funkcij kopula na štirih praktič nih primerih. Prikazani so rezultati bivariatnih ter trivariatnih verjetnostnih analiz visokovodnih valov, rezultati trivariatnih analiz konic pretokov, volumnov valov ter koncentracij suspendiranih snovi (SSC) in rezultati ocenjevanja vrednosti suspendiranih snovi na podlagi podatkov o pretokih in padavinah. Bivariatne verjetnostne analize visokovodnih valov Za izvedbo bivariatnih verjetnostnih analiz visokovodnih valov smo uporabili podatke z vodomerne postaje Litija na reki Savi. Uporabili smo 58 vrednosti letnih maksimumov. Č e želimo določ iti vrednosti volumnov in č asov trajanja valov, moramo najprej izloč iti bazni odtok. Za izloč anje baznega odtoka smo uporabili grafič no tro-toč kovno metodo (Šraj & Bezak, 2013; Šraj et al., 2013b). Naredili smo bivariatne analize parov podatkov: konica pretoka - volumen vala (Q-V), konica pretoka - č as trajanja vala (Q-D) ter volumen vala - č as trajanja vala (V-D). Za modeliranje konic pretokov in č asov trajanja visokovodnih valov smo na podlagi statistič nih ter grafič nih testov izbrali log-Pearsonovo III porazdelitev, za opis volumnov valov pa Pearsonovo III porazdelitev. Za ocenjevanje parametrov univariatnih porazdelitev smo uporabili metodo momentov L (Hosking & Wallis, 1997), za izbiro najustreznejših porazdelitvenih funkcij pa različ ne statistič ne ter grafič ne teste (Šraj et al., 2013b). Preglednica 2 prikazuje vrednosti Pearsonovega, Kendallovega ter Spearmanovega koeficienta korelacije za obravnavane pare podatkov. Korelacijski koeficienti določ ajo odvisnost med dvema spremenljivkama, torej č e je vrednost koeficienta blizu vrednosti 1 to pomeni, da sta spremenljivki odvisni, č e je koeficient enak 0, to pomeni neodvisnost spremenljivk. Korelacijske koeficiente za opis odvisnosti se lahko uporabi tudi za primer, kjer imamo več kot 2 spremenljivki, vendar se izrazi za izrač un korelacijskih koeficientov spremijo in zato je bolje za opis odvisnosti uporabiti parameter kopule, ki ga lahko v bivariatnem primeru ocenimo na podlagi Kendallovega koeficienta korelacije ali pa npr. s pseudo metodo največ jega verjetja, ki jo je nato lažje uporabiti v multivariatnem primeru. Od vrednosti korelacije oz. dejanske odvisnosti, ki jo opiše korelacijski koeficient, je odvisno katero funkcijo kopula bomo uporabili za analize. Nekatere kopule so primerne le za pozitivne vrednosti (npr. Gumbel- Hougaard) odvisnosti, spet druge le za korelacije blizu 0 (npr. Ali-Mikhail-Haq kopula). Torej izrač unana vrednost korelacijskega koeficienta določ a, katere funkcije kopula lahko uporabimo, nadalje pa tudi kakšen bo odziv modeliranih spremenljivk, ki sta medsebojno (ne)odvisni. Iz preglednice 2 lahko vidimo, da je korelacija med konicami pretokov in č asi 12 trajanja visokovodnih valov negativna. To pomeni, da za izvedbo analize (za ta par) ne moremo uporabiti vseh funkcij kopula, ki smo jih uporabili pri drugih dveh parih, saj so nekatere primerne le za pozitivne vrednosti odvisnosti. Obravnavali smo različ ne funkcije iz Arhimedove, eliptič ne družine ter družine kopul ekstremnih vrednosti. Na podlagi CramØr von-Mises ter Kolmogorov-Smirnov statistič nih testov (Genest & Favre, 2007) ter različ nih grafič nih testov smo kot najprimernejšo funkcijo določ ili kopulo Gumbel-Hougaard iz Arhimedove družine kopul. Ta kopula spada tudi v družino kopul ekstremnih vrednosti. Te ugotovitve veljajo za para Q-V ter V-D, medtem ko smo pri paru Q-D kot najprimernejšo določ ili Student-t kopulo iz eliptič ne družine kopul. Za ocenjevanje parametrov kopul smo uporabili metodo momentov (Kendallov koeficient korelacije). Definicije teoretič nih kopul in izraze, ki povezujejo parametre kopul s koeficienti korelacije (Kendall ali Spearman) lahko najdemo v osnovni literaturi (Nelsen, 2006; Salvadori et al., 2007). Preglednica 3 prikazuje izrač unane skupne povratne dobe za OR ter AND primera. Izrač unali pa smo tudi univariatne (obič ajne) povratne dobe za pretoke, volumne ter č ase trajanja za mediane in maksimalne vrednosti (preglednica 3). Preglednica 2: Vrednosti Pearsonovega, Kendallovega ter Spearmanovega koeficienta korelacije za obravnavane pare podatkov Pearson Kendall Spearman Q-V 0,52 0,39 0,54 V-D 0,68 0,48 0,63 Q-D -0,15 -0,08 -0,14 Preglednica 3: Izrač unane vrednosti univariatnih in skupnih bivariatnih povratnih dob [leta] za tri funkcije kopula za vsak posamezen par spremenljivk Univariatni primer Bivariatni primer Spr. -. Par Kopula , , Mediana Maks. Mediana Maks. Mediana Maks. Q 2,0 69,5 Q-V Gumbel 1,2 40,7 2,8 132,8 V 1,9 56,4 Frank 1,5 32,1 2,7 1035,5 D 2,1 70,7 Tawn 1,5 40,7 2,8 133,0 V-D Gumbel 1,6 43,6 2,8 111,7 Frank 1,6 32,6 2,7 801,4 Normal 1,6 36,4 2,8 224,6 Q-D Clayton 1,3 35,2 4,5 5700,3 AMH 1,3 35,1 4,6 8125,7 Student-t 1,3 36,8 4,5 708,6 Trivariatna verjetnostna analiza visokovodnih valov Tudi pri trivariatnih analizah visokovodnih valov smo uporabili podatke o letnih maksimumih z vodomerne postaje Litija na reki Savi. Analizirali smo maksimalne konice pretokov in pripadajoč e vrednosti volumnov in č asov trajanja visokovodnih valov. To pomeni, da so le konice pretokov maksimalne v vseh letih, volumni ter č asi trajanja pa so bili določ eni glede na vrednosti pretokov in niso nujno letni maksimumi. Uporabili smo kopule Gumbel-Hougaard, Frank, Joe ter Clayton iz Arhimedove družine in normalno ter 13 Student-t kopulo iz eliptič ne družine. Za ocenjevanje parametrov kopul smo uporabili metodo največ jega verjetja (Nelsen, 2006; Salvadori et al., 2007), medtem ko so robne porazdelitvene funkcije ostale enake kot pri bivariatnih analizah. Slika 1: Grafič na testa za preverjanje ustreznosti kopule Gumbel-Hougaard pri trivariatni analizi visokovodnih valov Save v Litiji Slika 1 prikazuje rezultate grafič nega testa za ustreznost kopule Gumbel-Hougaard iz Arhimedove družine kopul, kjer rdeč i znaki x označ ujejo dejanske vrednosti letnih maksimumov, č rne oznake pa so generirane vrednosti, ki so bile določ ene na podlagi popolnoma definiranega modela. Pri izrisu slike 1b smo uporabili tudi robne porazdelitvene funkcije. Na podlagi grafič nih ter statistič nih testov (CramØr von-Mises) smo kot najprimernejše določ ili kopule Gumbel-Hougaard, Frank (Arhimedova družina) ter normalno (eliptič na družina). Izrač unali smo tudi dve skupni povratni dobi imenovani OR ter AND (enač ba 2), ki pa se zaradi tega, ker obravnavamo tri spremenljivke (trivariatni primer) hkrati, malce spremenita (Ganguli & Reddy, 2013) Slika 2 prikazuje vrednosti pogojne povratne dobe OR, ki označ uje primer, ko se zgodi vsaj ena izmed obravnavanih spremenljivk. Prikazani so rezultati za kopuli Frank ter normalno. AND primer pa pomeni, da se hkrati zgodijo vse tri vrednosti spremenljivk. Poleg rezultatov za povratno dobo OR so na sliki 2 prikazane tudi vrednosti letnih maksimumov, ki smo jih uporabili v analizah. Opazimo lahko, da razlika med prikazanima kopulama ni izrazita. Vidimo lahko tudi, da ima več ina letnih maksimumov povratno dobo OR manjšo od 1,5 let, le za en dogodek pa je ta več ja od 3,5 let, kar si lahko razlagamo tako, da je verjetnost, da se zgodi vsaj ena izmed spremenljivk (konice, volumni ali č asi trajanja) relativno velika in posledič no je multivariatna vrednost skupne povratne dobe OR dokaj majhna. 14 Slika 2: Prikaz izmerjenih podatkov in skupnih povratnih dob (OR primer) za normalno ter Frank kopulo za trajanje visokovodnega vala (D) 15 dni za vodomerno postajo Litija na reki Savi Trivariatna verjetnostna analiza konic pretokov, volumnov in koncentracij suspendiranih snovi Naslednji problem, ki smo se ga lotili z uporabo funkcij kopula, je analiza konic pretokov (Q), pripadajoč ih volumnov valov (V) ter pripadajoč ih vrednosti koncentracij suspendiranih snovi (SSC) (Bezak et al., 2013b). Za izvedbo teh analiz smo uporabili podatke z vodomerne postaje Gornja Radgona na reki Muri (29 vrednosti letnih maksimumov), ki ima najdaljši niz meritev suspendiranih snovi, ki jih je izvajala Agencija RS za okolje (ARSO, 2013). Postopek analize je enak kot pri trivariatnih analizah visokovodnih valov (prejšnje poglavje). Za izloč anje baznega odtoka smo tokrat uporabili dvo-toč kovno metodo. Kot robno porazdelitev konic smo izbrali Gumbelovo porazdelitev, za modeliranje volumnov smo določ ili log-normalno porazdelitveno funkcijo, za opis koncentracij suspendiranih snovi pa generalizirano Pareto porazdelitev. Parametre univariatnih porazdelitev smo ponovno ocenili z uporabo metode momentov L (Hosking & Wallis, 1997). Izrač unali smo vrednosti Kendallovega koeficienta korelacije za pare podatkov: Q-V 0,56; Q-SSC 0,50 in V-SSC 0,40. Ker so vse tri vrednosti korelacij statistič no znač ilne (stopnja znač ilnosti 0,05) in razlike med njimi niso izrazite, smo se odloč ili za izvedbo verjetnostnih analiz z uporabo simetrič nih trivariatnih funkcij kopula (alternativna možnost bi bile asimetrič ne kopule). Za ocenjevanje parametrov smo tokrat uporabili pseudo metodo največ jega verjetja (Salvadori et al., 2007; Ma et al., 2011). Z uporabo grafič nih ter statistič nih testov (CramØr-von Mises) smo tudi tokrat poskušali izbrati najustreznejšo kopulo. Sliki 3 in 5 prikazujeta dva grafič na testa, kjer so prikazani rezultati za kopulo Frank iz Arhimedove družine kopul. Z rdeč imi znaki x so ponovno označ eni letni maksimumi, s sivo barvo pa so označ eni generirani podatki, ki so bili določ eni na podlagi definiranega kopula modela, ki smo ga določ ili na podlagi podatkov o letnih maksimumih. Z uporabo grafič nih in statistič nih testov smo kot najustreznejšo določ ili kopulo Gumbel-Hougaard iz Arhimedove družine. Slika 4 prikazuje rezultate skupne (angl. joint) povratne dobe OR za kopulo Gumbel- Hougaard iz Arhimedove družine. Izrisani so rezultati za vrednost volumna visokovodnega vala 2,4* 8 m 3 . Pogojna povratna doba OR ponovno označ uje primer, ko se zgodi vsaj ena izmed obravnavanih spremenljivk. 15 Slika 3: Grafič ni test I za kontrolo ustreznosti kopule Frank za podatke z vodomerne postaje Gornja Radgona na reki Muri Slika 4: Skupna povratna doba OR za kopulo Gumbel-Hougaard za vrednost pretoka 2,4*10 8 m 3 na postaji Gornja Radgona na reki Muri 16 Slika 5: Grafič ni test II za kontrolo ustreznosti kopule Frank za podatke z vodomerne postaje Gornja Radgona na reki Muri Ocenjevanje vrednosti koncentracij suspendiranih snovi na podlagi vrednosti pretokov in padavin Kot zadnji praktič ni primer smo obravnavali model, s katerim lahko ocenimo vrednosti koncentracij suspendiranih snovi (SSC) na podlagi vrednosti pretokov in padavin. Ker so meritve koncentracij suspendiranih snovi precej redkejše kot meritve pretokov in padavin, se pogosto pojavi potreba po oceni manjkajoč ih podatkov. Agencija RS za okolje je v letu 2012 (zač asno) prenehala izvajati monitoring te hidrološke spremenljivke (Bezak et al., 2013b). Najpogosteje se za ocenjevanje vrednosti koncentracij suspendiranih snovi uporabljajo različ ne regresijske funkcije, kot so linearna funkcija, polinomi višjih redov ali eksponentna funkcija. Uporaba teh preprostih modelov je smiselna, ko je korelacija med vrednostmi pretokov in koncentracijami relativno visoka. V primeru, ko je korelacija majhna oz. blizu vrednosti 0, to pomeni, da visoke (nizke) vrednosti pretokov ne pomenijo nujno visokih (nizkih) vrednosti koncentracij suspendiranih snovi. V tem primeru je uporaba regresijskih krivulj vprašljiva. V predstavljenem primeru smo za ocenjevanje manjkajoč ih SSC vrednosti uporabili meseč ne vsote treh merjenih spremenljivk (pretoki, koncentracije suspendiranih snovi ter padavine). Možna je uporaba modela tudi z upoštevanjem dnevnih serij podatkov, vendar je avtokorelacija v tem primeru izrazitejša, zato je potrebno za opis posameznih robnih porazdelitev uporabiti t.i. č asovne modele (angl. time series model), kot je npr. ARMA model (angl. autoregressive moving average model). Z uporabo teh modelov v bistvu transformiramo izhodišč ne dnevne vrednosti 17 spremenljivk tako, da se znebimo avtokorelacije v podatkih. Pogoj za izvedbo vseh analiz s funkcijo kopula ter tudi za izvedbo verjetnostnih analiz je, da podatki niso avtokorelirani. To pomeni, da vrednost pretoka v nekem dnevu ni odvisna od vrednosti pretokov v prejšnjih dnevnih. Zato smo se odloč ili, da bomo model za ocenjevanje koncentracij suspendiranih snovi najprej poskusili narediti na meseč nih vsotah. Tak model bi npr. lahko uporabili za ocenjevanje masne bilance lebdeč ih plavin. Uporaba vrednosti padavin je smiselna, ker smo uporabili podatke z vodomerne postaje Ranca na reki Pesnici, ki ima prispevno površino veliko 84 km 2 , in ker padavinska postaja Polič ki vrh leži blizu hidrološke postaje, poleg tega pa je korelacija med meseč nimi vsotami padavin in meseč nimi vsotami koncentracij suspendiranih snovi več ja kot korelacija med pretoki ter koncentracijami. Za izvedbo analize smo uporabili 45 meseč nih vsot pretokov, koncentracij suspendiranih snovi ter meseč nih vsot padavin (preglednica 1). Najprej smo preverili, ali je v vzorcu prisotna avtokorelacija. Z uporabo Box-Pierce testa smo ugotovili, da je v seriji koncentracij suspendiranih snovi prisotna avtokorelacija, smo te podatke najprej transformirali z uporabo Box-Cox transformacije. Za preverjanje stacionarnosti in homogenosti vzorca smo uporabili še Mann-Kendall ter SNHT (angl. standard normal homogeneity test) testa. Ugotovili smo, da so vsi vzorci primerni za izvedbo analiz. Z uporabo grafič nih in statistič nih testov smo kot robne porazdelitve določ ili naslednje funkcije: generalizirana porazdelitev ekstremnih vrednosti za pretoke, Gumbelova porazdelitev za padavine ter log-normalno porazdelitev za koncentracije suspendiranih snovi. Za ocenjevanje parametrov smo ponovno uporabili metodo momentov L (Hosking & Wallis, 1997). Z uporabo statistič nega testa CramØr-von Mises in grafič nih testov smo kopulo Gumbel-Hougaard določ ili kot najprimernejšo za izvedbo analize. Nato smo generirali 10000 vzorcev na podlagi enakomerne porazdelitve. Z Newtonovo metodo smo poiskali nič le funkcije, zapisane z enač bo 4. Slika 6: Primer raztrosa možnih rešitev enač be 4 v obliki histograma z izrisano neparametrič no gostoto verjetnosti Slika 6 prikazuje raztros možnih rešitev enač be 4, ki smo jih dobili s simulacijo, na sliki pa je izrisana tudi neparametrič na gostota verjetnosti, na podlagi katere smo določ ili 18 najverjetnejšo vrednost porazdelitvene funkcije (maksimum neparametrič ne gostote verjetnosti). Za vsak mesec smo tako dobili tak raztros vrednosti (slika 6), na podlagi teh rezultatov pa smo določ ili tudi empirič ne intervale zaupanja. Slika 7: Primerjava rezultatov kopula modela (levo) za ocenjevanje vrednosti koncentracij suspendiranih snovi z izrisanimi 20 % intervali zaupanja in rezultati obič ajnih regresijskih modelov (desno) Primerjava med modelom kopula ter med obič ajnimi regresijskimi modeli je prikazana na sliki 7. Pri modelu kopula so izrisani tudi empirič ni intervali zaupanja. Nato smo z uporabo testov primerjali 4 predlagane modele. Odloč ili smo se za uporabo Nash-Sutcliffe, RMSE (angl. root mean square error) ter MAE (angl. mean absolute error) testov, poleg tega smo izrač unali tudi povpreč no razliko med izmerjenimi in napovedanimi (modeliranimi) vrednostmi, preverili pa smo tudi, kakšne so razlike med minimalnimi, maksimalnimi in srednjimi vrednostmi. Rezultati so prikazani v preglednici 4. Opazimo lahko, da pri modelu kopula dobimo boljše rezultate kot pri regresijskih modelih. Rezultati tudi z uporabo modela kopula niso najboljši (slika 7), so pa boljši kot tisti, do katerih bi prišli z uporabo metod, ki se več inoma uporabljajo za ocenjevanje vrednosti koncentracij suspendiranih snovi. Iz preglednice 4 lahko vidimo, da regresijski modeli podcenijo maksimalne vrednosti in precenijo minimalne vrednosti koncentracij. Tudi v tem primeru smo z uporabo modela kopula dosegli izboljšanje rezultatov. Cilj je dobiti č im manjše vrednosti testov RMSE in MAE, ter vrednost testa NSE č im bližje vrednoti 1. Regresijski modeli so dali boljše rezultate le pri oceni povpreč ne vrednosti, kar pa je verjetno posledica dejstva, da je varianca teh modelov manjša. 19 Preglednica 4: Primerjava rezultatov različ nih modelov za ocenjevanje vrednosti koncentracij suspendiranih snovi Podatki Model kopula Linearni model Eksponentni model Polinom 2. stopnje Nash-Sutcliffe 1 0,22 0,16 0,08 0,16 RMSE 0 1015 1053 1101 1048 MAE 0 682 773 798 762 Povpreč na razlika [%] 0 6,7 28,6 29,5 28,6 Minimalna napoved [g/m 3 ] 604 918 1243 997 1364 Srednja napoved [g/m 3 ] 1724 1520 1725 1712 1725 Maksimalna napoved [g/m 3 ] 4939 4058 3037 2421 3360 Analiza rezultatov in razprava Pri bivariatnih verjetnostnih analizah visokovodnih valov smo ugotovili, da razlike med posameznimi kopulami iz Arhimedove družine kopul, kot so Gumbel-Hougaard ali Frank, niso izrazite (Šraj et al., 2013a). Za modeliranje para konice pretokov in pripadajoč ih volumnov visokovodnih valov smo izbrali kopulo Gumbel-Hougaard, kar je v skladu z ugotovitvami nekaterih drugih raziskovalcev (Poulin et al., 2007; Karmakar & Simonovic, 2009). Uporaba kopul, ki nimajo poudarka na zgornjem robu porazdelitve (angl. upper tail dependence), kar je recimo kopula Frank iz Arhimedove družine kopul, lahko pripelje do tega, da podcenimo vrednosti, ki pripadajo določ eni povratni dobi (Poulin et al., 2007). Pri primerjavi univariatnih (obič ajnih) in skupnih bivaraitnih povratnih dob smo ugotovili (Šraj et al., 2013a), da velja zveza , < -. < , , kar je v skladu z ugotovitvami drugih raziskovalcev (Salvadori et al., 2007). Pri trivariatnih verjetnostnih analizah smo uporabili različ ne v literaturi pogosto uporabljene kopule (Ma et al., 2011; Ganguli & Reddy, 2013). Z uporabo statistič nih ter v literaturi pogosto uporabljenih grafič nih testov (Genest & Favre, 2007) smo določ ili najustreznejšo funkcijo kopula. Ponovno smo izbrali kopulo Gumbel-Hougaard. Izrač unali pa smo tudi skupni povratni dobi OR ter AND (Ganguli & Reddy, 2013). Tudi pri trivariatnem primeru smo ugotovili, da velja zveza ,, < -. < ,, . Vrednosti povratnih dob so bile odvisne od lastnosti funkcij kopul, predvsem od obnašanja na zgornjem robu porazdelitve (Poulin et al., 2007). Pri trivariatnih verjetnostnih analizah konic pretokov, volumnov valov ter vrednosti koncentracij suspendiranih snovi smo uporabili enak postopek kot pri trivariatnih analizah visokovodnih valov. Z uporabo grafič nih in statistič nih testov (Genest & Favre, 2007; Genest et al., 2009; Kojadinovic et al., 2011) smo določ ili kopulo Gumbel-Hougaard kot najprimernejšo za izvedbo analiz. To kopulo smo izbrali tudi zaradi dejstva, ki smo ga upoštevali že pri bivariatnih analizah, in sicer da lahko nekatere kopule podcenijo ocenjene vrednosti spremenljivk, ki pripadajo izbrani skupni povratni dobi (Poulin et al., 2007). Predstavljena metodologija se je izkazala za uporabno, predvsem pa bi rezultati verjetnostnih analiz morali biti bolj zanesljivi, saj v analizah hkrati upoštevamo tri spremenljivke, medtem ko pri obič ajnih univariatnih verjetnostnih analizah upoštevamo le eno spremenljivko. 20 Ocenjene vrednosti koncentracij suspendiranih snovi, ki so rezultat kopula modela, bi lahko označ ili za bolj zanesljive od tistih, ki jih lahko izrač unamo z regresijskimi modeli (slika 7), saj poleg vrednosti pretokov pri ocenjevanju upoštevamo tudi vrednosti padavin, kar predstavlja izboljšavo do sedaj v praksi največ krat uporabljenih metod. Postopek izvedbe analize je enak kot pri prvih treh problemih, ki smo jih obravnavali, s to razliko, da nas tukaj ne zanimajo povratne dobe, ampak želimo izvesti simulacije z uporabo enač be 4. Teoretič no bi lahko v model vpeljali še več spremenljivk in tako poskušali izboljšati konč ne rezultate modela, vendar bi v tem primeru postopek vseboval dodatne korake (Salvadori et al., 2007). Število uporabljenih parametrov pri kopula modelu je nekoliko več je od števila parametrov pri regresijskih modelih, saj moramo oceniti parameter kopule, ki opisuje odvisnost (vse uporabljene kopule imajo 1 parameter, obstajajo pa tudi kompleksnejše kopule z več parametri), poleg tega pa ocenjene vrednosti določ a tudi pasovna širina (angl. bandwidth) izbrane neparametrič ne porazdelitvene funkcije. Določ iti pa je potrebno še parametre porazdelitvenih funkcij, ki opisujejo modelirane spremenljivke. Zaključ ki V prispevku smo prikazali uporabo funkcij kopula na štirih praktič nih primerih. Obravnavali smo bivariatne ter trivariatne analize visokovodnih valov reke Save v Litiji, naredili smo trivariatne verjetnostne analize konic pretokov, volumnov valov ter koncentracij suspendiranih snovi, za kar smo uporabili podatke s postaje Gornja Radgona na reki Muri. Poleg tega smo definirali model, s katerim lahko na podlagi vrednosti pretokov ter padavin ocenimo vrednosti lebdeč ih plavin (SSC) in obenem tudi definiramo empirič ne intervale zaupanja. Za izgradnjo tega modela smo sicer uporabili meseč ne vsote hidroloških spremenljivk, vendar bi bila aplikacija možna tudi na dnevnih vrednostih, kar pa presega okvire tega prispevka. Med izvajanjem analiz smo ugotovili, da je za hidrološke probleme v več ini primerov najprimernejša kopula Gumbel-Hougaard, ki spada tako v Arhimedovo družino kot v družino kopul ekstremnih vrednosti. Za ocenjevanje parametrov kopul smo uporabili različ ne metode, kot so metoda momentov, metoda največ jega verjetja ter pseudo metoda največ jega verjetja. Pri verjetnostnih analizah smo izrač unali skupne povratne dobe imenovane OR ter AND, pri ocenjevanju parametrov pa smo s simulacijami ocenili vrednosti koncentracij suspendiranih snovi. Ker smo model za ocenjevanje SSC vrednosti definirali sami, rezultatov analiz nismo mogli primerjati z rezultati iz literature, ampak le z rezultati regresijskih modelov, kjer smo ugotovili, da je predstavljena metodologija za ocenjevanje vrednosti koncentracij suspendiranih snovi primerna. Možna pa bi bila tudi razširitev modela, tako da bi upoštevali dnevne vrednosti spremenljivk. Uporaba funkcij kopula na praktič nih primerih je pokazala, da so kopule uporabno matematič no orodje, ki ga lahko uporabimo pri številnih praktič nih problemih v hidrologiji, kjer nastopa več medsebojno bolj ali manj odvisnih spremenljivk. Zahvala Zahvaljujemo se Agenciji RS za okolje za posredovane podatke, ki smo jih uporabili za izvedbo analiz. Del rezultatov raziskave je nastal v okviru temeljnega raziskovalnega projekta J2-4096, ki ga financira Javna agencija za raziskovalno dejavnost Republike Slovenije. Del rezultatov pa je prispevek dela UL FGG na mednarodnem raziskovalnem 21 projektu SedAlp, ki ga financira Evropska unija v okviru programa Alpine Space. Zahvaljujemo se tudi obema recenzentoma za njune koristne komentarje. Literatura ARSO. 2013. Podatki, Hidrološki arhiv, Arhiv površinskih voda. Dostopno na: http://vode.arso.gov.si/hidarhiv/pov_arhiv_tab.php. Bardossy, A., Li, J. (2008). Geostatistical interpolation using copulas, Water resources research 44(7). Bezak, N., Brilly, M., Šraj, M. (2013a). Comparison between the peaks over threshold method and the annual maximum method for flood frequency analyses, Hydrological Science Journal, doi: 10.1080/02626667.2013.831174. Bezak, N., Šraj, M., Mikoš, M. (2013b). Pregled meritev vsebnosti suspendiranega materiala v Sloveniji in primer analize podatkov, Gradbeni vestnik (v tisku). De Michele, C., Salvadori, G., Canossi, M., Petaccia, A., Rosso, R. (2005). Bivariate Statistical Approach to Check Adequacy of Dam Spillway, Journal of Hydrological Engineering 10(1), 50- 57. Favre, A. C., Aldouni, E.L., Perreault L., Thiemonge, N., Bobee, B. (2004). Multivariate hydrological frequency analysis using copulas, Water resources research 40(1). Ganguli, P., Reddy, M. J. (2013). Probabilistic assessment of flood risks using trivariate copulas, Theoretical and Applied Climatology 111(1-2), 341-360. Genest, C., Favre, A. C. (2007). Everything you always wanted to know about copula modelling but were afraid to ask, Journal of Hydrological Engineering 12, 347-368. Genest, C., Remillard, B., Beaudoin, D. (2009). Goodness-of-fit tests for copulas: A review and a power study, Insurance: Mathematics and Economics 44, 199-213. Grimaldi, S., Serinaldi, F. (2006). Asymmetric copula in multivariate flood frequency analysis, Advances in Water Resources 29(8), 1155-1167. Graler, B., van den Berg, M. J., Vandenberghe, S., Petroselli, A., Grimaldi, S., De Baets, B., Verhoest, N. E. C. (2013). Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation, Hydrology and Earth System Sciences 17(4), 1281-1296. Hofert, M., Mä chler, M. (2011). Nested Archimedean Copulas Meet R: The nacopula Package, Journal of Statistical Software 39(9). Hosking, J. R. M., Wallis, J. R. (1997). Regional frequency analysis: an approach based on L- moments. Cambridge, Cambridge University Press, 224 p. Joe, H. (1997). Multivariate Models and Multivariate Dependence Concepts. Chapman & Hall/CRC Monographs on Statistics & Applied Probability, 424 p. Karmakar, S., Simonovic, S. P. (2009). Bivariate flood frequency analysis. Part 2: a copula based approach with mixed marginal distributions, Journal of Flood Risk Management 2, 32-44. Kojadinovic, I., Yan, J., Holmes, M. (2011). Fast large-sample goodness-of-fit for copulas, Statistica Sinica 21(2), 841-871. Ma, M., Song, S., Ren, L., Jiang, S., Song, J. (2011). Multivariate drought characteristics using trivariate Gaussian and Student t copulas, Hydrological Processes 27(8), 1175-1190. Nelsen, R. B. (2006). An Introduction to Copulas. Springer Series in Statistics, 272 p. Poulin, A., Huard, D., Favre, A. C., Pugin, S. (2007). Importance of tail dependence in bivariate frequency analysis, Journal of Hydrological Engineering 12, 394-403. Salvadori, G., De Michele, C. (2004). Frequency analysis via copulas: Theoretical aspects and applications to hydrological events, Water resources research 40(12). Salvadori, G., De Michele, C., Durante, F. (2011). On the return period and design in a multivariate framework, Hydrology and Earth System Sciences 15, 3293-3305. 22 Salvadori, G., De Michele, C., Kottegoda, N., Rosso, R. (2007). Estremes in nature: An Approach Using Copulas. Water Science and Technology Library, 292 p. Shiau, J. T., Feng, S., Nadaraiah, S. (2007). Assessment of hydrological droughts for the Yellow River, China, using copulas, Hydrological Processes 21(16), 2157-2163. Sklar, A. (1959). Fonctions de repartition a n dimensions et leurs marges, Publ. Inst. Statist. Univ. Paris 8, 229-231. Šraj, M., Bezak, N. (2013). Analiza visokovodnih valov Save v Litiji, Ujma 27, 228-235. Šraj, M., Bezak, N., Brilly, M. (2013a). Bivariate flood frequency analyses using copula function, European Geosciences Union 2013. Šraj, M., Bezak, N., Brilly, M. (2013b). Vpliv izbire metode na rezultate verjetnostnih analiz konic, volumnov in trajanj visokovodnih valov Save v Litiji, Acta Hydrotechnica (v tisku). Wong, G., Lambert, M. F., Leonard, M., Metcalfe, A. V. (2010). Drought Analysis Using Trivariate Copulas Conditional on Climate States, Journal of Hydrological Engineering 15, 129- 141. Zhang, L., Singh, V. P. (2007). Bivariate rainfall frequency distributions using Archimedean copulas, Journal of Hydrology 332(1-2), 93-109.