ZASTOSOWANIE MAPY WEKTOROWEJ POLSKI W SKALI 1:50 000 (VMAPL2) DO BUDOWY MODELI PRZEPŁYWU WÓD PODZIEMNYCH

(1)

ZASTOSOWANIE MAPY WEKTOROWEJ POLSKI W SKALI 1:50 000 (VMAPL2) DO BUDOWY MODELI PRZEP£YWU WÓD PODZIEMNYCH

APPLICATION OF THE VECTOR MAP OF POLAND 1:50 000 (VMAPL2) TO THE DEVELOPMENT OF NUMERICAL GROUNDWATER FLOW MODEL

DARIUSZKASZTELAN¹, MAREKMARCINIAK²

Abstrakt. W trakcie budowy matematycznych modeli przep³ywu wód podziemnych jednym z najwa¿niejszych zadañ jest dok³adne od- wzorowanie ukszta³towania powierzchni modelowanego obszaru. Dla terytorium Polski nie opracowano cyfrowego modelu terenu DTM w du¿ej skali. Dotychczas dla potrzeb modelowania przep³ywu wód podziemnych przeprowadzano schematyzacjê powierzchni terenu z mapy topograficznej na siatkê dyskretyzacyjn¹ modelu. Inn¹ drog¹ do uzyskania modelu rzeŸby jest digitalizacja rastrowych map topograficznych. W artykule omówiono procedurê konwersji warstw informacyjnych mapy wektorowej poziomu drugiego VMAPL2 do bazy danych programu Surfer. Pokazany zosta³ tak¿e wp³yw kroku siatki dyskretyzacyjnej na dok³adnoœæ odwzorowania ukszta³towania powierzchni terenu.

S³owa kluczowe: modelowanie matematyczne, cyfrowy model powierzchni terenu, mapa wektorowa.

Abstract. During the development of numerical groundwater flow models the accurate surface terrain reproduction is one of the most important tasks. No large-scale digital terrain models (DTM) of Poland have been developed yet. The numerical groundwater flow models were hitherto based on a digital terrain model schematically transformed from topographic maps. The other way to obtain DTM was digitizing topographic maps. This paper shows the procedure of thematic layer conversion from a vector map (VMAPL2) to a database of Surfer soft- ware. Moreover, the influence of different discretization network steps on the accuracy of the surface terrain reproduction is presented.

Key words: numerical modelling, digital terrain model, vector map.

WSTÊP

Mapa wektorowa poziomu drugiego (VMAPL2) posiada pokrycie arkuszowe dla ca³ego kraju. Powsta³a na podstawie wojskowych map topograficznych w skali 1:50 000. Opraco- wana zosta³a w formacie VPF (Vector Product Format), któ- ry jest standardem NATO. Mapa sk³ada siê z 13 grup tema- tycznych, w obrêbie których znajduj¹ siê warstwy informacyjne. Geometria poszczególnych obiektów mapy zapisana jest w uk³adzie wspó³rzêdnych WGS-84. W sk³ad mapy wektorowej wchodz¹ obiekty punktowe, liniowe i obszaro-

we wraz z przypisanymi tym obiektom informacjami zapisa- nymi w bazie danych.

Do modelowania matematycznego mo¿na wykorzystaæ dwie grupy tematyczne: powierzchniê terenu oraz hydrogra- fiê. W grupach tych znajduj¹ siê informacje o ukszta³towaniu powierzchni terenu w formie odwzorowanych poziomic wraz z rzêdnymi, rzêdnych punktów charakterystycznych oraz rzêdnych wód powierzchniowych. Elementy hydrogra- ficzne to m.in. rzeki, rowy, kana³y, tereny podmok³e, jeziora.

1Uniwersytet im. Adama Mickiewicza, Wydzia³ Nauk Geograficznych i Geologicznych, Instytut Geologii, ul. Maków Polnych 16, 61-606 Poznañ 2Uniwersytet im. Adama Mickiewicza, Wydzia³ Nauk Geograficznych i Geologicznych, Instytut Geografii Fizycznej i Kszta³towania Œrodowiska

Przyrodniczego, ul. Dziêgielowa 27, 61-680 Poznañ; e-mail: mmarc@amu.edu.pl

