Karolina Pazio-Urbanowicz
Ocena możliwości wykorzystania wybranych metod geostatystycznych do generowania powierzchni modelowych
w analizach przepływu wód podziemnych
Praca magisterska wykonana pod kierunkiem: dr Katarzyny Ostapowicz
Praca magisterska złożona zgodnie z wymogami programu studiów magisterskich UNIGIS (Master of Science)
w zakresie "Geographical Information Science & Systems".
Uniwersytet Jagielloński w Krakowie, Uniwersytet Paris Lodron w Salzburgu
2015
2
Oświadczam, że wszystkie źródła wykorzystane w pracy zostały wymienione zgodnie z regułami cytowania. Niniejszą pracę wykonałam samodzielnie. Praca ta nie była i nie będzie składana jako praca dyplomowa w innym miejscu.
Data………... Podpis………
3
Spis treści
1. Wstęp 5
2. Cel pracy 7
3. Stan badań 8
3.1.Geostatystyka – krótka historia 8
3.2.Geostatystyka – podstawowe pojecia 12
3.3.Przegląd wybranej literatury z zakresu modelowania przepływu wód
podziemnych 16
3.4.Zasady schematyzacji hydrogeologicznej 17
4. Obszar testowy 21
4.1.Geologia obszaru 23
4.2.Model konceptualny obszaru badań 24
5. Dane 26
5.1.Dane z Centralnego Banku Danych Hydrogeologicznych HYDRO 26 5.2.Dane z Centralnego Ośrodka Dokumentacji Geodezyjnej
i Kartograficznej 28
6. Metody 29
6.1.Wstępne przetworzenie danych 29
6.2.Wybór próby punktów weryfikacyjnych 32
6.3.Ekploracyjna analiza danych 32
6.4.Interpolacja 34
6.4.1.Metoda odwrotnych odległości 35
6.4.2.Triangulacja 36
6.4.3.Kriging 37
6.5.Wariogram jako miara ciągłości przestrzennej 40
6.6.Ocena poprawności interpolacji 44
4
7. Wyniki 47
7.1.Wyniki wstępnego przetworzenia danych 47
7.2.Wyniki wyboru próby punktów weryfikacyjnych 49
7.3.Wyniki analizy ekspoloracyjnej danych 50
7.4.Interpolacja 54
7.4.1.Metoda odwrotnych odległości 54
7.4.2.Triangulacja 59
7.4.3.Kriging 61
7.5.Wyniki oceny poprawności interpolacji 75
7.5.1.Metoda odwrotnych odległości 75
7.5.2.Triangulacja 77
7.5.3.Kriging 78
8. Dyskusja 84
9. Wnioski 89
10. Spis rycin 91
11. Spis tabel 93
12. Literatura 94
5
1. Wstęp
Badania dotyczące przepływu wód podziemnych stanowią obecnie podstawowe narzędzie pracy hydrogeologa. Zakres i tematyka, w której mogą być wykorzystane są bardzo szerokie – od wyznaczania wielkości zasobów wód podziemnych w ujęciu lokalnym dla pojedynczych studni czy regionalnym dla całych zlewni, przez symulacje transportu zanieczyszczeń w warstwie wodonośnej, po projektowanie odwodnień na potrzeby budownictwa lądowego (Dąbrowski i in. 2010; Kresic, Mikszewski 2011).
Najczęściej badania modelowe wykorzystywane są w opracowaniach regionalnych przy określeniu odnawialności systemów wodonośnych, określeniu wielkości zasobów dyspozycyjnych (np. Oficjalska i in. 2010; Rodzoch i in. 2012) czy wyznaczaniu obszarów ochronnych Głównych Zbiorników Wód Podziemnych (np. Rodzoch i in. 2012). Modelowanie matematyczne wykorzystuje się także podczas przeprowadzania oceny wpływu istniejących lub projektowanych inwestycji na środowisko gruntowo-wodne, szczególnie w odniesieniu do obszarów ekosystemów zależnych od wód gruntowych (Dąbrowski i in. 2010). Badania modelowe wykorzystywane są również do prognozowania wpływów odwodnień górniczych i budowlanych na system wodonośny. Poprawne wykorzystanie metod modelowania numerycznego wymaga odpowiedniego rozpoznania systemu wodonośnego oraz stworzenia pojęciowego modelu tego systemu wraz z określeniem warunków brzegowych, rozpoznaniem parametrów hydrogeologicznych obiektu jak i rozkładu ciśnień piezometrycznych (Dowgiałło 2002).
System wodonośny można zdefiniować jako zespół poziomów wodonośnych znajdujących się w kontakcie hydraulicznym, ograniczony ściśle zdefiniowanymi przestrzennie i dynamicznie granicami (Dowgiałło 2002). Model natomiast możemy zdefiniować jako układ fizyczny lub opis matematyczny o właściwościach analogicznych do właściwości obiektu modelowanego, działający analogicznie do tego obiektu (Dowgiałło 2002). Numeryczny model przepływu wód podziemnych jest to algorytm numeryczny pozwalający na uzyskanie przybliżonego rozwiązania równania różniczkowego cząstkowego określonego procesu, w tym przypadku ogólnego równania przepływu (Dowgiałło 2002). Najczęstszymi metodami służącymi do uzyskania rozwiązania tego równania są metody różnic skończonych oraz elementów skończonych. W pierwszej z nich, algorytm obejmuje rozwiązanie układu równań liniowych rozpisanych dla wszystkich węzłów zdyskretyzowanego obszaru, będących
6
jednocześnie analogami różnicowymi wyjściowego równania różniczkowego.
W drugiej algorytm opiera się na poszukiwaniu minimalnej wartości funkcji wysokości hydraulicznej. W obu przypadkach warunkiem uzyskania rozwiązania jest znajomość:
formy przestrzennej obszaru, jego struktury, rozkładu parametrów, warunków brzegowych i początkowych (Dowgiałło 2002).
Niezależnie od celu, jaki chcemy osiągnąć konstruując numeryczny model przepływu wód podziemnych, najważniejszym etapem pracy jest właściwa schematyzacja warunków hydrogeologicznych. Obejmuje ona wyznaczenie granic obszaru modelowego, wydzielenie poszczególnych warstw wodonośnych oraz określenie ich parametrów hydrogeologicznych (Staśko i in. 2010). Zrozumienie wzajemnych zależności pomiędzy czynnikami kształtującymi geometrię środowiska gruntowo-wodnego a elementami wyrażającymi stan tego systemu wymaga nie tylko gruntownej znajomości opracowań archiwalnych dostępnych dla analizowanego obszaru, ale i przeprowadzenia szeregu analiz przestrzennych dotyczących m.in.
rozkładu parametrów hydrogeologicznych, takich jak: współczynnik filtracji, przewodność, głębokość występowania warstwy wodonośnej czy miąższość nadkładu słaboprzepuszczalnego. Niezwykle przydatnym narzędziem pozwalającym na analizę danych oraz tworzenie map wejściowych rozkładu poszczególnych parametrów, którymi dysponujemy w procesie schematyzacji systemu hydrogeologicznego, jest analiza przestrzenna. Zastosowanie metod geostatystycznych pozwala na interpretację oraz przewidywanie (estymację) rozkładu przestrzennego badanego zjawiska. Jedną z wielu metod z zakresu analizy przestrzennej jest interpolacja przestrzenna. Procedura ta pozwala na szacowanie wartości cechy w nieopróbowanych punktach na obszarze objętym istniejącymi pomiarami i wykorzystywana jest powszechnie do uzupełnienia siatki pomiarowej lub też konwersji pomiarów punktowych na poligony wypełniające obszar badawczy. Obecnie dysponujemy szeregiem metod wyznaczających funkcje interpolacyjne. W niniejszej pracy szerzej zostaną omówione metody najczęściej wykorzystywane w praktyce hydrogeologicznej – kriging, metoda odwrotnych odległości i triangulacji.
7
2. Cel pracy
Model numeryczny w konkretnych przypadkach stanowi narzędzie opisu systemu wodonośnego i jego zachowania w warunkach wymuszeń. Poprawność działania modelu zależy przede wszystkim od prawidłowego wyznaczenia jego granic.
W praktyce geologicznej wyznaczanie granic, zarówno „pionowych” wydzielających modelowany system wodonośny z pewnego większego continuum, jak i granic
„poziomych” stanowiących powierzchnie poszczególnych warstw modelu, odbywa się często w sposób analogowy, poprzez ręczne kreślenie izolinii, a następnie ich digitalizację. W większości przypadków ilość informacji geologicznej jest wystarczająca do przeprowadzenia procesu kreślenia map powierzchni w sposób zautomatyzowany.
Celem niniejszej pracy jest ocena możliwości wykorzystania wybranych metod geostatystycznych do przygotowania mapy powierzchni geologicznej stanowiącej jedną z powierzchni modelowych numerycznego modelu przepływu odwzorowującego wybrany system wodonośny wraz z oceną dokładności otrzymanych wyników. W pracy zostaną przedstawione wyniki interpolacji za pomocą metody odwrotnych odległości, krigingu oraz triangulacji. Ocena stopnia dokładności uzyskanych w wyniku interpolacji map powierzchni modelowej zostanie przeprowadzona na podstawie analizy wartości błędów oraz poprzez porównanie uzyskanych powierzchni. Jednocześnie, w pracy zostanie przedstawiona próba zebrania ogólnych cech, jakie powinny charakteryzować powierzchnie modelowe.
8
3. Stan badań
W niniejszym rozdziale przedstawione zostały zagadnienia z zakresu geostatyki oraz matematycznego modelowania przepływu wód podziemnych. Dokonany został przegląd literatury z zakresu modelowania przepływu wód podziemnych. Wyjaśnione zostaną zasady przeprowadzania schematyzacji warstw wodonośnych na potrzeby modeli numerycznych.
3.1. Geostatystyka – krótka historia
Statystyka przestrzenna to zbiór metod analiz danych, które oprócz cech obiektów uwzględniają ich lokalizację w przestrzeni i/lub czasie (Stach 2008).
O pierwszym przypadku wykorzystania analizy przestrzennej możemy mówić w odniesieniu do badań angielskiego lekarza, prekursora epidemiologii, Johna Snow.
W 1854 roku, analizując przyczyny epidemii cholery w Londynie, sporządził mapę przedstawiającą miejsca zgonów oraz miejsca poboru wody. Analizując rozmieszczenie tych miejsc, doszedł do wniosku, że woda z jednego z ujęć stanowi prawdopodobne źródło epidemii (Longley i in. 2008). Można przyjąć, że Snow był nieświadomym prekursorem zarówno wykorzystania systemów informacji geograficznej (GIS), jak i budowy modeli konceptualnych.
Geostatystka jest zbiorem narzędzi statystycznych uwzględniających w analizie danych ich przestrzenną lokalizację. Analizy te oparte są o teorię funkcji losowych (Wackernagel i Bertino 2005). Wyniki pomiarów traktowane są w geostatystyce jako pewne realizacje funkcji losowych. Wyodrębnienie geostatystyki jako odrębnej gałęzi statystyki przestrzennej zajmującej się badaniem zjawisk losowych z uwzględnieniem wymiaru przestrzennego nastąpiło w latach sześćdziesiątych XX wieku (Zawadzki 2011). Podwaliny geostatystyki powstawały jednocześnie w kilku miejscach.
W Afryce Południowej za sprawą prac złożowych Daniela Krige’a (1951), który przedstawił ideę przewidywania zawartości złota w złożu na podstawie punktowego opróbowania złoża. Rosyjski matematyk Kołmogorow (1941) przedstawił funkcję, która jednocześnie przedstawia korelację przestrzenną (obecny wariogram) oraz metodę interpolacji (kriging) zmiennej. W Szwecji Bertil Matérn pokazał dopuszczalne metody analiz opisujących kowariancję. W Wielkiej Brytanii Geoffrey Jowett (1955) przedstawił opis podstaw obecnego wariogramu. Podstawy teoretyczne geostatystyki
9
zgromadził i zebrał Georges Matheron, założyciel Centrum Geostatystyki w Fointainebleau (Francja), w latach 1962-1963 publikując prace zawierające teoretyczne podstawy metod opisu zmiennych wykorzystywanych przy szacowaniu wielkości zasobów złóż (Agterberg 2014).
Zastosowanie metod geostatystycznych przyczyniło się do zmiany sposobu szacowania zasobów i jakości złóż. Wprowadzenie instrumentu matematycznego pozwalającego na uwzględnienie struktury zmienności parametrów w przestrzeni (semiwariogramu) pozwoliło na dokładne oszacowanie wielkości złoża w miejscach nieobjętych opróbowaniem i odejście od założenia o całkowitej losowości rozkładu parametrów złożowych. Procedurę szacowania parametrów złoża – jego wielkości i zasobności, nazwano krigingiem na cześć ww. Daniela Krige’a. Rozpowszechnienie krigingu wywołało szereg uwag krytycznych co do metody, a to z kolei zaowocowało powstaniem kriging uniwersalnego (Huijbregts, Matheron 1971 za Mucha, Wasilewska- Błaszczyk 2010) uwzględniającego występowanie trendu w rozkładzie przestrzennym parametrów. Dalszy rozwój geostatystyki związany jest z opracowaniem metod nieliniowych (tab. 3.1) oraz innych metod interpolacji, np. metody odwrotnych odległości (Harvard, lata 60. XX wieku). Podstawy triangulacji Delunay zostały opracowane przez rosyjskiego matematyka Borysa Delone w 1934 roku.
Tab. 3.1. Twórcy wybranych rodzajów krigingu (Źródło: Mucha, Wasielewska-Błaszczyk 2010)
Rodzaj krigingu Twórca metody Rok opisania metody
Kriging rozłączny Matheron 1976
Kriging indykatorowy Journel 1983
Kriging probabilistyczny Sullivan 1984
Ze względu na pierwsze kierunki wykorzystania geostatystyka do końca lat 80.
XX wieku była używana głównie w pracach geologicznych, ze szczególnym naciskiem na geologię złożową. Z czasem jednak, wraz ze zwiększeniem dostępności komputerów osobistych oraz rozwojem komercyjnego oraz darmowego oprogramowania, zaczęła być coraz częściej wykorzystywana w innych naukach przyrodniczych: w kartografii gleb m.in. do generowania map glebowych w oparciu o wyniki badań terenowych (Malczyk, Długosz i Radke 2006), do charakteryzowania pokrywy glebowej (Długosz, Jaworska, Ockerman 1999); w meteorologii m.in. do określania zmian temperatury globalnej (Cowtan, Way 2014); w oceanologii do opracowywania map siedlisk,
10
w ekologii do określania skażenia środowiska (Pająk, Szostak, Sławiński 2012).
Narzędzia geostatystyczne wykorzystywane są między innymi przy:
- identyfikacji struktury zmienności przestrzennej i/lub czasowej cechy bądź zjawiska;
- szacowaniu wartości cechy ilościowej w nieopróbowanym punkcie w przestrzeni i/lub momencie czasu;
- szacowaniu wartości średniej cechy ilościowej dla określonej powierzchni i/lub okresu czasu;
- ocenie błędu szacunku wartości punktowej i/lub obszarowej cechy ilościowej;
- szacowaniu wartości ekstremalnych cechy możliwych w danym punkcie, czy obszarze;
- obliczeniach prawdopodobieństwa przynależności danego punktu lub obszaru do określonej kategorii;
- obliczeniach prawdopodobieństwo przekroczenia w danym punkcie lub obszarze wartości progowej cechy ilościowej;
- filtrowaniu składowych przestrzennych i/lub czasowych cechy lub zjawiska;
- separacji tzw. „szumu” (składowej losowej);
- optymalizacji próbkowania i projektowania sieci monitoringowych.
Zainteresowanie geostatystyką wzrastało lawinowo z roku na rok, co odnotowano w liczbie artykułów naukowych z nią powiązanych, i obecnie utrzymuje się na stałym (i dosyć wysokim) poziomie (Hengl, Minasny, Gould 2009). Istnieje wiele przeglądowych publikacji na temat geostatystyki m.in. Deutsch (2002), Gòmez- Hernández (2006). Publikacja autorstwa Issacs’a i Srivatsatavy (1989) An Introduction to Applied Geostatistics zawiera przystępny i przejrzyście zaprezentowane opis podstaw geostatystyki. natomiast podręcznik Webster’a i Oliver’a (2008) Geostatistics for Environmental Scientists omawia również wiele aspektów z zakresu nauk o Ziemi.
Bardzo uboga natomiast jest literatura z tego zakresu w języku polskim. Do tej pory w zasadzie można mówić o dwóch istotnych pozycjach literatury: Geostatystyka. Teoria i praktyka Barbary Namysłowskiej-Wilczyńskiej (2006) oraz Metody geostatystyczne dla kierunków przyrodniczych i technicznych Jarosława Zawadzkiego (2011), obu prezentujących podstawowe zagadnienia z zakresu geostatystyki. Stosunkowo wiele publikacji w formie artykułów dotyczących geostatystyki popełnił Stach. Dotyczyły one analizy struktury pola opadów atmosferycznych na terenie Polski (2005, 2007b, 2009),
11
teorii analizy przestrzennej danych jakościowych (2006), analizy autokorelacji przestrzennej i czasowej danych paloeekologicznych (2007a) czy składu petrograficznego osadów fluwioglacjalnych (2008 z Górska-Zabielska).
Oczywiście przedstawione publikacje nie wyczerpują literatury tematu. Szereg informacji na temat praktycznego wykorzystania metod geostatystycznych można również znaleźć w instrukcjach użytkowania programów komputerowych jak choćby Surfer firmy Golden Software (Surfer. User’s Guide 2002) czy modułu Geostatistical Analyst oprogramowania ArcGIS (Johnston i in. 2001). Szereg podstawowych informacji można znaleźć na stronach internetowych, m.in. Centrum Geostatystyki na Paryskiej Szkole Górniczej (http://www.geosciences.mines-paristech.fr) czy Centrum Szacowania Złóż na Uniwersytecie Stanforda (https://pangea.stanford.edu/re searchgroups/scrf/). Na stronie www.spatialanalysisonline.com dostępna jest również wersja on-line podręcznika autorstwa zespołu: Michael de Smith, Michael Goodchild i Paul Longley stanowiącego wszechstronny przewodnik po zagadnieniach analiz geosprzestrzennych.
Wraz ze wzrostem zainteresowania wykorzystaniem metod geostatystycznych w różnych dziedzinach nastąpił wzrost dostępności tych metod poprzez różnorodne rodzaje oprogramowania komputerowego. Podstawowe funkcje geostatystyczne dostępne są jako wbudowane procedury lub moduły dodatkowe w GIS-owych programach komercyjnych takich jak ArcGIS firmy ESRI (www.esri.com), MapInfo firmy PitneyBowes Software (www.mapinfo.com) czy Geomedia firmy Intergraph (www.intergraph.com), a także programach typu open source: SAGA (www.saga- gis.org), R (www. r-project.org) czy Quantum GIS (www.qgis.org). Szczególnie warte uwagi jest specjalistyczne oprogramowanie dla geologów, zawierające szereg wbudowanych funkcji geostatystycznych - GMS firmy Aquaveo (www.aquaveo.com), GoCAD firmy Paradigm (www.pdgm.com) czy wykorzystywany nie tylko w geologii program Surfer firmy Golden Software (www.goldensoftware.com).
Opracowania geostatystyczne dotyczące zagadnień geologicznych związane są przede wszystkim z geologią złożową i problematyką szacowania wielkości złóż oraz ich geometrii. Prace przeglądowe z tego zakresu publikowane były przez Davisa (1986), czy Journel i Huijbregts (2004), a w Polsce przez Muchę (1994). Pierwsze polskojęzyczne prace z zakresu geostatystyki górniczej dotyczyły opisu zmienności zawartości rud cynku i ołowiu w złożu rud Zn-Pb Bolesław za pomocą
12
semiwariogramów (Mucha 1978, 2010). W kolejnych latach zainteresowanie geostatystyką było umiarkowane z powodu braku odpowiedniego sprzętu obliczeniowego oraz oprogramowania, a w Polsce dodatkowo dzięki ograniczonemu dostępowi do literatury światowej przedmiotu. Dopiero upowszechnienie komputerów osobistych na przełomie lat 80. i 90. XX wieku doprowadziło do zwiększonego zainteresowania geostatystyką, przede wszystkim w geologii górniczej. Niemniej jednak liczba prac z tego zakresu wciąż jest niewielka i skupia się głównie na szacowaniu wartości zasobów i parametrów w złożu różnymi metodami (Kokesz 2006; Kokesz, Nieć 1992; Mucha 2002; Namysłowska-Wilczyńska 1993).
Z niewielu prac choćby częściowo związanych z hydrogeologią można wymienić opracowanie Helsena i Hirsha (1992) Statistical methods in water resources czy Deutscha (2002) Geostatistical reservoir methods. Bardzo nieliczne są również opracowania stricte hydrogeologiczne zawierające elementy analizy geostatystycznej.
Są to przeważnie artykuły problemowe z zakresu analizy jakości wód autorstwa Kani (2011), Zawadzkiego, Kucharka i Traichela (2004) czy artykułu dotyczące numerycznego modelowania przepływu wód podziemnych, jednak w tych ostatnich analiza geostatystyczna jest zwykle jedynie wymieniona jako jedno z zastosowanych narzędzi i nie jest szczegółowo opisywana i analizowana. W zasadzie można powiedzieć, że brak jest szczegółowej literatury, szczególnie polskojęzycznej, dotyczącej sposobu wykorzystania metod geostatystycznych w hydrogeologii, a w szczególności przy budowie modeli przepływu.
3.2. Geostatystyka – podstawowe pojecia
Podstawowym pojęciem geostatystyki, zdefiniowanym przez Matherona (1963) jako funkcja, która przyjmuje określoną wartość w każdym punkcie przestrzeni trójwymiarowej, jest pojęcie zmiennej zregionalizowanej. Zmienne te charakteryzują się wzajemną korelacją związaną z przestrzennym rozkładem danego zjawiska, a jednocześnie mogą być traktowane lokalnie jako zmienne losowe o określonym rozkładzie prawdopodobieństwa. Zmienne zregionalizowane wyraża funkcja losowa zmiennej przestrzennej Z(x) charakteryzująca się tym, że:
- lokalnie wartość z(x) funkcji losowej Z(x) jest) traktowana jest jako realizacja funkcji zmiennej losowej;
13
- w dowolnej parze punktów xi i xi+h zmienne losowe Z(xi) oraz Z(xi+h) związane ze sobą zależnością korelacyjna wynikającą z przestrzennej ciągłości badanego zjawiska (Zawadzki 2011).
Za pomocą zmiennej zregionalizowanej można przedstawić np.: zanieczyszczenie pewną substancją na określonym obszarze, zawartość pierwiastka w rudzie w danym złożu, temperaturę morza na danym obszarze, głębokość zalegania wybranych utworów geologicznych i wiele tym podobnych.
Istnieje przynajmniej kilkanaście różnorodnych metod interpolacyjnych, które można klasyfikować jako:
- metody lokalne i metody globalne;
- metody wierne i metody wygładzające;
- metody zakładające ciągłość powierzchni oraz te dopuszczające ich nieciągłość;
- metody deterministyczne i metody stochastyczne.
Wybór metody globalnej sprowadza się do zastosowania pojedynczej funkcji matematycznej do danych ze wszystkich punktów pomiarowych, co w praktyce daje powierzchnie o „wygładzonym” przebiegu, pozbawione lokalnych szczegółów.
Wykorzystanie metod lokalnych uwzględnia wielokrotne wykorzystanie pojedynczej funkcji w wydzielonym lokalnym podzbiorze, a następnie połączenie wyników lokalnej interpolacji na całym obszarze, co w efekcie daje szczegółowy obraz zmienności przestrzennej zjawiska. Przykładem metod estymacji lokalnej (punktowej) jest kriging zwyczajny, metoda odwrotnych odległości, metoda najbliższego sąsiedztwa.
Przykładem metod globalnych jest kriging blokowy (Urbański 2010).
Metody wierne uwzględniają wszystkie dane pomiarowe, tak, że znajdą się one dokładnie na wyinterpolowanej powierzchni. Metody ta powinna być stosowana przy absolutnej pewności poprawności danych, którymi dysponujemy. Metody wygładzające z kolei nie uwzględniają ściśle danych pomiarowych i są wskazane w sytuacji niepewności co do jakości danych pomiarowych. Przykładem metody wiernej jest metoda odwrotnych odległości, metoda najbliższego sąsiedztwa oraz kriging zwykły (Urbański 2010).
Metody ciągłe dają wygładzone powierzchnie między punktami i wskazane są do interpolacji danych charakteryzujących się małą zmiennością lokalną natomiast
14
metody nieciągłe dają powierzchnie o charakterze tarasowym i wskazane są do interpolacji danych charakteryzujących się dużą zmiennością lokalną. Przykładem metody nieciągłej jest metoda najbliższego sąsiedztwa. Przykładem metody ciągłej jest kriging zwykły (Urbański 2010).
Metody deterministyczne pozwalają na tworzenie modelu jako jednoznacznie określonej powierzchni matematycznej, natomiast metody stochastyczne umożliwiają uwzględnienie w interpolowanej powierzchni zmienności losowej. Przykładem metody deterministycznej jest metoda odwrotnych odległości, metoda najbliższego sąsiedztwa.
Przykładem metody stochastycznej jest kriging (Urbański 2010).
Tworząc reprezentację środowiska za pomocą metod interpolacji nie można zapominać, że:
- tworzona reprezentacja, dotyczy pojedynczych punktów przestrzeni, w których dokonano pomiarów/wykonano badania;
- reprezentacja jest wybiórcza, a zatem niepełna;
- w tworzeniu reprezentacji wykorzystujemy model środowiska przyrodniczego w postaci pól ciągłych lub przestrzeni wypełnionej przez fizycznie istniejące i dobrze zdefiniowane obiekty;
- właściwości sąsiedztwa są podstawą zrozumienia zmienności przestrzennej oraz metod uzupełniania brakujących danych przestrzennych;
- na poprawność stworzonej reprezentacji wpływ ma zagadnienie skali i stopnia szczegółowości;
- każda stworzona reprezentacja ze względu na fakt, że dotyczy punktów i jest wybiórcza, jest również niepewna.
Wybór metody interpolacji powinien być poprzedzony szczegółową analizą zbioru danych którymi dysponujemy. Pierwszym krokiem jest wizualny przegląd dostępnych danych, a następnie wizualna próba znalezienia wartości odstających.
Przeprowadzenie eksploracyjnej analizy danych przestrzennych (exploratory spatial data analysis ESDA) określającej ilościowe cechy danych pozwoli na zrozumienie trendów jakim podlega opracowywany zbiór. Podstawowymi zadaniami ESDA jest określenie: rozkładu danych, globalnych i lokalnych wartości odstających, zróżnicowania danych oraz autokorelacji. Można to wykonać sporządzając histogram rozkładu, wykresy kwanty-kwantyl (qq-plot), analizując trendy zbioru danych,
15
sporządzając diagram Voronoi (Voronoi map) czy wariogram. Nie wszystkie z wymienionych kroków są konieczne do określenia cech dostępnego zestawu danych.
Ostatecznie na wybór metody interpolacji znaczący wpływ powinien mieć rozkład danych. W przypadku nierównomiernie rozmieszczonych danych istotny jest również stopień homogeniczności rozprzestrzenienia danych, liczba punktów na jednostkę powierzchni, wariancja populacji oraz cel przeprowadzania interpolacji (Goldsztejn, Skrzypek 2004). Ocena wiarygodności uzyskanej interpolacji powinna opierać się na analizie wielkości błędów interpolacyjnych oraz ocenie czytelności, realistyczności i wiarygodności uzyskanej mapy. Szczegółowy opis rodzajów błędów interpolacyjnych zostanie przedstawiony w rozdziale 6.6.
Najpopularniejszymi metodami interpolacji są: triangulacja lub metoda trójkątów, metoda minimalnej krzywizny, metoda najbliższego sąsiada, metoda naturalnych sąsiadów, metoda odwrotnych odległości, interpolacja wielomianowa, metoda Sheparda oraz kriging.
Triangulacja (Thiessen (Voronoi) polygons method) generuje powierzchnię zbudowaną z trójkątów, co w przypadku niewystarczającej liczby punków lub ich nierównomiernego rozkładu tworzy powierzchnię niewygładzoną. Metoda ta nie generuje wartości spoza zakresu danych. Wadą tej metody jest skokowa zmiana wartości zmiennej na granicy wieloboków. Metoda minimalnej krzywizny (minimum curvature) generuje powierzchnie o określonym stopniu wygładzenia (metoda wygładzająca). W metodzie tej możliwe jest uzyskanie wartości spoza zakresu danych, szczególnie w przypadku obszarów pozbawionych pomiarów możliwe jest powstanie znaczących artefaktów. Metoda najbliższego sąsiada (nearest neighbour method) jest szczególnie przydatna w przypadku regularnie lub prawie regularnie rozmieszczonych danych. Nie generuje ona wartości spoza zakresu. Metoda naturalnych sąsiadów (natural neighbour method) daje dobre wyniki przy interpolacji nierównomiernie rozmieszczonych danych. Nie daje wyników w obszarach pozbawionych danych i spoza zakresu. Metoda odwrotnych odległości (inverse distance weighting method) nie generuje wartości spoza zakresu. Jest także interpolatorem dokładnym (metoda wierna).
Jej wadą jest tworzenie koncentrycznych stref wokół punktów (tzw. bull’s eye).
Interpolacja wielomianowa (polynomial regression method) jest wykorzystywana do analizy trendu, dobrze radzi sobie z dużą ilością danych, ale daje powierzchnie o wysokim stopniu generalizacji. Może generować wartości spoza zakresu. Metoda
16
Sheparda (Shepard’s method) jest metodą zbliżoną do metody odwrotnych odległości jednakże nie generuje artefaktów. Możliwe jest za to uzyskanie wartości spoza zakresu.
Kriging jest jedną z bardziej elastycznych metod interpolacji. Przeważnie daje wiarygodne i satysfakcjonujące wyniki interpolacji. Możliwe jest uzyskanie wartości interpolowanych spoza zakresu danych (Galon 2014).
3.3. Przegląd wybranej literatury z zakresu modelowania przepływu wód podziemnych
Badania modelowe stały się jedną z metod rozwiązywania problemów hydrogeologicznych w połowie lat 60. XX wieku. Początkowo stosowano metody analogowe, potem metody modelowania matematycznego (Dąbrowski i in. 2010)
Obecnie modelowanie numeryczne przepływu wód podziemnych stanowi potężne narzędzie hydrogeologiczne i jest wykorzystywane w większości prac związanych z tą dziedziną. Literatura polska, jak i zagraniczna, obfituje w artykuły związane z szeroko pojętym modelowaniem przepływu wód podziemnych, jednak dominuję przedstawienie wyników modelowania i wniosków, głownie dotyczących ilości zasobów wód podziemnych czy ich rozmieszczenia (np. Kania, Oszczypko, Witczak 2009; Marciniak i in. 2013; Niedbalska 2013). Nieliczne są prace z zakresu sposobu przeprowadzania schematyzacji warstw, budowy modelu konceptualnego czy techniczne aspekty budowy samego modelu (np. Gurwin, Serafin 2008, 2010; Kania , Oszczypko, Witczak 2009). Opis podstaw procesu numerycznego modelowania przepływu wód podziemnych, sposobu działania poszczególnych warunków brzegowych można znaleźć w publikacjach Pazdry, Kozerskiego (1990), Szymanki (1977, 1980), Michalaka i in. (2010), Kulmy i Zdechlika (2009). Najnowsze publikacje z tego zakresu to publikacje przeglądowe w których podjęto próbę przybliżenia metod modelowania matematycznego w badaniach hydrogeologicznych (Dąbrowski i in. 2010) oraz próbę przedstawienia sposobów schematyzacji warunków hydrogeologicznych na potrzeby numerycznego modelowania przepływu wód podziemnych w odniesieniu do jednolitych części wód podziemnych (Michalak, Nawalany i Sadurski 2011). Obie publikacje stanową polskojęzyczną odpowiedź na kompleksowy podręcznik modelowania przepływu Applied groundwater modelling (Anderson 1992). Pomimo tego, że wszystkie wymienione publikacje w mniejszym lub
17
większym stopniu skupiają się na metodyce modelowania przepływu wód podziemnych, w żadnej z nich nie zamieszczono konkretnych wskazówek dotyczących sposobu konstruowania warstw modelowych.
3.4. Zasady schematyzacji hydrogeologicznej
Modelowanie matematyczne jest powszechnie stosowanym narzędziem w realizacji zadań badawczych oraz przy wykonywaniu dokumentacji hydrogeologicznych. Najczęściej badania, w których wykorzystywane jest modelowanie matematyczne wykonywane są na potrzeby tzw. dokumentacji zasobowych i obejmują ustalenie zasobów ujęć wód podziemnych, określania obszarów zasilania ujęcia lub zlewni czy wielkości przewidywanego zasięgu obniżenia zwierciadła wód podziemnych. Prawidłowe sformułowanie zadania badawczego decyduje o zasadności budowy modelu, jego obszarze, a także stopniu szczegółowości.
Ma zatem bezpośrednie przełożenie na wymagany stopień rozpoznania budowy geologicznej modelowanego obszaru, a następnie przełożenie tej budowy na warunki modelowe. Poprawne wykonanie modelu matematycznego wymaga właściwego rozpoznania systemu wodonośnego, budowy modelu konceptualnego oraz określenia przestrzennej zmienności parametrów hydrogeologicznych i zdefiniowania warunków brzegowych. System wodonośny można określić jako pewnego rodzaju ideę całościowego przedstawienia formy występowania i dynamiki zbiorowiska wód podziemnych, jego składu chemicznego i własności fizyko-chemicznych oraz związków przyczynowo-skutkowych występujących pod wpływem czynników naturalnych i antropogenicznych w granicach obiektu przestrzennego tworzącego strukturę hydrogeologiczną (Szymanko 1980). Jest to jedna z wielu stosowanych definicji systemu wodonośnego. Jej zaletą jest możliwość przyjęcia, że system wodonośny stanowi pewien wydzielony i względnie odosobniony obiekt w którym następują przemiany indukowane pod wpływem zmian w otoczeniu oraz na to otocznie wpływające. Stworzenie modelu matematycznego systemu wodonośnego wymaga dokonania odpowiednich uproszczeń tego systemu, a więc budowy modelu konceptualnego będącego uproszczoną reprezentacją modelowanej rzeczywistości.
Model konceptualny jest opisem struktur systemu wodonośnego oraz przebiegu procesów w nim zachodzących z uwzględnieniem kontaktów z otoczeniem (Dąbrowski i in. 2010). Obecnie model pojęciowy (konceptualny) zastępuje się
18
terminem schemat warunków hydrogeologicznych. W literaturze polskiej wciąż brakuje spójnego sposobu określania czym dokładnie różnią się od siebie model koncepcyjny, model konceptualny, model pojęciowy i schemat warunków hydrogeologicznych.
Ryc. 3.1. Uproszczony schemat konstrukcji modelu numerycznego
W niniejszej pracy przyjęto więc że model konceptualny to inaczej model hydrogeologiczny, będący wynikiem schematyzacji rzeczywistych warunków hydrogeologicznych, natomiast jego dalsza schematyzacja prowadzi do powstania modelu matematycznego (Dąbrowski i in. 2010) Na model konceptualny systemu hydrogeologicznego składają się następujące parametry:
- budowa przestrzenna systemu hydrogeologicznego;
- zestaw parametrów hydrogeologicznych tj. wartości współczynnika filtracji w otworach, wydajności studni, znajomość litologii, rozmieszczenie wymuszeń tj. rzek, zbiorników wodnych oraz obiektów antropogenicznych: studni, zapór itd.
W niniejszej pracy nacisk zostanie położony na samą budowę struktury systemu wodonośnego i jej przekształcenie w układ warstw modelowych. Schematyzacja systemu wodonośnego w pierwszym etapie będzie więc opierała się na wyborze koncepcji tworzenia warstw modelowych. Po rozpoznaniu rozprzestrzenienia utworów
Model matematyczny
(wyznaczenie poszczególnych powierzchni modelowych, zadanie warunków brzegowych oraz warunku początkowego, wprowadzenie parametrów hydrogeologicznych)
Model konceptualny
(określenie granic modelowanego systemu, rozpoznanie parametrów hydrogeologicznych takich jak współczynnik filtracji, przewodność warstw
wodonośnych, określenie wymuszeń tj. drenażu czy zatłaczanie wód) Modelowana rzeczywistość (system wodonośny)
(analiza materiałów archiwalnych, w tym: map, dokumentacji regionalnych i lokalnych)
19
wodonośnych, słaboprzepuszczalnych i praktycznie nieprzepuszczalnych oraz granic tych utworów, jak i ich wieku podjęta powinna zostać decyzja o sposobie wydzielania warstw modelowych. Warstwy modelowe mogą być wydzielane stratygraficznie lub być zagregowane według litologii, a zatem i przepuszczalności skał. Źródłem danych do powstania modelu konceptualnego są dane otworowe zgromadzone w Centralnym Banku Danych Hydrogeologicznych (CBDH HYDRO), na Szczegółowej mapie geologicznej Polski w skali 1 : 50 000 i w archiwalnych dokumentacjach hydrogeologicznych. Tworzenie warstw wejściowych opiera się na generowaniu regularnych siatek wartości tzw. gridów, wykorzystywanych w obliczeniach metodą różnic skończonych lub generowaniu siatki trójkątów wykorzystywanej w obliczeniach prowadzonych metodą elementu skończonego. Wszystkie kolejne warstwy generowane są w sposób analogiczny zachowując jednolity układ współrzędnych danych wejściowych, krok siatki oraz punkt początkowy modelu.
Podstawą budowy modelu jest przestrzenne wyodrębnienie elementów składowych systemu dla którego jest on budowany. Schematyzacja warunków hydrogeologicznych jest najważniejszym etapem budowy modelu hydrogeologicznego (Michalak, Nawalany, Sadurski 2011). Niemniej jednak warto zauważyć, że jej efekt uzależniony jest w znacznej mierze od doświadczenia i poziomu warsztatu zawodowego autora, a więc jest w wysokim stopniu subiektywny.
Generalizacja ilościowa stosowana przy budowie schematu warunków hydrogeologicznych obejmuje uproszczenie przebiegu linii konturowych ograniczających wyróżnienia hydrogeologiczne i może dotyczyć zasięgu występowania poszczególnych utworów wyznaczających przebieg powierzchni stropowych i spągowych. Pominięcie niewielkich struktur, jak choćby soczewek piasków, jeśli ma to uzasadnienie merytoryczne, bo nie ma znaczenia dla analizowanego bilansu wód, jest w pełni uzasadnione. Generalizacja jakościowa opiera się na grupowaniu wyróżnień o podobnych parametrach (agregacji). Niemniej istotna jest generalizacja przestrzenna prowadzona na potrzeby modelu komórkowego, gdzie przestrzeń ciągła jest zastąpiona poprzez teselację siecią komórek i/lub dyskretyzację skończoną liczbą węzłów. Proces ten jest konieczny przy opracowywaniu modelu przepływu z zastosowaniem metody różnic skończonych.
Szczególnie istotne jest wydzielenie pierwszego od powierzchni terenu poziomu wodonośnego, który decyduje o funkcjonowaniu ekosystemów zależnych od wód
20
podziemnych (Dąbrowski i in. 2010). Obecnie zakłada się, że mniejszy błąd metodyczny popełnia się modelując pierwszy poziom wodonośny jako warstwę jednorodną niż pomijając ten poziom lub agregując go z poziomami położonymi głębiej.
Schematyzacja modelu powinna zostać poprzedzona ścisłym określeniem granic systemu hydrogeologicznego, a następnie określeniem warstwowości modelu, tj. zastąpieniem przestrzennie zorientowanego strumienia wód podziemnych odpowiadającym mu strumieniem płaskim w wydzielonych warstwach. Rozpoznanie lokalnych i regionalnych układów krążenia wód oraz ustalenie stopnia więzi hydraulicznej między poszczególnymi warstwami, więc rozpoznanie warstw wodonośnych-przewodzących i warstw izolujących, pomaga w przeprowadzeniu schematyzacji.
Nieprawidłowe wyznaczenie powierzchni stropowych/spągowych warstw może być przyczyną nieprawidłowości w obliczonym położeniu zwierciadła wód podziemnych. Gdy zwierciadło zostanie obniżone poniżej spągu warstwy przewodzącej, w kolejnym kroku iteracyjnym da to ujemną wartość przewodnictwa, co nie jest możliwe (Dąbrowski i in. 2010). Podobna sytuacja ma miejsce, gdy obliczone zwierciadło ma rzędną tylko niewiele ponad spąg warstwy wodonośnej. W takim przypadku przewodnictwo warstwy jest bardzo małe, a co za tym idzie przy kolejnych iteracjach rzędne zwierciadła będą zbliżały się do rzędnej spągu. Proces ten będzie się powtarzał i prawdopodobnie nie uzyskamy rozwiązania równania przepływu.
Struktura modelu matematycznego, czyli przestrzenny rozkład warstw wodonośnych biorących udział w krążeniu wód, określana jest podczas procesu schematyzacji. Model hydrogeologiczny (konceptualny) opisuje ciągłe pole filtracji, natomiast model matematyczny przekształca model o parametrach ciągłych na model o parametrach dyskretnych, określonych w konkretnych punktach. Analogicznie rozwiązanie również określone jest tylko w tych punktach. Struktura modelu matematycznego powinna odzwierciedlać wielowarstwowość systemu hydrogeologicznego. Nie można zastąpić rzeczywistego wielowarstwowego układu krążenia wód podziemnych jednowarstwowym modelem o uśrednionych parametrach, gdyż trudno wtedy o uzyskanie realnych wyników modelowania i określenie stref spływu czy stref zasilania.
21
4. Obszar testowy
Na obszar testowy wybrano zlewnie Białej Przemszy i Przemszy położone na terenie województw śląskiego i małopolskiego, zajmujące powierzchnię 2212,8 km2 (ryc. 4.1).
Ryc. 4.1. Obszar testowy
Zgodnie z podziałem fizycznogeograficznym według Kondrackiego (2011) obszar badań należy do trzech prowincji: na północnym zachodzie – Niżu Środkowoeuropejskiego (31), na południu – Karpat Zachodnich z Podkarpaciem Zachodnim i Północnym (51) oraz w części centralnej Wyżyn Polskich (34).
22
Przeważająca część obszaru badań należy do makroregionów: Wyżyna Śląska (341.1) i Wyżyna Woźnicko-Wieluńska (341.2). Położenie obszaru, zarówno w obrębie jednostek wyżynnych, jak i niżu, warunkuje urozmaiconą rzeźbę terenu. Centralna część obszaru, należąca do makroregionu Wyżyny Śląskiej, cechuje się rzeźbą strukturalną, silnie nawiązującą do budowy geologicznej podłoża. Liczne zręby tworzące płaskowyże, progi strukturalne oraz kotliny zapadliskowe porozcinane są dolinami rzek: Brynicy oraz Białej i Czarnej Przemszy. Wysokości bezwzględne kształtują się w przedziale od 250 do 300 m n.p.m., malejąc gwałtownie w kierunku Wyżyny Krakowsko-Częstochowskiej (341.3). Charakterystycznym elementem w morfologii jest Pustynia Błędowska, o powierzchni około 33 km2, będąca największym w Polsce obszarem piasków lotnych. Obszar Wyżyny Katowickiej (341.13), niemal w całości pokrywający się z obszarem konurbacji katowickiej, pierwotnie miał charakter zrębowy. Liczne płaskowyże rozdzielone były wypełnionymi materiałem czwartorzędowym kotlinami. Wschodnia część obszaru badań, obejmująca obszar Wyżyny Olkuskiej (341.32) i Wyżyny Częstochowskiej (341.31), wznosi się na wysokość od 350 do około 500 m n.p.m. Obszar ten zbudowany jest przeważnie z wapieni jurajskich, w których obserwować można dobrze rozwinięte zjawiska krasowe (ryc. 4.2.).
Zgodnie z regionalizacją Paczyńskiego (1995) praktycznie cały obszar zlewni Białej Przemszy i Przemszy należy do regionu śląsko-krakowskiego (XII). W obrębie regionu wydziela się trzy subregiony: triasu śląskiego (XII1), górnośląski (XII2) i jurajski (XII3). Dodatkowo w obrębie subregionu triasu śląskiego, na dokumentowanym obszarze, wydzielić można następujące rejony: gliwicki (XII1B), bytomski (XII1C) oraz chrzanowski (XII1D). Tylko niewielki, wschodni kawałek zlewni Białej Przemszy, należy do regionu nidziańskiego (XI).
23
Ryc. 4.2. Ukształtowanie powierzchni terenu obszaru badań (Źródło: NMT, CODGiK)
Dominującym sposobem zagospodarowania terenu według stanu na 2006 rok są obszary rolne (CORINE Land Cover). Zajmują one około 41% powierzchni omawianego obszaru, 29% obszaru stanowią tereny leśne. Pozostałe 30% obszaru testowego pokryte jest zwartą zabudową miast, rozległymi strefami przemysłowymi oraz kopalniami odkrywkowymi. Tereny te charakteryzują się znacznym stopniem przekształcenia naturalnego krajobrazu i zmienioną w wysokim stopniu morfologią.
4.1. Geologia obszaru
Obszar testowy charakteryzuje się skomplikowaną budową geologiczną. Niemal w całości położony jest w obrębie Bloku Górnośląskiego. Niewielki, południowy fragment obszaru położony jest w obrębie Bloku Małopolskiego (Żelaźniewicz 2011).
Osady podłoża paleozoicznego występują na całym omawianym obszarze na różnych głębokościach. Przeważnie są dolomity i wapienie dewonu, przykryte miąższym kompleksem kulmu lub wapienia muszlowego (karbon). W okolicy Jaroszewic i Wolbromia, występujące bezpośrednio pod utworami czwartorzędowymi. Utwory
[m n.p.m.]
24
permu, wyłącznie w facjach lądowych, odsłaniają się w części centralnej obszaru badań i zapadają w kierunku północno-wschodnim pod utwory triasu i jury, aż po strefę uskokową Kraków-Lubliniec. Wyżej ległe utwory mezozoiczne występują praktycznie na całym obszarze badań, z wyjątkiem Zapadliska Górnośląskiego. Są to piaskowce oraz wapienie triasu, iły, łupki, piaskowce jury oraz piaski i margle kredy. Osady te we wschodniej części omawianego obszaru tworzą Monoklinę Śląsko-Krakowską.
Tektonika obszaru ma charakter zrębowo-fałdowy, co wpływa na zróżnicowane głębokości występowania utworów, zróżnicowane miąższości i duże deniwelacje podłoża.
Osady czwartorzędowe występują powszechnie na analizowanym obszarze z wyjątkiem części wschodniej, gdzie na powierzchni odsłaniają się utwory triasu i jury.
Miąższość ich jest bardzo zróżnicowana i zmienia się od kilku metrów do ok. 60 metrów w dolinie kopalnej Przemszy (Górnik 2005a) oraz niemal 100 metrów w rejonie rowu Chrzanów-Dąb (Gajowiec 2005; Górnik 2006a, 2006b). W dolinach rzek: Małej Panwi, Brynicy, Warty, Mitręgi i Bobrka miąższość utworów czwartorzędu waha się w przedziale od 5 do 20 metrów (Górnik 2005b; Krawczyk, Zdechlik 2006).
Na obszarze Pustyni Błędowskiej, na pograniczu Wyżyny Śląskiej i Wyżyny Olkuskiej, w okolicach Kluczy, miąższość utworów czwartorzędu osiąga 50 metrów, w części dystalnej zaś 20 metrów (Gorczyca, Zdechlik 2006a; Gorczyca, Zdechlik 2006b).
Zwiększone miąższości, do około 20-25 metrów, obserwuje się również w rejonie rzeki Sztoły, na terenie obecnych kopalni piasków podsadzkowych DB Schenker Rail Polska Szczakowa S.A. i CTL Maczki Bór S.A. (Gajowiec, Siemiński 1997). Pod względem litologicznym wśród utworów plejstoceńskich dominują piaski, żwiry fluwioglacjalne oraz gliny pylaste i zwałowe. Holocen wykształcony jest w postaci osadów akumulacji rzecznej: piasków, żwirów i namułów.
4.2. Model konceptualny obszaru badań
Obszar badan charakteryzuje wysoki stopień komplikacji warunków hydrostrukturalnych oraz hydrodynamicznych. Dlatego też w celu symulacji warunków występowania użytkowych poziomów wodonośnych wód słodkich (w tym przypadku do ok. 200 m p.p.t.) i ustalenia ich zasobów, z podziałem na poszczególne piętra wodonośne, system wodonośny odwzorowano za pomocą siedmiu warstw modelowych
25
(Rodzoch i in. 2012; Pazio-Urbanowicz i in. 2012). Wydzielone warstwy mają charakter ciągły, a zmienność litologiczna skał została przedstawiona za pomocą zróżnicowanego współczynnika filtracji. Szczegółowy sposób wydzielenia warstw nie jest przedmiotem niniejszej pracy i nie będzie przedstawiany. Ze względu na zakres niniejszej pracy przedstawiona zostanie charakterystyka pierwszej warstwy modelowej. Obejmuje ona czwartorzędowy poziom wodonośny, a w przypadku wychodni skał starszych, ich stropową zawodnioną część. W tej warstwie powinny się znaleźć wszystkie utwory czwartorzędowe. Spąg warstwy powinien znajdować się na głębokości od kilku do 20 metrów p.p.t., a w dolinach rzek, dolinach kopalnych (pradolina Przemszy) oraz lokalnie, w obrębie Pustynie Błędowskiej, na głębokości do 60 metrów p.p.t.
Powierzchnia stropowa powinna być stosunkowo wygładzona i pozbawiona ostrych krawędzi, gdyż utwory czwartorzędowe zalegają na częściowo zdenudowanym podłożu.
26
5. Dane
W pracy wykorzystano dane dotyczące profili otworów wiertniczych zgromadzone w CBDH HYDRO (PIG-PIB) oraz dane wysokościowe pochodzące z numerycznego modelu terenu omawianego obszaru. Ze względu na charakter pracy, w tym brak konieczności dokładnego odwzorowania przebiegu struktur geologicznych, zbioru informacji geologicznych nie uzupełniano o profile wierceń ze Szczegółowej mapy geologicznej Polski, Mapy hydrogeologicznej Polski czy innych map i archiwaliów.
5.1. Dane z Centralnego Banku Danych Hydrogeologicznych HYDRO Centralny Bank Danych Hydrogeologicznych HYDRO jest bazą danych hydrogeologicznych prowadzoną przez Państwową Służbę Hydrogeologiczną (PSH).
Gromadzone są w niej między innymi podstawowe dane o odwiertach i ujęciach wód podziemnych z obszaru Polski. Choć początki bazy sięgają lat siedemdziesiątych ubiegłego stulecia, to dopiero w 2000 roku ostatecznie wdrożono system wspólnej platformy dla Centralnego Banku Danych Hydrogeologicznych (CBDH) i sieci Regionalnych Banków Danych Hydrogeologicznych (RBDH), które do tej pory niezależnie gromadziły dane. Następnie system zintegrowano z bazą danych GIS Mapy hydrogeologicznej Polski w skali 1 : 50 000 oraz Mapy geologiczno-gospodarczej Polski w skali 1 : 50 000. Aktualnie (kwiecień 2014) baza danych zawiera informacje o około 130 tysiącach obiektów hydrogeologicznych. Aktualizacja bazy danych odbywa się w trybie ciągłym.
Na informacje o obiekcie hydrogeologicznym, dostępne w Banku HYDRO, składają się m.in: dane geologiczne (profil geologiczny otworu), dane hydrogeologiczne (informacje charakteryzujące poziom wodonośny, w tym wydajności otworu, wyniki badań hydrogeologicznych, analizy fizykochemiczne wody, analizy izotopowe itp.) oraz dane techniczne otworu (położenie, głębokość otworu, przeznaczenie otworu itp.).
Dostęp do danych Banku HYDRO zapewniony jest poprzez platformę Państwowej Służby Hydrogeologicznej znajdującą się pod adresem http://spdpsh.pgi.gov.pl/PSHv7/. Podstawowa funkcjonalność systemu przetwarzania i dostępu do informacji Skarbu Państwa (SPD PSH) daje możliwość przeglądania, na warunkach nieodpłatnego dostępu, zasobu bazy. Funkcjonalność rozszerzona pozwala
27
na uzyskiwanie bezpośrednio za pomoc serwisu ściśle oznaczonych informacji hydrogeologicznych (Regulamin Banku HYDRO). Aktualnie (kwiecień 2014) wersja on-line Banku HYDRO umożliwia jedynie wyszukiwanie atrybutowe i przeglądanie w serwisie mapowym wybranych rodzajów danych. Pozostałe funkcjonalności, według zapewnień właściciela serwisu, będą dostępne od roku 2015 (www.psh.gov.pl). Na potrzeby komercyjne oraz naukowe dane z Banku HYDRO udostępniane są, w zależności od wymagań zamawiającego, jako wydruki i pliki .pdf (karty otworów) lub wydruki oraz pliki .xls zestawień tabelarycznych dowolnie zdefiniowanych atrybutów. Odchodzi się natomiast od udostępnia baz danych w formacie .mdb wraz z aplikacją Hydro2000plus v.Access pozwalającą na samodzielne pozyskanie danych.
Na potrzeby niniejszej pracy wykorzystano właśnie tę metodę pozyskania informacji o otworach hydrogeologicznych. W aplikacji Hydro2000plus (licencja HYDROEKO- BPiOW) utworzono zestaw zapytań SQL w wyniku którego otrzymano następujące informacje o otworach hydrogeologicznych:
- podstawowe dane identyfikacyjne obiektu hydrogeologicznego obejmujące jego nazwę, unikalny numer w bazie CBDH Bank HYDRO;
- informacje o położeniu obiektu (w tym województwo, powiat, gminę, miejscowość oraz ewentualny adres działki na której położony jest obiekt hydrogeologiczny;
współrzędne geograficzne obiektu w układzie 1942 oraz/lub WGS 84 i PUWG 1992 wraz z informacją w którym z tych układów pozyskano dane o położeniu oraz rzędną terenu;
- informacje techniczne o obiekcie w tym m.in. głębokość całkowitą otworu, przeznaczenie obiektu, stan obiektu, datę jego odwiercenia, okres i epokę spągu obiektu oraz informacje o profilu litostratygraficznym;
- informacje o numerze i nazwie ujęcia do którego przynależy obiekt.
Wiarygodność danych zgromadzonych w Banku HYDRO jest różna w zależności od kategorii danych. Za wiarygodne można uznać informacje dotyczące poziomu wodonośnego ujętego poprzez dany obiekt hydrogeologiczny, ponieważ pochodzą one bezpośrednio z badań hydrogeologicznych prowadzonych podczas odwiercania otworu i tuż po jego zafiltrowaniu. Parametry utworów nadległych, oparte są wyłącznie na makroskopowym opisie litologicznym przeprowadzonym podczas wiercenia otworu oraz właściwościach analogicznych utworów i ogólnie cechują się mniejszą wiarygodnością (Herbich i in. 2013).
28
5.2. Dane z Centralnego Ośrodka Dokumentacji Geodezyjnej i Kartograficznej
Dane o wysokości terenu pobrano z Numerycznego Modelu Terenu (NMT) udostępnionego przez Centralny Ośrodek Dokumentacji Geodezyjnej i Kartograficznej (CODGiK) w Warszawie w formacie ESRI TIN. NMT jest numeryczną, dyskretną reprezentacją wysokości terenu tworzoną poprzez ciągi współrzędnych x, y, z odpowiadającym rozproszonym punktom terenowym. Elementem składowym NMT jest algorytm interpolacyjny umożliwiający odtworzenie kształtu powierzchni terenu w danym obszarze. Model ten został opracowany na podstawie zdjęć lotniczych, skaningu lotniczego oraz map topograficznych. Wykonany został w układzie współrzędnych płaskich prostokątnych PUWG-1992, a dane wysokościowe odnoszą się do wysokości normalnych Kronsztadt 86. Pliki w formacie ESRI TIN zawierają rozproszone punkty wysokościowe tworzące nieregularną siatkę trójkątów, utworzone na podstawie danych pomiarowych Bazy Danych Topograficznych (BDT) – ASCII TBD. Zasób pierwotnie składa się z dziewięciu plików tekstowych zestawionych w tabeli 5.1.
Tab. 5.1. Zawartość plików wchodzących w skład zasobu NMT CODGiK
L.p. Oznaczenie Opis
1 P Punkty siatki
2 J Linie położone wokół obszarów planarnych (np. jezior)
3 C Cieki
4 K Punkty wysokościowe
5 O Obiekty inżynieryjne
6 Pz Punkty na obszarach wydzielonych 7 S Linie nieciągłości
8 Sz Linie nieciągłości w obszarach wydzieleń
9 z Punkty położone wokół tzw. obszarów wyłączonych
Szczegółowe parametry NMT opracowywanego na potrzeby TBD przestawiono w Wytycznych technicznych TBD – wersja 1.0 (uzupełniona). Specyfikacja danych zasobu podstawowego TBD. Na obszarach, gdzie nie było możliwe wykonanie pomiarów fotogrametrycznych dopuszcza się uzupełnienie zbioru danych o dane wysokościowe zawarte w dostępnych mapach topograficznych w skali 1 : 10 000 (warstwice, koty wysokościowe).
Średni błąd NMT w formacie ESRI TIN zawiera się w przedziale od 0,8 do 2,0 metrów, a aktualność danych to rok 2009 (codgik.gov.pl).
29
6. Metody
Analizy przeprowadzone w pracy, obejmowały kilka etapów: eksploracyjną analizę danych, analizę strukturalną i modelowanie wariogramu oraz interpolację.
W niniejszym rozdziale opisano zastosowane metody.
Wstępne przetworzenie danych przeprowadzono z wykorzystaniem oprogramowania Arc GIS 10.1 (Esri). Eksploracyjną analizę danych przeprowadzono w środowisku R (www.r-project.org) oraz w programie Microsoft Excel. Modelowanie wariogramów oraz interpolacje wykonano w programie Surfer 12 (Golden Software).
6.1. Wstępne przetworzenie danych
Do konstrukcji powierzchni modelowej zgodnie z przyjętymi w rozdziale 4.2 założeniami użyteczne będą tylko te otwory, w których przewiercono utwory czwartorzędu. Dlatego też, w pierwszym etapie pracy, ze wszystkich zgromadzonych otworów wyselekcjonowano otwory przewiercające utwory czwartorzędu. Rozkład przestrzenny wybranych otworów przedstawiono na rycinie 7.1. Następnie mając na uwadze fakt, że informacje pochodzące z Banku HYDRO obarczone są często błędną lokalizacją, dokonano weryfikacji lokalizacji. Początkowo porównano rzędne terenu otworów uzyskane z NMT wraz z rzędnymi znajdującymi się w Banku HYDRO.
W przypadku, gdy różnica wysokości wynosiła poniżej jednego metra założono, że otwór jest poprawnie zlokalizowany. W pozostałych przypadkach lokalizację weryfikowano na podstawię map topograficznych i Mapy hydrogeologicznej Polski 1 : 50 000 wykorzystując informacje o właścicielu i użytkowniku obiektu oraz danych technicznych tj. głębokości czy wydajności studni. W razie konieczności dokonywano korekty położenia obiektu i przypisywano rzędną terenu zgodną z NMT. Obiekty w przypadku których określenie prawidłowej lokalizacji nie było możliwe oraz obiekty zduplikowane usunięto ze zbioru danych.
Schemat postępowania przy wstępnym przetwarzaniu danych, na potrzeby niniejszej pracy, przestawiono na rycinie 6.1.
30
Ryc. 6.1. Schemat postępowania – wstępne przygotowanie danych
Podczas weryfikacji danych szczególną uwagę poświęcono obiektom, w których położenie spągu utworów czwartorzędowych odbiegało wyraźnie od obiektów otaczających. Analiza rozmieszczenia takich obiektów pozwoliła na znalezienie tych, których lokalizacja była błędna oraz na określenie ogólnego oraz lokalnych trendów danych. W szczególności otwory te, po dodatkowej kontroli głębokości przewidywanej powierzchni spągowej utworów czwartorzędowych na podstawie Szczegółowych map geologicznych Polski w skali 1 : 50 000, jeśli nie znaleziono wystarczających przesłanek do poprawienia lokalizacji, usunięto ze zbioru danych. Usunięcie wartości anomalnych było możliwe ze względu na bardzo dobre rozpoznanie geologiczne obszaru. Założono, że wszelkie wartości anomalne pojawiające się do tej pory w otworach i dotyczące występowania utworów czwartorzędowych byłyby wyjaśnione na Szczegółowej mapie geologicznej Polski lub podczas tworzenia map z serii Mapy hydrogeologicznej Polski 1 : 50 000, pierwszy poziom wodonośny, występowanie i hydrodynamika, skupiającej się w dużej mierze na występowaniu utworów czwartorzędowych.
Uzyskanie zbioru punktów przygotowanego do przeprowadzenia analizy statystycznej i interpolacji Usunięcie otworów o błędnej lub niezweryfikowanej lokalizacji
Wstępna weryfikacja lokalizacji otworów, wprowadzenie poprawek Wybór otworów, które przewiercają utwory czwartorzędu
Zgromadzenie dostępnych dla analizowanego obszaru danych otworowych z Banku HYDRO
31
Do przeprowadzenia powyższych analiz wykorzystano 16 arkuszy Szczegółowej mapy geologicznej Polski 1 : 50 000 oraz 14 arkuszy Mapy hydrogeologicznej Polski 1 : 50 000. Zestawienie wykorzystanych arkuszy zamieszczono w tabelach 6.1 oraz 6.2.
Tab. 6.1. Arkusze Szczegółowej mapy geologicznej Polski 1 : 50 000 (SmgP) wykorzystane w pracy
Nazwa i numer
arkusza mapy Autor mapy
Rok wydania
mapy
Autor objaśnień do mapy
Rok wydania objaśnień do
mapy
Kalety (877) Wyczółkowski 1969 Wilamowski 2002
Koziegłowy (878) Wyczółkowski 1962 Romanek 2009
Żarki (879) Heliasz i in. 1964 Heliasz i in. 1994
Bytom (910) Biernat 1955 Wilamowski 2001
Wojkowice (911) Biernat 1957 Wilamowski 2009
Zawiercie (912) Kotlicki 1966 Kotlicki 1967
Ogrodzieniec (913) Bednarek, Kaziuk,
Zapaśnik 1978 Bednarek, Kaziuk,
Zapaśnik 1978
Wolbrom (914) Bukowy 1964 Bukowy 1964
Gliwice (941) Żero 1957 Haisig 2009
Zabrze (942) Wyczólkowski 1960 Wilamowski 2001
Katowice (943) Biernat, Krysowska 1970 Wilamowski 2001
Jaworzno (944) Kurek, Paszkowski,
Preidl 1990 Kurek, Paszkowski,
Preidl 1992
Olkusz (945) Kurek, Preidl 1990 Kurek, Preidl 1993
Skała (946) Płonczyński 1995 Płonczyński 2000
Oświęcim (970) Biernat, Krysowska 1958 Wilamowski 2002
Chrzanów (971) Żero 1956 Płonczyński 2005
Tab. 6.2. Arkusze Mapy hydrogeologicznej Polski 1 : 50 000, pierwszy poziom wodonośny, występowanie i hydrodynamika (MhP PPW-WH) wykorzystane w pracy
Nazwa i numer arkusza MhP
PPW-WH Autor mapy i objaśnień Rok wydania
Kalety (877) Kempa, Pękała 2005
Koziegłowy (878) Górnik 2005
Żarki (879) Krawczyk, Zdechlik 2006
Bytom (910) Kempa, Bielewicz 2005
Wojkowice (911) Górnik 2005
Zawiercie (912) Gorczyca, Zdechlik 2006
Ogrodzieniec (913) Gorczyca, Zdechlik 2006
Zabrze (942) Górnik 2006
32
Nazwa i numer arkusza MhP
PPW-WH Autor mapy i objaśnień Rok wydania
Katowice (943), Cudak, Razowska-Jaworek, Brodziński 2005
Jaworzno (944) Gajowiec, Siemiński 1997
Olkusz (945) Gorczyca, Gągulski 2005
Skała (946) Krawczyk, Gągulski 2006
Oświęcim (970) Górnik 2006
Chrzanów (971) Gajowiec 2005
6.2. Wybór próby punktów weryfikacyjnych
Spośród wszystkich obiektów o zweryfikowanej lokalizacji dokonano wyboru otworów, które mają posłużyć do oceny dokładności interpolacji powierzchni. Punkty te dalej nazywane będą punktami kontrolnymi. Wybór punktów był arbitralny, przy czym starano się zachować ich równomierny rozkład.
6.3. Ekploracyjna analiza danych
Przeprowadzenie eksploracyjnej analizy danych pozwala na charakterystykę rozkładu analizowanych zmiennych, w tym zbadanie czy w zbiorze danych znajdują się wartości odstające, określenie istnienia trendu w danych czy poznanie ich rozkładu statystycznego (Urbański 2010). Aby zrozumieć strukturę danych zwyczajowo przeprowadza się szereg analiz oraz obliczeń mających na celu klasyfikację danych, a więc ich grupowanie na podstawie określonych cech. Odpowiednia klasyfikacja danych powinna pozwolić na dokładną aproksymację powierzchni statystycznej oraz wskazanie charakterystycznych cech rozmieszczenia danych. Liczba wyznaczanych klas powinna pozwalać na precyzyjne odwzorowanie danych i ukazanie dominujących trendów zaznaczających się w ich zbiorze. Uzyskana w wyniku wizualizacji klas mapa powinna pokazywać charakterystyczne cechy rozmieszczenia zjawiska (Mucha 1994).
Ustalenia granic przedziałów klasowych można dokonać na podstawie graficznej lub matematycznej analizy zbioru danych.
W pracy wykorzystano jedną z metod graficznych - sporządzenie histogramu (wykresu częstości), który służy do analizy liczebności wartości zmiennych, pokazuje częstości występowania wartości lub ich zakresów, ich rozmieszczenie, stopień koncentracji wartości wokół wartości średniej oraz pozwala wstępnie scharakteryzować
33
rozkład. Wyznaczenie przedziałów o równej wartości opiera się na podziale różnicy najmniejszej i największej wartości danego zjawiska na przyjętą arbitralnie liczbą klas.
Dla cechy ciągłej można również zastosować szereg rozdzielczy z przedziałami klasowymi. Liczbę klas k wyznacza się na podstawie liczby obserwacji n (Mucha 1994) według poniższych wzorów:
𝑘 ≈ √𝑛 (6.1)
𝑘 ≈ 1 + 3,222𝑙𝑜𝑔𝑛 (6.2)
W kolejnym etapie eksploracyjnej analizy danych wyznaczono miary położenia:
średnią, medianę, modę oraz kwantyle. Parametry te charakteryzują średni lub typowy poziom wartości cechy, a wokół nich skupiają się wszystkie pozostałe wartości analizowanej cechy.
Kwantyle są wartościami cechy badanego zbioru, które dzielą go na określone części pod względem liczby jednostek. Kwantyl p-rzędu dzieli zbiór na dwie części zawierające p*100% i (1-p)*100% elementów tego zbioru. Kwartyl pierwszy (Q1) to wartość próby, która dzieli zbiór tak, że 75% prób ma od niej wartości większe. Kwartyl drugi (Q2) to inaczej mediana. Kwartyl trzeci (Q3) to wartość próby, która dzieli zbiór tak, że 25% prób ma od niej wartości nie mniejsze. Rozstęp kwartylowy to różnica pomiędzy trzecim a pierwszym kwartylem. Oznacza on przedział zmienności cechy, w której znajduje się 50% „środkowych” obserwacji.
Następnie scharakteryzowano zmienność danych wykorzystując wariancję, odchylenie standardowe oraz rozstęp – miary opisujące rozproszenia danego zjawiska.
Rozstęp jest to różnica pomiędzy wartością maksymalną a minimalną występującą w danym zbiorze obliczana ze wzoru:
𝑅 = 𝑥𝑚𝑎𝑥− 𝑥𝑚𝑖𝑛 (6.3)
Wariancja określa stopień rozproszenia wartości zmiennej losowej wokół wartości średniej (oczekiwanej). Cecha ta jest wrażliwa na występowanie wartości anomalnych. Oblicza się ją ze wzoru:
𝜎2 =𝑛1∑𝑛𝑖=1(𝑥𝑖 − 𝑥̅)2 (6.4)
gdzie: xi - wartość i-tego pomiaru w próbie (i=1,2…,n) 𝑥̅ – wartość średniej arytmetycznej
