ZÁKLADNÍ ODTOK A JEHO SOUVISLOSTI S VYHODNOCOVÁNÍM MĚRNÝCH KŘIVEK PRŮTOKŮ Ladislav Budík Český hydrometeorologický ústav, pobočka Brno Kroftova 43, 616 67 Brno
[email protected]
Abstrakt: Příspěvek popisuje nově navrženou metodu sestavení měrné křivky průtoků založenou na vyjádření vztahu mezi stavem a průtokem pomocí jediné funkce jedné nebo více proměnných. Nastiňuje jiný způsob interpolace hydrogramu, než lineární, pomocí řešení Bernoulliovy rovnice. Zkoumá také dopad takovéto změny způsobu převodu stavu na průtok na možnosti výpočtu tzv. základního odtoku. Klíčová slova: základní odtok, hydrogram, výtoková křivka, měrná křivka průtoků, Bernoulliova rovnice 1. Úvod Hydrologové měří stavy hladiny na tocích nepřetržitě a několikrát za rok také měří průtok. Z těchto několika měření získají vztah mezi stavem a průtokem. - měrnou křivku průtoků (MKP). Ta se nejčastěji aproximuje pomocí splajnů.. Takto získaná MKP však není hladká a obtížně se s ní pracuje. Interpolace hydrogramu (průběh průtoku v čase) v místech, kde je to nutné, se provádí lineárně ve stavech, což plně neodpovídá realitě. Tím se do výtokových křivek (poklesové části křivek) a hydrogramu vůbec vnáší nespojitosti, což prakticky znemožňuje použít matematické postupy pro hydrogeologicky zdůvodněný výpočet základního odtoku (ZO). Za základní odtok se považuje odtok z podzemních vod do povrchových toků. V podstatě ho zaznamenáváme v době sucha, kdy delší dobu neprší a v tocích již teče pouze voda vyvěrající z podzemí. Změřit ho nelze, protože je obvykle smíchán s odtokem povrchovým a podpovrchovým.Všechny dosud používané metody pro stanovení ZO (existuje jich asi dvacet) jsou velmi hrubými odhady reality a jejich výsledky se navzájem dosti liší. V ČR je momentálně nejpoužívanější metoda Kliner-Kněžek [1]. 2. Odhad základního odtoku 2.1. Metoda klouzavých minim V publikaci Budík-Budíková [2] je uveden způsob, jak zlepšit odhady ZO získané Castanyho metodou (ta je popsána v [1]). Castanyho metoda je založena na tom, že za ZO je považován nejnižší měsíční průtok. Průběh ZO je pak znázorněn skokovou funkcí času. Zavedl jsem minima počítaná v klouzavém okně délky n:
Yi = min (Q i , K, Q n +i −1 ) , i = 1, …, m – n + 1
(1)
kde Yi lze považovat za první aproximaci odhadu základního odtoku, Q1, …, Qm jsou průtoky (v m3s-1). Tímto způsobem dostaneme skokovou nepravidelnou funkci, jejíž průběh pro stanici Popov na Vláře je uveden na obr.1
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
10
9
8
7
6 2
3
1
0
.
/
,
-
+
(
)
denní průměry průtoků
5
klouzavá minima, okno 31
*
'
&
4
3
2
1
0 5000
5100
5200
5300
5400
5500
!
"
#
y
5600
5700
5800
5900
6000
%
Obr. 1: Průběhy průtoků klouzavých 31 denních minim ve stanici Popov na Vláře Získaná skoková funkce je pak vyhlazena geometrickým průměrem podle vzorce
ZOj = h Yj ⋅K⋅ Yh + j−1 , j = 1, …, m – (h+n) + 2
(2)
Kde ZO je aproximace základního odtoku, n je délka okna, Yi … Ym-n jsou minima, z nichž počítáme ZO. Průběhy odhadu základních odtoků v Popově vidíme na obr. 2:
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
y
e
10
9
8
7
6 3
4
2
1
/
0
-
.
,
)
*
denní průměry průtoků
5
ZO zvyrovnaných klouzavých minim, okno 31
+
(
'
4
3
2
1
0 5000
5100
5200
5300
5400
5500
!
"
#
$
5600 y
5700
5800
5900
6000
&
Obr. 2: Průběhy odhadu základních odtoků ve stanici Popov na Vláře Vyhlazovací okno může být symetrické i asymetrické. Jedná se v podstatě o jednoduché jádrové vyhlazení. Jako nejvhodnější pro odhad ZO se ukázalo okno symetrické shodné délky pro oba navazující výpočty (klouzavá minima i geometrický průměr) s délkou 31 nebo 51 dnů. Při tomto postupu nemůže dojít k tomu, že by vypočtený odhad ZO byl vyšší než naměřený průtok. Takto volená okna dávají dlouhodobé průměry ZO podobné výsledkům některé z dosud užívaných metod (31 dnů Kliner Kněžek, 51 dnů metoda výtokových křivek). V detailech průběhů se od zmíněných dvou liší. Upozorňuji, že z výsledku vůbec nevyplývá, která z metod více odpovídá realitě. Oproti ostatním metodám má tato metoda vyhlazených klouzavých minim značnou výhodu v rychlosti a opakovatelnosti výpočtu. Reálnost výsledku ale neřeší. 2.2. Metoda výtokových křivek Metoda výtokových křivek je založena na fyzikální podstatě problému. Je ovšem zapotřebí rozhodnout, co ještě výtoková křivka je a co už ne a jaký má tvar. Tvar výtokové křivky se předpokládá v exponenciálním tvaru: Q(t ) = e − k⋅t
(3)
kde Q(t) je průtok v čase t, k je reálný parametr. Jelikož z tvaru výtokové křivky by mělo jít odvodit ZO, vyšel jsem z této skutečnosti a pokusil se najít algoritmus a následně napsat počítačový program tak, aby tato metoda odhadu ZO byla co nejobjektivnější a výsledek byl reprodukovatelný. Především je nutno rozlišit ty části hydrogramu, které představují výtokové křivky. Musí to být část poklesová – tedy první derivace (podle času) té části hydrogramu, která může být výtokovou křivkou, musí být záporná. To ovšem může splnit i ta část křivky, kde trochu zaprší a pokles se jen zpomalí nebo naopak zrychlí. Tedy druhá derivace musí být kladná – první derivace musí stabilně růst (stává se méně zápornou). Abychom zcela rozlišili, kdy se pokles zpomalí prudčeji tím, že trochu zaprší, musí třetí derivace být záporná, tj. druhá deri8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
vace se musí stále zmenšovat. Data máme ovšem nespojitá, a to buď se vzorkovacím intervalem jeden den nebo jedna hodina. Místo derivací musíme tedy použít diference příslušného stupně. V první fázi jsem použil denní průměry. Charakteristiky většího množství vybraných výtokových křivek vyneseme do tečkových diagramů, kde na vodorovné ose je ln(dQ/dt) a na svislé ose je lnQ (Q je průtok v čase t) viz obr 3.
e
e
e
4
2 y = 1.05x - 1.3 R2 = 1
0
-8
-6
-4
0
-2
2
4
-2
-4
-6
tečkový diagram
!
"
#
hranice lnQ při daném ln(dQ/dt)
Lineární (hranice lnQ při daném ln(dQ/dt))
Obr. 3: Dvourozměrný tečkový diagram závislosti logaritmu průtoků na 1. diferenci logaritmu průtoků – denní hodnoty, využití prvních tří diferencí Tento tečkový diagram budu dále nazývat stopou říčního profilu (stručně stopou profilu). Název stopa proto, že je pro každý profil zcela jedinečný. Odráží se v něm hydrogeologická skladba povodí a jeho klimatické a vegetační podmínky, jež všechny dohromady utvářejí výtokové křivky. Za přesně stejných podmínek by měla být výtoková křivka vždy stejná, tedy stejná by měla potom být i stopa. V tomto tečkovém diagramu by měly být jako přímky zobrazeny jednotlivé vypreparované výtokové křivky. Vidět to lze jen stěží a s velkou fantazií. V podstatě se ukazuje, že diference denních průměrů vykazují velké skoky oběma směry a že třetí diferenci pro rozlišení už vůbec nejde použít. Nezbudou nám skoro žádné výtokové křivky, což neodpovídá teoretickým předpokladům. Proto sestrojíme tečkové diagramy jen s rozlišením podle 1. a 2. diference. S tímto zjednodušením se musíme pro další analýzy spokojit. Výsledek je znázorněn na obrázku 4:
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
e
e
e
4 y = -0.0034x 3 + 0.0078x 2 + 0.6018x + 0.6191 R2 = 0.6536 2 y = 1.05x - 1.3 R2 = 1
0
-8
-6
-4
-2
0
2
4
-2
-4
-6
tečkový diagram
!
"
#
Lineární (hranice lnQ při daném ln(dQ/dt))
Polynomický (tečkový diagram )
hranice lnQ při daném ln(dQ/dt)
Obr. 4: Dvourozměrný tečkový diagram závislosti logaritmu průtoků na 1. diferenci logaritmu průtoků – denní hodnoty, využití prvních dvou diferencí Stejný postup jsem použil pro vyhodnocení hodinových průtoků - viz obr. 5:
$
%
&
$
'
(
'
4
2
0
-9
-7
-5
-3
1
-1
3
5
data
směrnice čar -2
-4
-6
!
"
#
Obr. 5: Dvourozměrný tečkový diagram závislosti logaritmu průtoků na 1. diferenci logaritmu průtoků – hodinové hodnoty, využití prvních dvou diferencí
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
V tomto tečkovém diagramu lze jednotlivé výtokové křivky vidět bez velké fantazie, ale opět použití 3. diference pro důsledné rozlišení výtokových křivek je téměř nemožné, jak je patrno z obr. 6:
$
%
&
$
'
(
'
4
2
0
-9
-7
-5
-3
1
-1
3
5
data
směrnice čar -2
-4
-6
!
"
#
Obr. 6: Dvourozměrný tečkový diagram závislosti logaritmu průtoků na 1. diferenci logaritmu průtoků – hodinové hodnoty, využití prvních tří diferencí Na obr. 5 jsou také vidět dlouhé soubory hodinových průtoků, kde průtok klesá lineárně se stále stejnou diferencí. K tomu dochází v místech, kde se provádí z důvodu obtížné rozlišitelnosti charakteru poklesu lineární interpolace stavu (dříve i průtoku). Zjistil jsem, že do vyhodnocení se různě ručně zasahuje pro nemožnost příslušné procesy matematicky vyjádřit. Jak již bylo uvedeno, MKP se získává pomocí splajnů 2. řádu, což vlastně řád diferencí také rozbíjí. MKP kolísá v ostrých vlnkách. Tuto určitou kritiku zpracování průtoků, není třeba brát příliš vážně. Popsaný postup pro běžná zpracování průtoků v povrchových vodách je naprosto vyhovující a zatím nikomu nevadil. Pouze znemožňuje spolehlivě matematicky separovat výtokové křivky. To vadí spíše těm, kdo zpracovávají podzemní vody. Závěr tedy je, že hlavní operace, které poškozují výtokové křivky, jsou tři. • Nedefinovatelné a později neopakovatelné ruční zásahy do převodu stavů na průtoky (některé jsou nutné, např. změny MKP, tvoření ledových bariér apod.). Musí být matematicky definovány nebo alespoň zaznamenány. • MKP nejsou zapsány ve tvaru funkce, ale jsou vyjádřeny jako splajny, tudíž mění průběh derivací (diferencí). • Lineární interpolace stavů (ani stav ani průtok se obecně při použití konstantního časového vzorkování nemění lineárně.
3. Návrh zpřesnění postupů při zpracování hydrologických dat 3.1. Využití tečkových diagramů závislosti logaritmu průtoků na 1. diferenci logaritmu průtoků Vraťme se nyní zpět k dvourozměrným tečkovým diagramům závislosti logaritmu průtoků na 1. diferenci logaritmu průtoků. Zkusme proložit těmito tečkovými diagramy metodou 8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
nejmenších čtverců regresní přímky. Zjistíme, že mají většinou index determinace kolem 0,5, ale jejich směrnice je dobře vyjádřena, protože tečkový diagram sestává fakticky z téměř rovnoběžných čar jednotlivých výtokových křivek. To zjistíme poté, kdy sestrojíme tečkové diagramy pro jednotlivé výtokové křivky jedné stanice. Pro názornost uvádím na obr. 7 výsledek nového zpracování hodinových průtoků v roce 2003 pro profil Popov na Vláře, kde na rozdíl od zpracování klasického čáry výtokových křivek vidět jsou.
4
2
0
#
!
"
-12
-10
-8
-6
-4
-2
0
2
4
data okraj podzemního odtokur
-2
-4
-6
Obr. 7: Stopa profilu Popov na Vláře, hodinový krok, nový způsob Obecně každý říční profil má svůj tečkový diagram, jehož tvar a směrnice jsou závislé na hydrogeologických podmínkách území (podzemních i povrchových), které odvodňuje. Kdybychom vynesli směrnice přímek prokládaných tečkovými diagramy do grafu v závislosti na průměrném sklonu povodí, řadily by se do přímky, kde by velikost směrnice rostla s velikostí průměrného sklonu povodí. Se sklonem povodí narůstá logicky rychlost výtoku vody. Tohle platí pouze u toků z území s podkladem karpatského flyše. U toků z ostatního území je nedostatečná homogenita hydrogeologických poměrů. Nyní napíšeme rovnici přímky v tečkovém diagramu: dQ ln (Q ) = a ⋅ ln +b dt
(4)
Je to diferenciální rovnice Bernoulliova a její analytické řešení je známo. Uvedeme řešení pouze pro b > 0:
a b − a ⋅ ln ( t + C ) ⋅ (a − 1) Q = exp a −1
(5)
kde Q je průtok, a je směrnice přímky ve stopě profilu, b je úsek na ose Y vyťatý přímkou o směrnici a, t je čas.
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
Tento vztah popisuje tvar výtokových křivek povodí. Není to tedy exponenciála (i když tvar výtokové křivky vcelku dobře aproximuje), ale výtoková křivka má tvar daný řešením Bernoulliovy rovnice. Dá se říci, že je to v různé míře kombinace exponenciály a hyperboly. To je velmi důležité pro další řešení problému ZO a vlastně kvality průtokových hydrogramů vůbec.
3.2. Měrná křivka průtoků jako funkce více proměnných Kvalita MKP je naprosto zásadní pro převod stavu na průtoky V současné době se začínají prosazovat ultrazvuková měření průtoků. Ta mají střední chybu 1 – 2 % oproti současným měřením vrtulí (střední chyba okolo 5%). Stavy lze již měřit v podstatě kontinuálně (zápis po 10 minutách). Podmínky pro to, abychom znali nejen průtok, ale ve vcelku slušné kvalitě i tvar křivek, a tudíž bylo možno z tvaru hydrogramu vyhodnocovat i jiné chování povodí než jen průtok, se zlepšují. Pro matematická vyhodnocení výtokových křivek je nutné, aby MKP byla spojitá ve všech bodech a aby byla potom spojitá i křivka hydrogramu. Chceme docílit, aby se falešné signály z MKP nepřenášely do křivky hydrogramu. Na základě analogie s řešením Bernoulliovy rovnice jsem navrhl tvar regresní funkce pro stanovení MKP jako modifikovanou exponenciální funkci s možností lineárního posunu celé křivky: sign (A ) ⋅ A
Q(H ) =
⋅s + u
d⋅H + r
(6)
1 je logistická funkce, H je stav (v cm), Q je průtok 1 + i ⋅ j k ⋅H (v m3s-1), d, r, s, u, c, a, b, i, j, k jsou reálné parametry, které získáme iteračním postupem jako regresní koeficienty metodou nejmenších čtverců (MNČ). Parametry s a u mohou celou MKP posouvat a naklánět. Jejich konkrétní využití bude předvedeno později. Na obr. 8 je znázorněn tvar MKP ve stanici Ptáčov na Jihlavě získané popsanou metodou:
kde A = c + e H
a
+b
, funkce F(H ) =
510 460 410 360
"
310
MKP splajny
260
MKPexp
Výsledek
210 160 110 60 10 0.000
50.000
100.000
150.000 "
8.
NÁRODNÍ KONFERENCE
200.000
250.000
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
Obr. 8: MKP ve stanici Ptáčov na Jihlavě Nastanou-li rozlivy v inundaci (tj. voda teče mimo koryto), dochází k lomu MKP, což nelze postihnout funkcí (6). Proto navrhuji extrapolovat základní křivku i do části MKP, kde nevyhovuje funkce (6) a rozdíl od skutečného průběhu MKP pak proložit další funkcí stejného tvaru, jen s jinými koeficienty. Tyto dvě různé křivky se spojí logistickou funkcí podle vzorce
Q(H ) = Q1 (H ) + F1 (H ) ⋅ Q 2 (H )
(7)
kde Q1(H) je funkce (6) s parametry získanými MNČ ze základní části MKP, Q2(H) je funkce (6) s parametry získanými MNČ z rozlivové části MKP a F1(H) je logistická funkce s parametry získanými MNČ tak, aby součet Q1(H) + Q2(H) dal hladký průběh celé MKP. Tímto způsobem lze spojit i více částí MKP. Na obr. 9 je vidět detail rozdílu včetně lomu splajnu. Největší rozdíl mezi náhradou splajnem a funkční MKP je v ohybu na přechodu do rozlivu, a to 2,5%.
150 140 130 120
"
110
100
MKP splajny
Výsledek 90 80 70 60 50 0.000
2.000
4.000
6.000
8.000
10.000 "
12.000
14.000
16.000
18.000
20.000
Obr. 9: Detail MKP v Ptáčově na Jihlavě V případech, kdy se nemění tvar břehů, ale jen tvar dna nebo průtočnost koryta (např. profil zarůstá v létě vegetací), stačí pouze vyhodnotit míru změny průtočnosti koryta (např. přímým měřením průtoku nebo jiným způsobem porovnání). Není třeba stanovovat nové všechny parametry MKP. Stačí odhadnout nové parametry s a u a celá MKP se posune a pootočí do nové polohy, jak je znázorněno na obr. 10:
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
y
e
e
150
140
130
data bez z áros tu funk č ní M K P č ervený m i daty
,
+
)
*
(
&
funk č ní M K P m odrý m i daty
120
s tandardní M K P
'
%
data s podobný m č as em
P os unutá funk č ní M K P (původněč ervený m i body ) 110
data s e z áros ty
100
90 0
0.2
0.4
0.6
0.8
1
1.2
!
"
#
1.4
1.6
1.8
2
$
-
#
.
Obr. 10: MKP ve Velkých Pavlovicích na Trkmance Pro tuto operaci by mělo stačit, že se nemění tvar břehů a máme k dispozici jedno přímé měřění průtoku s odchylným výsledkem než dosud. Parametry s a u lze také stanovit také jako funkce času a MKP bude svou polohu v čase měnit plynule a tudíž nevzniknou žádné skoky, jež by bylo nutno ručně vyrovnávat. Také změna průtoku bude věrnější. Protože lze takto zadanou měrnou křivku měnit plynule v čase, máme k dispozici ještě jeden mocný nástroj - možnost, že MKP bude dvourozměrná. Průtok může záviset na stavu ve dvou různých recipientech. Mějme profil, který měří průtok nad přehradou. Je-li v přehradě příliš mnoho vody, pak dojde ke vzdutí vody v měrném profilu (lze si představit i jiný způsob vzdouvání např. na soutoku dvou toků). Budeme-li kromě stavu ve vlastním měrném profilu znát stav vody v druhém recipientu, je možné konstruovat dvojrozměrnou MKP, ve funkční MKP je vlastně počítána oprava stavu v měrném profilu na vzdutí. Takto konstruovaná MKP bude plynule reagovat na změny stavu v obou recipientech zároveň, do Q a A ve vzorci (6) dosadíme za H výraz (8): H = H'−a ⋅ (H R − H`o ) ⋅ (H`− H`o ) (8) kde H je opravená výška hladiny v měrném profilu, HR je výška hladiny v ovlivňujícím recipientu, H`je stav vody v měrném profilu, a, Ho jsou parametry zjišťované iterací, přičemž Ho je v podstatě hladina v recipientu, při níž již nedochází k ovlivnění hladiny v měrném profilu. Na grafu závislosti průtoku na stavu se pak zřetelně zobrazí hladká hysterezní křivka, kterou klasickou metodou v takové kvalitě nezískáme. Na obr. 11. je znázorněna hysterezní křivka pro použitou dvourozměrnou plynule se měnící MKP. Pro porovnání je na obr. 12 týž výpočet s dvěma klasickými MKP měněnými skokem. V místě přerušení bude navázání provedeno ručně.
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
5 4.5 4 3.5
"
3
2.5 2 1.5 1 0.5 0 0
50
100
150
200
250 "
300
350
400
450
Obr. 11: Hysterezní smyčka na stanici Podhradí na Dyji – funkční MKP
y
e
e
y
y
e
5 4.5 4 3.5
"
3
2.5 2 1.5 1 0.5 0 0
50
100
150
200
250 "
300
350
400
450
Obr. 12: Hysterezní smyčka na stanici Podhradí na Dyji – klasické zpracování (autor Pavel Kupsa, ČHMÚ)
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
Hysterezní smyčka je získaná z povodně 2006 ve stanici Podhradí na Dyji, kterou ovlivňují vysoké stavy nádrže Vranov. Jedna nezávisle proměnná jsou stavy ve stanici Podhradí a druhá jsou stavy nádrže Vranov. Ultrazvuková měření průtoků umožňují rozlišit i velikost průtoku při stejném stavu na vzestupné větvi průtokové vlny a na větvi sestupné. Tento rozdíl vzniká různým sklonem hladiny toku. Na vzestupné větvi průtokové vlny do prázdného koryta vtéká vysoký průtok a na sestupné větvi do plného koryta menší průtok. O této hysterezi se pochopitelně ví, ale dosavadními metodami nebylo možno ji běžně vyhodnocovat. Zaměníme-li ve výše uvedeném vztahu pro dvourozměrnou MKP stav druhého recipientu změnou stavu hladiny v měrném profilu v posledním časovém intervalu měření (ohodnotíme, jak rychle vlna stoupá nebo klesá), získáme možnost vyhodnotit i rozdíl v průtoku mezi nástupnou a poklesovou větví průtokové vlny. Pro momentální nedostatek vhodných dat nelze zatím předvést ukázku.
3.3. Interpolace a vyrovnání výsledných průtoků pro výpočet ZO – interpolace hydrogramů Na měření v profilu Popov na řece Vláře se pokusím ukázat výhody interpolace s využitím řešením Bernoulliovy rovnice. Použiji k tomu hodinové průtoky získané interpolací průtoků v tzv. lomových bodech, které se archivují. Nabízí se řada možností zpracování. a) Lineární interpolace stavů: V tomto případě zbývá pouze převést pomocí MKP stavy na průtoky. Na obr.13 je vidět, jak dopadnou výtokové křivky v grafu s osami ln(dQ/dt) a lnQ. Jsou nespojité a přerušované, křivky jsou lomené.. Pro konstrukci ZO nevyhovují.
e
e
e
10 y = 0.4201x + 2.4547 R2 = 0.9333
5
0 -20
-15
-10
-5
0
5
lineárně interpolovaná data
směr přímek stop -5
Lineární (lineárně interpolovaná data)
-10
-15
Obr. 13: Stopa výtokové větve při lineární interpolaci stavů b) Průtoky s pomocí MKP spočítáme pouze v lomových bodech a interpolaci provedeme až v průtocích, a to lineárně. Výsledek je podobný jako na obr. 13. Pro jednu výtokovou křivku je vidět, že se jedná o lomenou čáru. Pro konstrukci ZO je tento způsob opět nevyhovující.
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
e
e
e
10 y = 0.5719x + 2.7221 R2 = 0.7668
5
0 -20
-15
-10
-5
0
5
interpolace lineární v průtocích
sklon obrazu lineární interpolace
Lineární (interpolace lineární v průtocích)
-5
-10
-15
Obr. 14: Stopa výtokové větve při lineární interpolaci průtoků
c) Průtoky s pomocí MKP spočítáme pouze v lomových bodech a interpolaci provedeme až v průtocích, a to exponenciálně, jak předpokládá dosud platná teorie. Na obr. 15 je vidět, že výtokové křivky vypadají mnohem lépe. Nicméně stále jsou v grafu ln(dQ/dt) lnQ přerušované a nespojité více, než je přijatelné.
e
e
e
e
10 y = 0.4407x + 2.5429 R2 = 0.9514
5
0 -20
-15
-10
-5
0
5
interpolace exponenciální funkcíí
sklon obrazu lineární interpolace
-5
Lineární (interpolace exponenciální funkcíí)
-10
-15
Obr. 15: Stopa výtokové větve při exponenciální interpolaci průtoků
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
d) Průtoky s pomocí MKP spočítáme pouze v lomových bodech a interpolaci provedeme až v průtocích, a to s pomocí řešení Bernoulliovy rovnice. Pro získání koeficientů stačí tři body křivky. Parametry lze získat iterací. Pro praktické použití je ale postup složitý. Na obr. 16 je vidět, že křivka v grafu lnQ ln(dQ/dt) je již více méně plynulá a mohla by být pro matematickou konstrukci ZO již vyhovující.
e
e
e
&
&
'
10 y = 0.4177x + 2.4153 R2 = 0.9622
5
0
-20
-15
-10
0
-5
5
interpolace Bernolliovou rovnicí
sklon obrazu lineární interpolace -5
Lineární (interpolace Bernolliovou rovnicí)
-10
-15
Obr. 16: Stopa výtokové větve při interpolaci průtoků Bernoulliovou rovnicí Na obrázcích 17 a 18 uvádím též pro srovnání vyhodnocení pro celý rok 2003 v hodinovém kroku pro klasický způsob s ručními úpravami a lineárními interpolacemi a nový způsob založený na Bernoulliově rovnici.
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
4
2
0
-9
-7
-5
-3
-1
1
3
5
data
směrnice čar -2
-4
-6
Obr. 17: Stopa profilu Popov na Vláře, hodinový krok, klasický způsob
4
2
0
"
!
-12
-10
-8
-6
-4
-2
0
2
4
data okraj podzemního odtokur
-2
-4
-6
Obr. 18: Stopa profilu Popov na Vláře, hodinový krok, nový způsob Na obr. 18 jsou zřetelně vidět letní křivky s nižším ZO v dolní části diagramu a zimní křivky s vyšším ZO v jeho horní části. Na obr. 17 jsou jejich stopy rozeznatelné při porovnání s obr. 18, při samostatném posouzení však ne. Na obr. 18 je i velmi pěkný zákres jednotlivých výtokových křivek se zřetelně vyjádřeným přechodem ze zimních křivek na letní v pravé části grafu. Přechod se v roce 2003 odehrával v květnu a v první polovině června (jedná se o horské povodí s velkým podílem listnatých lesů).
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
3.4. Odvození podzemního odtoku (ZO) z poklesové větve hydrogramu Podívejme se na tečkový diagram na obr. 19.
e
e
e
4 y = -0.0034x 3 + 0.0078x 2 + 0.6018x + 0.6191 2 R = 0.6536 2 y = 1.05x - 1.3 R2 = 1
0
-8
-6
-4
-2
0
2
4
-2
-4
-6
tečkový diagram
hranice lnQ při daném ln(dQ/dt)
Polynomický (tečkový diagram )
Lineární (hranice lnQ při daném ln(dQ/dt))
Obr. 19: Stopa profilu Vlára, Popov, klasické zpracování, denní krok Na mnoha profilech je i při klasickém zpracování a jednodenním kroku vidět útvar, který připomíná jakousi „hokejku“. Zvláště patrné je to na profilech s velkými zásobami podzemních vod. Spodní část „hokejky“ s malým sklonem zřejmě reprezentuje tu část výtokové křivky, kde do toku již vytékají jen podzemní vody (ZO), šikmá část reprezentuje kombinaci podzemního, podpovrchového a povrchového odtoku. Při jednodenním kroku nelze očekávat dobrou separovatelnost odtoku povrchového od podpovrchového. Lze předpokládat, že podzemní odtok (a tedy i ZO) může dosáhnout pouze určité maximální hodnoty a nemůže se dál zvyšovat, neboť pak už se stává odtokem podpovrchovým, případně, stoupne-li podzemní voda až k povrchu terénu, odtokem povrchovým. Proto po pravé straně tečkového diagramu na obr. 19 můžeme vést omezující šikmou čáru (na obr. 19 je nakreslena žlutě), za kterou by už podzemní (výše podpovrchový) odtok neměl jít. (Stále nesmíme zapomínat, že máme před sebou graf ln(dQ/dt), lnQ popisující Bernoulliovu rovnici.). Vybereme-li bod reprezentující celkový odtok, hodnota příslušejícího podzemního odtoku by měla přibližně ležet pod ním na rovnoběžce se žlutou čarou. Protože ale v tečkovém diagramu jsou pouze celkové změřené průtoky, nemá v něm příslušný podzemní odtok reprezentanta ve formě tečky. Výjimku tvoří pouze případy, kdy celkový odtok je složen pouze z podzemního odtoku v levé části tečkového diagramu.. Abychom mohli dále s tečkovým diagramem pracovat, musíme ho transformovat tak, aby šikmá žlutá omezující čára z obr. 19 byla svislá a tudíž sobě odpovídající celkový a podzemní odtok ležely na svislici nad sebou a tak je bylo možno snadno přiřadit k sobě. Tečkový diagram po této transformaci je na obr.20.
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
e
y
e
y = 1.1159x - 1.024 R2 = 0.1527
8 3
2
y = -0.0619x - 0.1942x + 1.4725x + 0.0402 R2 = 0.3538
transformovaný tečkový diagram
e
průměry skupin po intervalech osy x
3
posunutá charakteristická křivka pro výpočet ZO
1
Pravděpodobné čisté podzemní odtoky
+
e
-5
-4
-3
-2
-1
0
1
+
( .
-2
pravděpodobná nejvyšší hranice podzemního odtoku při daném ln(dQ/dt)
e
1 +
y = -0.04x 3 - 0.104x 2 + 1.5569x - 0.0004 R2 = 0.921
Polynomický (transformovaný tečkový diagram)
e
(
*
Polynomický (průměry skupin po intervalech osy x) -7
Lineární (Pravděpodobné čisté podzemní odtoky)
Lineární (pravděpodobná nejvyšší hranice podzemního odtoku při daném ln(dQ/dt))
-12
Obr.20: Transformovaná stopa profilu Vlára, Popov, klasické zpracování, denní krok
I podzemní odtok se řídí Bernoulliovou rovnicí, ovšem s jinými parametry než povrchový a podpovrchový odtok. Oblast ve stopě profilu, kde se vyskytují s vysokou pravděpodobností pouze tečky podzemního odtoku, jsou označena bleděmodře a také textem „pásmo podzemního odtoku“. Dále je označeno textem pásmo, kde se vyskytuje pouze povrchový a podpovrchový odtok a pod ním pásmo, kde je směs všech tří druhů odtoku. Nyní stačí v tomto transformovaném grafu pouze určit, jaká hodnota podzemního odtoku pravděpodobně patří k dané hodnotě celkového odtoku. Z horního omezení tohoto bleděmodrého pásma podzemního odtoku vychází směrem vpravo extrapolace čáry horního omezení možné velikosti podzemního odtoku. Nad touto čarou je již podzemní odtok nasycen a nemůže se dál zvyšovat, protože přejde v jiné formy odtoku. K určení hodnoty podzemního odtoku (ZO) bychom měli využít jednak celkovou hodnotu odtoku a jednak rychlost jejího poklesu (1. derivaci, v praxi 1. diferenci). K horním okrajovým hodnotám celkového odtoku patří na svislici horní okrajové hodnoty podzemního odtoku a k dolním celkovým zase dolní podzemní. V levé části grafu se celkové hodnoty průtoku již kryjí s hodnotami podzemních průtoků (označeny bleděmodře). U hodnot celkového odtoku by se měla odlišovat velikost směrnice poklesu výtokové křivky od směrnice poklesu výtokové křivky čistého podzemního odtoku (ta je nižší). Na tuto odlišnost by měla být zavedena redukce (u prudkých lijáků může být vysoký celkový odtok s prudkým poklesem, protože se jedná převážně o odtok povrchový, a malý odtok podzemní). To 24 hodinové průměry ani hodinové při standardním zpracování neumožňují provést. Diference silně kolísají kvůli drobným nepřesnostem převodu stavu na průtok. Proto výsledný ZO počítaný z výtokových křivek bude při vyšších celkových průtocích mírně nadhodnocen. Popsaným postupem získáme hodnoty ZO pouze v oblastech výtokových křivek. Chybějící průběh ZO ve vzestupných větvích a v oblastech kolísání průtoků doplníme prozatím konstantní hodnotou posledního stanoveného ZO (v době rychlého zvyšování průtoku roste výška hladiny v toku rychleji než výška hladiny vody v podzemních vodách a ZO tak může na krátký čas klesnout i do záporných hodnot – zásak vody z toku do podzemních vod) a výslednou lomenou čáru vyhladíme geometrickým průměrem. 8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
Na obr. 21 je uveden výsledek výpočtu ZO z výtokových křivek (fialově, vyhlazení oknem 19 dnů) a z vyhlazených klouzavých minim s oknem 19 dnů (žlutě).
y
e
10 9 8 7
3
6
4
2
1
denní průměry průtoků
0
/
-
.
,
)
*
5
vyrovnaný ZO z výtokových křivek
ZO zvyrovnaných klouzavých minim, okno 19
+
(
'
4
místo vyhovující podmíncevýtokové křivky
3 2 1 0 5000
5100
5200
5300
5400
5500
5600 y
5700
5800
5900
6000
Obr. 21: Vyhodnocení základního odtoku metodou klouzavých minim a metodou výtokových křivek Pokud by byly výtokové křivky tříděny i za pomoci 3. derivace, pravděpodobně by byl ZO z výtokových křivek nižší, především v průtokových vlnách. Mimo nejvyšší hodnoty průtokových vln by žádný velký rozdíl nebyl. Takto je ZO v oblasti průtokových vln pravděpodobně nadhodnocen. Tuto vadu by měl odstranit způsob zpracování pomocí funkčních měrných křivek. Pokusil jsem se také spočítat ZO z hodinových průtoků profilu Popov na Vláře pro rok 2003, jehož stopa pro nový způsob zpracování je na obr. 18. Výpočet byl prováděn iteračně a byl značně pracný a dá se říci, že iterace leckdy i selhávaly a nedávaly dobré výsledky. Jsou to ty úseky, které jsou konvexní (duté) místo konkávní (vypuklé). Grafická znázornění vidíme na obr. 22:
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
4 3.5 3 2.5
celkový průtok vypočtený ZO
2
1.5 1 0.5 0 0
2000
4000
6000
y
8000
Obr. 22: Výtokové části celkového odtoku a ZO v profilu Popov na toku Vlára, hodinový krok, rok 2003, nové zpracování Na obr. 22 je uveden nevyhlazený výsledek ZO jen v oblasti výtokových čar profilu, kde lze ZO přímo spočítat. Na počátku mezi 0 a 1000 hodin je vidět, jak krásně navazuje ZO na výtok vlny z deště pod následující průtokovou vlnou, která zvedla jen povrchový odtok, ale zřejmě ne podzemní. ZO je počítán postupně z průběhu modře kreslených celkových průtoků přímo nad ním. Pro posouzení kvality výpočtu ZO novou metodou uvádím na obr. 23 stopu ZO získanou ze ZO znázorněného fialovou čarou na obr.22.
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
y
4
2
0
#
!
"
-20
-15
-10
-5
data 0
omezení podzemního odtoku
-2
-4
-6
Obr. 23: Stopa ZO profilu Vlára, Popov, hodinové průtoky roku 2003, nové zpracování. U stopy hodinového ZO z obr. 23 je z hlediska hydrologa velmi užitečné porovnání se stopou celkového hodinového odtoku téhož profilu za týž rok 2003 na obr. 18. Je pěkně pozorovatelné odřezání strmějších větví mělčích druhů odtoku. Na tomto výpočtu založeném na hodinových hodnotách bude ještě třeba pracovat, protože postup je zatím pro automatizaci a tudíž i provozní užití nevyhovující a navíc iterační postupy mají tendenci selhávat. Je třeba také stanovit zdůvodněný postup výpočtu ZO pro jiné části hydrogramu, než jen pro výtokové větve.
4. Závěr Příspěvek popisuje některé postupy založené na statistickém i analytickém zpracování hydrologických dat, jejichž použití společně s nově zaváděnou technikou by mohlo vést ke zkvalitnění výstupů hydrologických zpracování jak povrchových tak podzemních vod. Z navrhovaných řešení některých problémů hydrologie je zatím prakticky použitelná vícerozměrná spojitě se měnící měrná křivka průtoků, která již byla využita při zpracování průtoků v Podhradí na Dyji, kde je stav vody ovlivňován při vysokých stavech vody v přehradě Vranov.Zde se model ukázal plně funkční a lze jej použít i na jiné případy proměnné MKP. Interpolace hydrogramu pomocí funkce, která je řešením Bernoulliovy rovnice, zatím nelze provozně použít (nedořešené některé problémy a nedostatek vhodných údajů), zatímco základní odtoky počítané alespoň z denních průměrů rychlou metodou klouzavými minimy by použitelné byly. Mnohem kvalitnější výpočet z hodinových průtoků je ovšem závislý na řešení problému interpolace, použití funkčních MKP a zpřehlednění převodu stavů hladin na průtoky. Literatura: [1] Kouřil Z.: Stanovení přírodních zdrojů podzemních vod podle metody G. Castanyho. In: Sborník přednášek. Výpočty využitelného množství podzemních vod. Brno, ČVTS 1975.
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006
[2] Budík L., Budíková M.: Statistické zpracování měsíčních a ročních srážkových a odtokových charakteristik povodí řeky Moravy. Práce a studie, sešit 28, ČHMÚ Praha, 2001.
8.
NÁRODNÍ KONFERENCE
STATISTICKÉ DNY
V
BRNĚ, BRNO 27. - 28.
ČERVNA,
2006