(2)

PROCEDURA IMPORTU WARSTW INFORMACYJNYCH Z MAPY WEKTOROWEJ VMAPL2 DO PROGRAMU VISUAL MODFLOW

Format VPF, w którym zosta³a opracowana mapa, jest formatem typowo wojskowym i rzadko wykorzystywanym przez oprogramowanie typu GIS. Pierwszym krokiem umo-

¿liwiaj¹cym dalsze wykorzystanie VMAPL2 do badañ modelowych jest przekonwertowanie danych do innego formatu. Konwersjê tak¹ wykonano u¿ywaj¹c programu MapInfo.

Potrzebne w procesie modelowania warstwy informacyjne zosta³y przekonwertowane do formatu SHP. Jednoczeœnie zmianie uleg³ uk³ad wspó³rzêdnych – ze wspó³rzêdnych k¹towych WGS-84 do wspó³rzêdnych kartezjañskich w uk³adzie 92. Konwersja taka umo¿liwi³a swobodne ope- rowanie danymi w œrodowisku GIS oraz w programie Visual MODFLOW ver. 4.1, w którym realizowano model prze- p³ywu wód podziemnych (Zheng, 1996). Jednym z najtrud- niejszych i czasoch³onnych zadañ na etapie konstruowania modelu jest dok³adne odwzorowanie ukszta³towania powierzchni modelowanego obszaru (Harbaugh i in., 2000).

Cel ten mo¿na osi¹gn¹æ poprzez przypisanie ka¿demu blo- kowi obliczeniowemu modelu odpowiedniej wartoœci rzêd- nej terenu oraz utworzenie mapy izoliniowej za pomoc¹ od- powiedniego programu komputerowego.

Ten drugi sposób jest dziœ czêœciej wykorzystywany, gdy¿ jest bardziej dok³adny oraz umo¿liwia wizualizacjê uzyskanej powierzchni w formie mapy izoliniowej albo blokdiagramu. Dodatkowo dok³adnoœæ tak uzyskanej mapy mo¿na ³atwo sprawdziæ, nak³adaj¹c j¹ na podk³ad topogra- ficzny. Mapê izoliniow¹ mo¿na te¿ wykorzystaæ do zobrazo- wania wyników badañ modelowych. Sposobem na wierne odwzorowanie powierzchni terenu jest zdigitalizowanie poziomic i punktów wysokoœciowych z mapy topograficznej, a nastêpnie utworzenie cyfrowego modelu rzeŸby terenu

(DTM). Do tego celu pos³u¿yæ mo¿e dobrze znany geolo- gom program Surfer. Digitalizacjê obiektów z rastrowych podk³adów topograficznych oraz ich przetworzenie do po- staci wektorowej (przebieg cieków, kszta³t jezior) czy bazy danych (wspó³rzêdne studni, piezometrów, punktów kalibra- cyjnych) mo¿na wykonaæ w niemal ka¿dym programie GIS.

Mo¿e to byæ wspomniany program Surfer lub program Did- ger, którego przyk³ady zastosowania w badaniach hydrogeo- logicznych opisa³ Kasztelan (2007). Aby uzyskaæ cyfrowy model powierzchni terenu w programie Surfer, niezbêdna jest baza danych, która bêdzie zawieraæ wspó³rzêdne x, y oraz rzêdn¹ z. Baza taka umo¿liwia otrzymanie po przelicze- niach statystycznych tablicy z danymi, która nastêpnie mo¿e byæ zaimportowana do programu komputerowego realizu- j¹cego modelowanie przep³ywu wód podziemnych. Aby skróciæ czas ¿mudnego digitalizowania poziomic i punktów wysokoœciowych, do budowy modelu wykorzystano warstwy informacyjne mapy wektorowej VMAPL2.

Do programu Didger zaimportowane zosta³y pliki *.SHP zawieraj¹ce: rzêdn¹ punktów wysokoœciowych, rzêdn¹ wód powierzchniowych, przebieg poziomic wraz z ich rzêdn¹.

Nastêpnie dane te zosta³y wyeksportowane do pliku *.DAT, w którym uzyskano wspó³rzêdne x, y oraz rzêdn¹ z. Dla ar- kusza VMAPL2, który sk³ada siê z dwóch arkuszy mapy 1:50 000, otrzymano bazê 257 801 punktów z przypisan¹ rzêdn¹. Baza ta pos³u¿y³a do stworzenia w programie Surfer cyfrowego modelu rzeŸby terenu (DTM), zapisanego w formacie pliku *.GRD. Format ten jest rozpoznawalny przez programy do modelowania matematycznego przep³ywu wód podziemnych, takie jak Visual MODFLOW, Processing MODFLOW czy te¿ GMS (McDonald, Harbaugh, 1988).

WP£YW SIATKI DYSKRETYZACYJNEJ MODELU NA DOK£ADNOŒÆ ODWZOROWANIA MORFOLOGII TERENU

Dobór kroku siatki dyskretyzacyjnej ma wa¿ne znaczenie dla dok³adnoœci odwzorowania modelowanego procesu (Franke i in., 1987). W modelach korzystaj¹cych z algoryt- mów obliczeniowych metody ró¿nic skoñczonych siatka dyskretyzacyjna powinna byæ zagêszczona w obszarach o du¿ej dynamice modelowanego procesu (Harbaugh, McDonald, 1996). Natomiast w obszarach gdzie zmiany modelowanego procesu s¹ niewielkie, krok siatki dyskretyzacyjnej mo¿e byæ wiêkszy. Zmiana kroku siatki dyskretyzacyjnej ma tak¿e wp³yw na dok³adnoœæ odwzorowania kszta³tu powierzchni terenu. A kolei kszta³t powierzchni terenu ma znacz¹cy wp³yw na warunki obiegu wody. Dlatego konieczne okaza³o siê sprawdzenie, jaki wp³yw na dok³adnoœæ odwzorowania powierzchni importowanej z mapy wektorowej ma krok siatki dyskretyzacyjnej modelu (Reilly, 2001).

Badania testowe przeprowadzono na przyk³adzie modelu filtracji wód podziemnych w rynnie gry¿yñsko-grabiñskiej po³o¿onej na Ziemi Lubuskiej. RzeŸba powierzchni tego obszaru jest bardzo urozmaicona (fig.1). G³ównym ciekiem od- wadniaj¹cym rynnê jest Gry¿ynka, która p³ynie z pó³nocy na po³udnie i uchodzi do Odry. Zlewnia Gry¿ynki w czêœci pó³nocnej i œrodkowej obejmuje obszar rynny gry¿yñsko-gra- biñskiej, natomiast w czêœci po³udniowej Pradolinê Warszaw- sko-Berliñsk¹. Po pó³nocno-wschodniej oraz pó³nocno-za- chodniej stronie zlewni Gry¿ynki wystêpuj¹ dwa obszary bezodp³ywowe, które zasilaj¹ rzekê sp³ywem podziemnym.

Rynna jest usytuowana po³udnikowo i ma 12 km d³ugoœci oraz oko³o 3 km szerokoœci. W œrodku rynny wystêpuj¹ dwa po³udnikowo usytuowane wa³y ozowe, które dziel¹ rynnê na trzy doliny: zachodni¹, œrodkow¹ i wschodni¹. Tworzy to ra-

(3)

zem bardzo urozmaicon¹ rzeŸbê terenu. Deniwelacje pomiê- dzy dnem dolin po³o¿onym na rzêdnych oko³o 50 m n.p.m.

a wysoczyzn¹ o rzêdnych oko³o 90 m n.p.m dochodz¹ do 40 m. RzeŸba obszarów bezodp³ywowych jest znacznie mniej urozmaicona, od oko³o 90 do oko³o 110 m n.p.m., maksy- malnie do 130 m n.p.m.

Odwzorowanie powierzchni terenu przeprowadzono na siatkach dyskretyzacyjnych o wymiarach: 200 × 200, 100 × 100, 50 × 50 oraz 25 × 25 m. Zazwyczaj w tym prze- dziale wymiarowym bloków obliczeniowych siatki dyskretyzacyjnej poszukuje siê kompromisu pomiêdzy niezbêdn¹ dok³adnoœci¹ odwzorowania modelowanego procesu a mo-

¿liwoœciami obliczeniowymi komputera.

Ocenê dok³adnoœci odwzorowania powierzchni terenu przeprowadzono na podstawie analizy dwóch przekrojów testowych. Pierwszy z nich (I – fig. 1) poprowadzono przez œro- dek rynny gry¿yñsko-grabiñskiej. Na przekroju tym wystê- puje szeœæ stoków o du¿ym nachyleniu, co pozwala dobrze oceniæ jakoœæ odwzorowania powierzchni terenu na siatkach dyskretyzacyjnych o ró¿nym kroku. Na figurze 2 pokazano przekrój testowy przy dyskretyzacji powierzchni siatkami

dyskretyzacyjnymi kwadratowymi o bokach: 200, 100, 50 oraz 25 m na tle zarysu rzeczywistego przekroju rynny.

Odchy³ki odwzorowania rzêdnych terenu obliczono ze wzoru:

D =ZT -ZM [1]

gdzie:

D – odchy³ka [m],

ZT – rzêdna terenu [m n.p.m.],

ZM– rzêdna terenu na modelu [m n.p.m.].

Wykresy przestrzennego rozk³adu odchy³ek dla ró¿nych kroków siatki dyskretyzacyjnej pokazano na figurze 3.

Drugi przekrój testowy (II – fig. 1) poprowadzono przez pó³nocny fragment wschodniego obszaru bezodp³ywowego.

Na przekroju tym ukszta³towanie powierzchni jest znacznie

³agodniejsze i tylko incydentalnie pojawia siê fragment wa³u kemowego o rzêdnych oko³o 130 m n.p.m. Na figurze 4 pokazano ten przekrój przy dyskretyzacji powierzchni siatkami dyskretyzacyjnymi kwadratowymi o bokach: 200, 100, 50 oraz 25 m na tle zarysu rzeczywistego przekroju rynny. Wy- kresy przestrzennego rozk³adu odchy³ek dla ró¿nych kroków siatki dyskretyzacyjnej przedstawiono na figurze 5.

Jak widaæ na figurach 3 i 5, odchy³ki maj¹ znak zarówno dodatni, jak i ujemny. Dlatego wartoœæ œredni¹ odchy³kiDœr

dla poszczególnych kroków siatki dyskretyzacyjnej obliczono jako:

[2]

D_œr ⁱ ^Ti ^Mi

n

Z Z

= n

å

-

gdzie:

n – liczba bloków obliczeniowych na przekroju testowym.

Tak¿e mediany obliczono na podstawie wartoœci bez- wzglêdnych odchy³ek dla poszczególnych kroków siatki dyskretyzacyjnej. W tabelach 1 i 2 zestawiono najwa¿niejsze elementy charakterystyki statystycznej. Elementy statystyki zbiorów odchy³ek pokazano równie¿ na wykresach (fig. 6).

Zmiennoœæ powierzchni terenu rynny gry¿yñsko-grabiñ- skiej by³a bardzo du¿a. Przy siatce dyskretyzacyjnej o boku 200 m maksymalne b³êdy odwzorowania powierzchni na obszarze rynny dochodzi³y do 12 m. Natomiast na obszarze wysoczyzny maksymalne b³êdy nieznacznie przekroczy³y 4 m. Wartoœæ œrednia odchy³ek w rynnie wynios³a oko³o 3 m, natomiast na wysoczyŸnie by³a o metr mniejsza. Media- na zbioru odchy³ek w ca³ym badanym obszarze by³a zbli¿o- na i wynios³a oko³o 1,5 m.

Kryteria oceny jakoœci odwzorowania powierzchni terenu powinny ka¿dorazowo uwzglêdniaæ nie tylko cel i skalê modelowania, ale tak¿e zró¿nicowanie powierzchni. Wa¿ne znaczenie ma lokalizacja i d³ugoœæ przekrojów testowych.

Przekroje testowe powinny byæ sytuowane zarówno w pod- obszarach o bardzo zró¿nicowanej rzeŸbie, jak i w podob- Fig. 1. Zlewnia Gry¿ynki wraz z przyleg³ymi

obszarami bezodp³ywowymi

Gry¿ynka river basin with surrounding blind drainage areas

(4)

Fig. 2. Dok³adnoœæ odwzorowania powierzchni terenu na przekroju testowym przez rynnê

The accuracy of surface terrain reproduction presented on the cross-section test through a postglacial channel

(5)

Fig. 3. Przestrzenny rozk³ad odchy³ek dla ró¿nych kroków siatki dyskretyzacyjnej dla przekroju testowego przez rynnê

Spatial distribution of deviation of different discretization network steps presented on the cross-section test through a postglacial channel

(6)

Fig. 4. Dok³adnoœæ odwzorowania powierzchni terenu na przekroju testowym przez wysoczyznê

The accuracy of surface terrain reproduction presented on the cross-section test through a till plain

(7)

Fig. 5. Przestrzenny rozk³ad odchy³ek dla ró¿nych kroków siatki dyskretyzacyjnej dla przekroju testowego przez wysoczyznê

Spatial distribution of deviation of different discretization network steps presented on the cross-section test through a till plain

(8)

szarach p³askich. Dla terenów o bardzo zró¿nicowanej rzeŸ- bie mo¿na przeprowadziæ kilka przekrojów testowych. D³u- goœæ ka¿dego przekroju powinna byæ tak dobrana, a¿eby przekrój oddawa³ lokalne cechy ukszta³towania powierzchni (Okoñska i in., 2004).

Dla badanego obszaru przyjêto, ¿e dok³adnoœæ odwzorowania powierzchni powinna spe³niaæ nastêpuj¹ce kryteria:

maksymalny b³¹d odwzorowania nie przekracza 2,0 m, œred-

nia i mediana zbioru odchy³ek s¹ mniejsze ni¿ 0,5 m. Tak sformu³owane kryteria na obszarze rynny spe³nia siatka dyskretyzacyjna o boku 50 m, natomiast na obszarze wysoczyzny siatka o boku 100 m ( tab. 1, 2, fig. 6).

Zmniejszenie kroku siatki dyskretyzacyjnej wi¹¿e siê ze znacznym wzrostem liczby bloków obliczeniowych, który tylko czêœciowo mo¿na ograniczyæ stosuj¹c siatkê dyskretyzacyjn¹ o zmiennym kroku.

Fig. 6. Wybrane statystyki zbiorów odchy³ek dla ró¿nych kroków siatki dyskretyzacyjnej

Selected statistics of deviation sets of different discretization network steps T a b e l a 1

Statystyki odchy³ek dla ró¿nych kroków siatki dyskretyzacyjnej dla przekroju przez rynnê Statistics of deviation of different discretization network steps presented on the cross-section test through a postglacial channel

Krok [m]

Minimum [m]

Maximum [m]

Œrednia [m]

Mediana [m]

25 –0,7039 0,8419 0,0986 0,0374

50 –1,8515 1,9932 0,3084 0,1317

100 –4,5888 5,6506 0,8941 0,5221

200 –11,9596 11,3424 2,7625 1,2925

T a b e l a 2 Statystyki odchy³ek dla ró¿nych kroków siatki dyskretyzacyjnej dla przekroju przez wysoczyznê Statistics of deviation of different discretization network steps

presented on the cross-section test through a till plain

Krok [m]

Minimum [m]

Maximum [m]

Œrednia [m]

Mediana [m]

25 –0,2570 0,2855 0,0437 0,0289

50 –0,5079 1,0122 0,1328 0,0848

100 –1,1413 1,3039 0,3678 0,2729

200 –3,7452 4,0518 1,7787 1,4854

(9)

MODEL PRZEP£YWU WÓD PODZIEMNYCH DLA RYNNY GRY¯YÑSKO-GRABIÑSKIEJ

Przy budowie modelu przep³ywu wód podziemnych dla rynny gry¿yñsko-grabiñskiej wykorzystano dwie warstwy informacyjne z VMAPL2: ukszta³towanie powierzchni oraz sieæ hydrograficzna. Wspó³rzêdne siatki modelu przyjêto zgodnie z siatk¹ kilometrow¹ uk³adu 92. Modelem objêto zlewniê Gry¿ynki oraz oba obszary bezodp³ywowe, ³¹cznie obszar o powierzchni 270 km². Po konwersji baza danych liczy³a 48 631 punktów wysokoœciowych, œrednia 180 punk- tów na km². Bazê punktów przeliczono w programie Surfer, stosuj¹c nastêpuj¹ce kroki siatki:

Dx = Dy = 200 m – 66 kolumn, 101 wierszy, czyli 6666 bloków;

Dx = Dy = 100 m – 131 kolumn, 201 wierszy, czyli 26 331 bloków;

Dx = Dy = 12,5 m – 1041 kolumn, 1601 wierszy, czyli 1 666 641 bloków.

Otrzymane mapy izoliniowe by³y nastêpnie porównywane z map¹ topograficzna w skali 1:50 000. Po analizie b³êdów odwzorowania na obu przekrojach testowych ostatecznie przyjêto siatkê dyskretyzacyjn¹ o zmiennym kroku (fig. 7).

Fig. 7. Siatka dyskretyzacyjna modelu przep³ywu wód podziemnych w rynnie gry¿yñsko-grabiñskiej

Discretization network of numerical groundwater flow model in the area of Gry¿ynka–Grabin postglacial channel

(10)

PODSUMOWANIE

Mapa wektorowa Polski w skali 1:50 000 mo¿e byæ bardzo efektywnym narzêdziem do zadawania ukszta³towania powierzchni terenu na modelach przep³ywu wód podziemnych. Z mapy tej mo¿na tak¿e wykorzystaæ inne warstwy informacyjne, np. sieæ hydrograficzn¹ czy rodzaj pokrycia powierzchni terenu. Dla zachowania wystarczaj¹co do- k³adnego odwzorowania zarówno powierzchni terenu, jak i innych elementów obiegu wody w modelowanym obsza-

rze konieczny jest w³aœciwy dobór kroku siatki dyskretyzacyjnej.

Przedstawion¹ w artykule procedurê mo¿na rozszerzyæ na tematyczne mapy wektorowe: hydrograficzn¹, sozolo- giczn¹, geologiczn¹ czy geoœrodowiskow¹. Zawsze jednak konieczne bêdzie sprawdzenie dok³adnoœci odwzorowania poszczególnych elementów œrodowiska przyrodniczego dla przyjêtego rodzaju siatki dyskretyzacyjnej.

LITERATURA

FRANKE O.L., REILLY T.E., BENNETT G.D., 1987 – Definition of boundary and initial conditions in the analysis of saturated ground-water flow systems – an introduction. U.S. Geological Survey Techniques of Water-Resources Investigations, Book 3, Chapter B5.

HARBAUGH A.W., BANTA E.R, HILL M.C., McDONALD M.G., 2000 – MODFLOW-2000, the U.S. geological survey modular.

Ground-water model – user guide to modularization. Concepts and the ground-water flow process. U.S. Geological Survey.

HARBAUGH A.W., McDONALD M.G., 1996 – User’s documen- tation for MODFLOW-96, an update to the USGS Modular Three-dimensional Finite-Difference Ground-Water Flow Mo- del. OFR 96-485.

KASZTELAN D., 2007 – Didger – zaawansowane narzêdzie do digitalizacji map. W: Informacja geograficzna w kszta³towaniu

i ochronie œrodowiska przyrodniczego (red. B. Medyñska-Gulij, L. Kaczmarek): 177–193. Bogucki Wyd. Nauk., Poznañ.

McDONALD M.G., HARBAUGH A.W., 1988 – A modular three- -dimensional finite-difference ground-water flow model.

USGS TWRI Chapter 6-A1.

OKOÑSKA M., KASZTELAN D., MARCINIAK M., 2004 – Uwa- runkowania dyskretyzacji obszaru dla modelu migracji znacz- ków przez próbkê gruntu. W: Hydrogeologia. Modelowanie przep³ywu wód podziemnych: 185–195. Wyd. UWroc., Wroc³aw.

REILLY T.E., 2001 – System and boundry conceptualization in ground-water flow simulation. USGS TWRI Book 3, Chapter B8.

ZHENG Ch., 1996 – MT3D99 a modular 3D multispecies transport simulator. User’s guide. S.S. Papadopulos & Associates, Inc.