Krajowa infrastruktura informacji przestrzennej jako źródło danych dla Systemu Informacji o Zlewni (RBIS)

ROCZNIKI GEOMATYKI 2013 m T XI m Z 1(58)

WP£YW ALGORYTMU OKREŒLANIA DRÓG

SP£YWU POWIERZCHNIOWEGO NA WYNIKI OCENY

ZAGRO¯ENIA GLEB EROZJ¥ WODN¥ W SKALI

ZLEWNI Z ZASTOSOWANIEM MODELU RUSLE*

THE INFLUENCE OF FLOW ROUTING ALGORITHM

ON THE RESULTS OF RUSLE-BASED

CATCHMENT-WIDE EROSION RISK ASSESSMENT

AGH Akademia Górniczo-Hutnicza w Krakowie, Wydzia³ Geodezji Górniczej i In¿ynierii Œrodowiska

S³owa kluczowe: erozja wodna gleb, RUSLE, modelowanie sp³ywu powierzchniowego, ocena wra¿liwoœci

Keywords: water erosion, RUSLE, flow routing, sensitivity assessment

Wprowadzenie

Dla oceny zagro¿enia gleb erozj¹ wodn¹ zaproponowano wiele zró¿nicowanych metod – od prostych podejœæ wskaŸnikowych po skomplikowane modele fizyczne. Do najbardziej znanych i najpowszechniej stosowanych nale¿¹ niew¹tpliwie model RUSLE (Revised Uni-versal Soil Loss Equation, Renard i in., 1997) wraz ze swoim poprzednikiem USLE (Univer-sal Soil Loss Equation, Wischmeier i Smith, 1978). Zosta³y one opracowane w celu szaco-wania strat gleby w skali pojedynczego pola, jednak po pewnych modyfikacjach stosowane s¹ równie¿ w skali zlewni lub regionu (Drzewiecki i in., 2013).

Model (R)USLE jest empirycznym równaniem matematycznym, w którym œredni roczny ubytek gleby szacowany jest na podstawie równania:

EA = R • K • LS • C • P (1)

gdzie:

E_A – masa gleby erodowanej z jednostki powierzchni w okreœlonym czasie [Mg×ha–1_×rok–1_]; R – czynnik erozyjnoœci opadu [MJ×ha–1_×mm×h–1_×rok–1_];

K – czynnik podatnoœci gleby na erozjê wodn¹ [Mg×h×MJ–1_×mm–1_];

LS – czynnik topograficzny – iloczyn czynników d³ugoœci (L) i nachylenia stoku (S) wymiarowy);

C – czynnik pokrycia i zarz¹dzania gruntami (bezwymiarowy); P – czynnik zabiegów przeciwerozyjnych (bezwymiarowy).

Pomimo powszechnoœci stosowania, modele USLE i RUSLE w niewielkim stopniu prze-badano w aspekcie ich wra¿liwoœci i propagacji niepewnoœci (Beven i Brazier, 2011). Doty-czy to zw³aszcza ich implementacji dla wiêkszych obszarów, realizowanych w œrodowisku systemów informacji geograficznej (GIS). Badania wra¿liwoœci modelu USLE w skali pola prowadzili np. Risse i in. (1993), wykazuj¹c i¿ najwiêkszy wp³yw na otrzymywane wyniki wywiera³y czynnik topografii (LS) oraz czynnik pokrycia i zarz¹dzania gruntami (C). Tet-zlaff i in. (2013) analizuj¹c wp³yw propagacji niepewnoœci danych wejœciowych na wyniki modelowania przeprowadzonego w œrodowisku GIS dla obszaru ca³ej Hesji (Niemcy) stwier-dzili, i¿ niepewnoœæ danych wejœciowych wywiera najwiêkszy wp³yw na okreœlenie warto-œci czynnika (LS), a w dalszej kolejnowarto-œci czynnika podatnowarto-œci erozyjnej gleb (K). Równie¿ Biesemans i in. (2000) stwierdzili w swoich badaniach, i¿ na niepewnoœæ koñcowego rezul-tatu modelowania przy u¿yciu modelu RUSLE w sposób decyduj¹cy wp³ywa niepewnoœæ okreœlenia czynnika topograficznego. Nale¿y podkreœliæ, i¿ w przypadku niektórych zastoso-wañ, zmiany wartoœci wskaŸników wchodz¹cych w sk³ad modelu, ograniczone nawet do stosunkowo niewielkich fragmentów modelowanego obszaru, mog¹ w znacz¹cy sposób wp³ywaæ na podejmowane w oparciu o wyniki modelowania decyzje (Drzewiecki i in., 2013). Niepewnoœæ okreœlenia wartoœci czynnika topograficznego nie jest wy³¹cznie rezultatem propagacji niepewnoœci danych wejœciowych, czyli niepewnoœci okreœlenia wysokoœci na numerycznym modelu terenu (NMT). W przypadku analiz prowadzonych przy u¿yciu GIS uzyskiwane wartoœci czynnika topograficznego zale¿eæ mog¹ równie¿ od algorytmu zasto-sowanego na etapie okreœlania dróg sp³ywu powierzchniowego (Mendicino, 1999; Drze-wiecki i Mularz, 2001). DrzeDrze-wiecki i Mularz (2008) oraz DrzeDrze-wiecki i in. (2008) stwierdzili, ¿e algorytmy zastosowane na etapie generowania œcie¿ek sp³ywu oraz sposób obliczania wartoœci czynnika topografii wp³ywaj¹ na mo¿liwoœæ skalibrowania modelu opartego na równaniu RUSLE w oparciu o obserwowane w ciekach ³adunki zawiesiny. Wp³yw wzorów przyjêtych do obliczania wartoœci czynnika topograficznego by³ w tym przypadku dominuj¹-cy, jednak wp³yw algorytmów generowania sp³ywu równie¿ by³ zauwa¿alny.

Celem badañ, prezentowanych w niniejszym artykule, by³o okreœlenie wp³ywu wyboru algorytmu generowania dróg sp³ywu powierzchniowego na wyniki oszacowania zagro¿enia gleb erozj¹ aktualn¹, przeprowadzonego w skali zlewni przy u¿yciu modelu RUSLE. Zagad-nienie to zosta³o zbadane zarówno w odniesieniu do wielkoœci, jak i rozk³adu przestrzennego szacowanych strat gleby.

Metodyka

Badania przeprowadzono na obszarze po³o¿onych na pó³noc od Krakowa zlewni Pr¹dnika i D³ubni. Obszar ten w czêœci zachodniej stanowi fragment Wy¿yny Olkuskiej nale¿¹cej do Wy¿yny Krakowsko-Czêstochowskiej, w czêœci wschodniej jest fragmentem Wy¿yny Mie-chowskiej, czêœæ po³udniowo-wschodnia natomiast nale¿y do P³askowy¿u Proszowickiego

– obie te jednostki wchodz¹ w sk³ad Wy¿yny Ma³opolskiej (Kondracki, 1988). Zarówno ze wzglêdu na wystêpuj¹c¹ w tym rejonie wykszta³con¹ na lessach pokrywê glebow¹, jak i ukszta³towanie powierzchni terenu, teren ten uwa¿any jest za potencjalnie zagro¿ony wystê-powaniem erozji wodnej gleb.

Obszar zlewni Pr¹dnika i D³ubni obejmuje tereny o charakterze wiejskim, podmiejskim i typowo miejskim (rys. 1). W czêœci po³udniowej znajduj¹ siê silnie zurbanizowane tereny miasta Krakowa. W œrodkowej czêœci obszaru badañ znajduje siê miasto Ska³a, najwiêksza poza Krakowem miejscowoœæ na analizowanym terenie. Na pó³noc od Ska³y mamy do czy-nienia z terenami typowo wiejskimi. Na obszarze pomiêdzy Krakowem a Ska³¹ krajobraz stopniowo ulega³ przekszta³ceniom, nabieraj¹c obecnie charakteru typowo podmiejskiego.

Dane

Modelowanie erozji wodnej gleb w œrodowisku systemów informacji geograficznej wy-maga zgromadzenia danych przestrzennych, umo¿liwiaj¹cych realizacjê wybranego modelu. W prezentowanych badaniach wykorzystano nastêpuj¹ce dane Ÿród³owe:

m Mapê czynnika erozyjnoœci opadu (Drzewiecki i in., 2013);

m Mapê gatunków gleb opracowan¹ na podstawie pozyskanej z IUNG cyfrowej wersji mapy glebowo-rolniczej w skali 1:25 000 (Kulas, 2011);

m Mapê u¿ytkowania terenu zlewni rzeki Pr¹dnik i D³ubnia o szczegó³owoœci tematycz-nej odpowiadaj¹cej czwartemu poziomowi szczegó³owoœci programu CORINE Land Cover, opracowan¹ na drodze fotointerpretacji ortofotomap lotniczych i satelitarnych (Drzewiecki, 2010);

m Numeryczny model terenu (NMT) opracowany poprzez rasteryzacjê modelu TIN, pochodz¹cego z Pañstwowego Zasobu Geodezyjno-Kartograficznego (Drzewiecki i in., 2013). Wykorzystano najlepszy dostêpny dla ca³oœci obszaru badañ model wyso-koœciowy (dok³adnoœæ wysokoœciowa w granicach 20-30 cm).

Wszystkie dane przetworzono, przekszta³caj¹c je do postaci map rastrowych o rozdziel-czoœci 30´30 m (1184´1241 komórek).

Sposób obliczania sk³adowych modelu USLE

W przeprowadzonych obliczeniach wykorzystano rozk³ad przestrzenny czynnika erozyj-noœci opadu opracowany dla obszaru Ma³opolski w oparciu o dane z lat 1981-2009 (Drze-wiecki i in., 2013). Wartoœci czynnika erozyjnoœci opadu oszacowane zosta³y w tym przy-padku dla punktów pomiarowych, na podstawie miesiêcznych sum opadów z lat 1981-2009, wed³ug wzoru zaproponowanego dla warunków polskich przez Licznara (2004).

Rozk³ad przestrzenny wartoœci czynnika podatnoœci gleb na erozjê wodn¹ uzyskano przy-pisuj¹c poszczególnym gatunkom gleb wartoœci czynnika zaproponowane dla gleb polskich, w oparciu o badania przeprowadzone w IUNG w Pu³awach (Wawer i in., 2005; Stuczyñski i in., 2010). W przypadku lessów i rêdzin zastosowano wartoœci zaproponowane przez Drze-wieckiego i in. (2013).

Na potrzeby okreœlenia wartoœci czynnika pokrycia i zarz¹dzania gruntami, wejœciow¹ mapê u¿ytkowania terenu zreklasyfikowano zgodnie z tabel¹ 1. Wartoœæ czynnika zabiegów przeciwerozyjnych przyjêto jako równ¹ 1, co jest to¿same z za³o¿eniem, ¿e zabiegi takie nie s¹ stosowane.

Wartoœæ czynnika topograficznego dla poszczególnych komórek rastra obliczono stosu-j¹c podejœcie zaproponowane przez Desmeta i Goversa (1996b). W metodzie tej d³ugoœæ stoku zastêpowana jest wielkoœci¹ jednostkowej powierzchni zasilania, co pozwala na lepsze modelowanie wp³ywu ukszta³towania terenu na zachowanie siê wody sp³ywaj¹cej po jego powierzchni. Je¿eli elementy zbocza reprezentowane s¹ poprzez komórki rastra, to jednost-kow¹ powierzchniê obszaru zasilania dla danej komórki otrzymamy dziel¹c pole powierzchni jej obszaru zasilania przez odleg³oœæ, jak¹ przep³ywaj¹ca woda przebywa przemieszczaj¹c siê wewn¹trz tej komórki (odleg³oœæ ta jest to¿sama z rozmiarem komórki tylko je¿eli sp³yw nastêpuje dok³adnie w kierunku pó³noc-po³udnie lub wschód-zachód).

Tabela 1. Wartoœci czynnika pokrycia i zarz¹dzania gruntami (C) przyjête dla poszczególnych klas u¿ytkowania C æ œ o tr a W Kalsau¿ytkowanai 0 zabudowazwatraweilkomeisjka,zabudowazwatrameisjka,zabudowaluznawielorodzinna , e w o ³ s y m e z r p y n e r e t , o g ei k sj ei m u p y t a n n i z d o r o l ei w a n z u l a w o d u b a z , o g e w o k o l b u p y t y r u t k u rt s a rf n i j e n z c y t si l a j c e p s y n e r e t , a w t ci n l o r h c y w o g u ³ s u i h c y w o ³ s y m e z r p ñ e z d ¹ z r u y n e r e t ,) e w o k sj o w . p n ( e n l a j c e p s y n e r e t ,j e n z ci l b u p i c œ o n z c e t y ¿ u y n e r e t , e w o l d n a h y n e r e t ,j e n z ci n h c e t y n e r e t i e j el o k ,¹ w o g o r d ¹ j c a k i n u m o k z e n a z ¹i w z y n e r e t i i g o r d , e w o tr o p s n a rt y z a b , u tl u k a c sj ei m -z c u t z s o a k si n t o l ,y ¿ a r a g y ³ o p s e z , e c ¹ z s y z r a w o t y n e r e t i a n j y c a k i n u m o k æ ei s a n o ¿ o ³ z , e w o j el o k ,y m o ³ o i n ei m a k , h c y t s al i i h c y w o h c u r k o w ó c w o r u s a k si b o r y w , h c y w o tr a t s g ó r d i n h c z r ei w a n j e n -a t s w ó d a p d o e w o ³ s y m e z r p a k si w o ³ a w z ,¹ i c œ o n n il œ o r e c ¹ j a t s a r a z i e n a w o w y tl u k e r a k si b o r y w -y s y w i a k si w o ³ a w z , e n l a n u m o k a k si p y s y w ,i k i n d a s o – h c y n n y ³ p w ó d a p d o a k si w o ³ a w z , h c y ³ ,y ³ a k s e t êi n o ³ s d o , e w o i n r al k z s y w a r p u , w ó n e s a b y ³ o p s e z , e c ¹ j a t s a r a z i e n a w o w y tl u k e r a k si p , e n b y r y w a t s , e w o r o p a z – e n j y c n e t e r i k i n r o i b z , a r o i z e j – e n d o w i k i n r o i b z e n l a r u t a n ,y ³ a n a k ,i k e z r ¹ d o w e n o i n ³ e p y w a k si l d a p a z i a k si b o r y w 1 0 0 , 0 cmentarze, alsy 2 0 0 , 0 zaelseinaiisamoseiw,ybagnaœród¹ldowe 5 0 0 , 0 zabudowaluŸnajednorodzinnatypumeisjkeigo,parkimeisjkeiiweisjkei,ogrodydydaktyczne, y b ê r z i ai n ei s el y w , e n œ el i k ³ ó k z s , ai n e z c a z r k a z i a k si w o s o z r w 1 0 , 0 zabudowajednorodzinnatypuweisjkeigo,zabudowa eltnsikowa,terenyspotrowe,³¹kizudzai -e n l a r u t a n a k si w t s a p i y w a r u m , w ó w e z r k i w e z r d m e ³ 2 0 , 0 zeielñce,skwer,y rtawnik,iterenyneiu¿ytkowanegospodarczo,terenywypoczynkowo-kem -w ó w e z r k i w e z r d u ³ ai z d u z e b i k ¹ ³ , e w o g n i p 5 0 , 0 ogrodydzai³kowe,sad,ypalntacjekrzewówowocowych,od³og,iterenyrolncizezprzewag¹³¹k ¹ g a w e z r p z e z ci n l o r y n e r e t , ñ ei w e z r k a z i ñ ei w e z r d a z ¹ g a w e z r p z e z ci n l o r y n e r e t , k si w t s a p i u ¿ o ³ d o p a n a n o z s o r p z o r æ œ o n n il œ o r , m y t s y z c z s ai p u ¿ o ³ d o p a n a n o z s o r p z o r æ œ o n n il œ o r , w ó g o ³ d o m y t si l a k s 1 , 0 z³o¿onesystemyuprawidzai³ekzdzai³kamima³ym,iterenyrolncizezeznacznymudzai³em³¹k, w ó g o ³ d o m e ³ ai z d u m y n z c a n z e z e z ci n l o r y n e r e t , ñ ei w e z r k a z i ñ ei w e z r d a z 5 1 , 0 z³o¿onesystemyuprawidzai³ekzdzai³kamidu¿ymiiœrednimi 2 , 0 gruntyorne 5 , 0 palcebudów 1 paiskipozasrtef¹brzegow¹morza

Algorytmy obliczania powierzchni zasilania

W ramach przeprowadzonych badañ przetestowano wp³yw siedmiu algorytmów okre-œlania dróg sp³ywu powierzchniowego.

Deterministic 8 (D8) (O'Callaghan i Mark, 1984)

Jest to jeden z najprostszych algorytmów s³u¿¹cych do wyznaczania kierunku sp³ywu wody na numerycznym modelu terenu. Kierunek sp³ywu jest wybierany w wyniku wylicze-nia ró¿nicy wysokoœci miêdzy rozpatrywan¹ komórk¹ oraz tymi, którymi jest otoczona. Sp³yw kieruje siê do komórki, dla której otrzymano najwiêksz¹ wartoœæ gradientu. Zasad-nicz¹ wad¹ algorytmu jest jego tendencja do generowania wielu równoleg³ych œcie¿ek sp³y-wu na zboczach o jednostajnym nachyleniu.

Random eight-neighbor (Rho8) (Fairfield i Leymarie, 1991)

Rho8 stanowi modyfikacjê algorytmu D8 polegaj¹c¹ na wprowadzeniu funkcji stocha-stycznej na etapie okreœlania gradientu. Rozwi¹zanie to pozwala na przezwyciê¿enie proble-mu równoleg³ych œcie¿ek sp³ywu. Element stochastyczny powoduje jednak, i¿ ka¿da realiza-cja algorytmu daje dla tych samych danych wejœciowych nieco odmienny rozk³ad prze-strzenny kierunków sp³ywu, co stanowi jego powa¿n¹ wadê. Dodatkowo w pewnych ob-szarach wprowadzaæ on mo¿e zafa³szowania (Jones, 2002). Podobnie jak w przypadku algorytmu D8 nie ma równie¿ mo¿liwoœci modelowania dyspersji strumienia sp³ywu. Kinematic Routing Algorithm (KRA) (Lea, 1992)

Sp³yw kierowany jest w dó³ stoku, wzd³u¿ œcie¿ki sk³adaj¹cej siê z prostoliniowych seg-mentów o kierunku zgodnym z kierunkiem ekspozycji poszczególnych komórek rastra. Ka¿-da komórka jest reprezentowana jako p³aszczyzna, która zosta³a wpasowana w punkty sta-nowi¹ce naro¿niki komórki. Ich wysokoœci s¹ szacowane poprzez interpolacjê wysokoœci œrodków s¹siednich komórek w siatce. Zalet¹ tej procedury jest okreœlanie kierunków sp³y-wu jako k¹tów od 0 do 2P. Podobnie jak w przypadku poprzednich algorytmów nie ma jednak mo¿liwoœci modelowania dyspersji.

Poniewa¿ do okreœlenia p³aszczyzny s¹ potrzebne tylko trzy punkty, Tarboton (1997) zauwa¿a, ¿e za³o¿enie p³aszczyzny dopasowanej lokalnie do ka¿dej komórki rastra wymaga przybli¿enia. Najlepiej dopasowana p³aszczyzna nie mo¿e na ogó³ przejœæ przez cztery punkty naro¿nikowe, co prowadzi do nieci¹g³ej reprezentacji powierzchni na krawêdziach komórek. P³aszczyzny dopasowane lokalnie do niektórych kombinacji punktów naro¿nikowych mog¹ prowadziæ do niespójnoœci lub przeciwnych do intuicyjnych kierunków przep³ywu, które s¹ problemem zarówno tej metody, jak i opisanego w dalszej czêœci algorytmu DEMON. Deterministic Infinity (D¥) (Tarboton, 1997)

Sp³yw z jednej komórki dzielony jest na dwie z oœmiu j¹ otaczaj¹cych. Wybór komórek odbywa siê poprzez formowanie oœmiu trójk¹tów, które tworzy siê przez po³¹czenie œrodka rozpatrywanej komórki ze œrodkami komórek j¹ otaczaj¹cych. Po utworzeniu trójk¹tów wy-brane zostaj¹ dwie komórki z najwiêkszym gradientem spadku. Sp³yw rozdzielany jest we-d³ug wzorów:

(2)

(3) gdzie:

P₁, P₂– udzia³ komórki w podziale sp³ywu;

a₁, a₂ – k¹t pomiêdzy odcinkiem ³¹cz¹cym œrodki komórek i wyznaczonym kierunkiem sp³ywu. Multiple Flow Direction (MFD) (Quinn i in., 1991)

Algorytmy FD8 oraz FRho8 czêsto opisywane s¹ pod wspóln¹ nazw¹ Multiple Flow Direction (MDF) i s¹ modyfikacjami algorytmów D8 oraz Rho8. W metodzie tej w przeci-wieñstwie do poprzednich, woda mo¿e byæ kierowana do wielu komórek. Czêœæ sp³ywu powierzchniowego przekazywana komórkom znajduj¹cym siê ni¿ej jest obliczana z wzoru:

(4) gdzie:

s_i – wielkoœæ nachylenia pomiêdzy wêz³ami centralnej komórki a s¹siaduj¹cymi komórkami; a – sta³a dodatnia (wartoœæ zalecana a=1,1).

Wad¹ tego algorytmu jest du¿e rozproszenie sp³ywu powierzchniowego. Triangular Multiple Flow Direction (MD¥) (Seibert i McGlynn, 2007)

Algorytm ten ³¹czy w sobie zalety algorytmów MFD i D¥. Podobnie jak w podejœciu Tarbotona (1997), do obliczenia lokalnych kierunków nachylenia oraz gradientów wokó³ badanej komórki, wykorzystuje siê osiem trójk¹tnych powierzchni, tworzonych przez po³¹-czenie œrodka rozpatrywanej komórki ze œrodkami dwóch s¹siednich. Dla ka¿dej z utworzo-nych w ten sposób p³aszczyzn oblicza siê kierunek przep³ywu przy najwiêkszym gradiencie. Po okreœleniu kierunków sp³ywu dla badanej komórki, nastêpuje rozdzielenie na nie jej obsza-ru zasilania. Powierzchnia obszaobsza-ru zasilania komórki (czyli obszaobsza-ru, z którego sp³ywa do niej woda po powierzchni terenu) dzielona jest pomiêdzy po³o¿one ni¿ej komórki, wed³ug proce-dury zaproponowanej przez Quinna (1991) (wzór 4).

W wielu przypadkach algorytmy MD¥ i D¥ daj¹ takie same wyniki (na przyk³ad na p³askich lub zbie¿nych stokach), ale rezultaty ró¿ni¹ siê, gdy wystêpuje wiêcej ni¿ jeden kierunek najwy¿szego lokalnego spadku z badanej komórki (na przyk³ad na rozbie¿nych stokach lub wzd³u¿ grzbietów).

DEMON (digital elevation model networks) (Costa-Cabral i Burges, 1994)

W przypadku tego algorytmu sp³yw jest inicjowany oddzielnie w ka¿dej komórce rastra, a nastêpnie kierowany jest wzd³u¿ spadku terenu do zag³êbienia terenu lub do krawêdzi NMT. Kierunek sp³ywu okreœlany jest wed³ug koncepcji przedstawionej przez Lea (1992), jednak modelowanie ma charakter dwuwymiarowy. Sp³yw kierowany jest w ca³oœci do pojedynczej po³o¿onej ni¿ej komórki tylko wtedy, jeœli okreœlony kierunek sp³ywu jest

wielo-    _D D_D 3      _D D_D 3 

¦

krotnoœci¹ 90 stopni. W innym przypadku nastêpuje podzia³ strumienia sp³ywu pomiêdzy dwie komórki oraz okreœlenie proporcji tego podzia³u. W efekcie dzia³ania algorytmu dla ka¿dej komórki rastra obliczana jest wartoœæ informuj¹ca, jaka czêœæ sp³ywu pochodz¹cego z komór-ki startowej do niej dociera. Obliczenia takomór-kie powtarzane s¹ dla œcie¿ek sp³ywu inicjowanych po kolei we wszystkich komórkach rastra. Wielkoœæ akumulacji sp³ywu dla pojedynczej komórki jest sum¹ wartoœci wykonanych dla wszystkich realizacji algorytmu. Ze wzglêdu na wyszuka-ne kodowanie algorytm ten jest rzadko realizowany w oprogramowaniu GIS.

Przebieg analizy

Pierwszym etapem przeprowadzonych badañ by³o wygenerowanie, w oparciu o NMT, map spadków i powierzchni obszarów zasilania, przy u¿yciu ka¿dego z wybranych do ana-lizy algorytmów generowania sp³ywu powierzchniowego. Na ich podstawie wygenerowano mapy czynnika topografii (LS). Obliczenia w tym zakresie zrealizowano wykorzystuj¹c opro-gramowanie SAGA GIS.

Uzyskane mapy zaimportowano do programu Idrisi, gdzie przeprowadzono dalsze anali-zy. Ka¿d¹ z map czynnika topograficznego przemno¿ono przez mapy pozosta³ych czynni-ków modelu RUSLE, uzyskuj¹c oszacowanie wielkoœci strat gleby w wyniku erozji aktual-nej. Wynikowe mapy przekszta³cono na mapy klas zagro¿enia erozyjnego wed³ug kryteriów przyjêtych w pracy Drzewieckiego i in. (2013).

Wyniki

W tabeli 2 przedstawiono wartoœci wspó³czynników korelacji (R2_{) pomiêdzy} wartoœcia-mi czynnika topograficznego (LS) okreœlonywartoœcia-mi z wykorzystaniem poszczególnych algoryt-mów. Analiza uzyskanych wyników pokazuje, i¿ najwiêkszym podobieñstwem cechuj¹ siê mapy czynnika LS utworzone przy u¿yciu algorytmów D¥, MD oraz MFD¥. Znaczne podobieñstwo wartoœci czynnika topografii wykazuj¹ równie¿ mapy powsta³e z wykorzy-staniem algorytmów KRA i DEMON. Wartoœæ R2 _{powy¿ej 0,9 uzyskano równie¿ dla} wyni-ków otrzymanych w oparciu o algorytmy D8 i MD¥. Najni¿szy poziom korelacji wartoœci czynnika topograficznego otrzymano dla map opracowanych z wykorzystaniem algoryt-mów D8 i DEMON lub KRA oraz Rho8 i DEMON lub KRA.

Analiza procentowego udzia³u klas zagro¿enia erozyjnego uzyskanych z wykorzystaniem poszczególnych map czynnika topograficznego (tab. 3) pokazuje, i¿ najwy¿sze zagro¿enie

Tabela 2. Wartoœci wspó³czynników korelacji (R2_{) pomiêdzy wartoœciami czynnika topografii}

wyznaczonymi z wykorzystaniem poszczególnych algorytmów

D¥ DEMON KRA MFD M ¥D Rho8 8 D 0,86 0,77 0,79 0,86 0,91 0,87 D¥ 0,82 0,85 0,97 0,96 0,83 N O M E D 0,95 0,81 0,83 0,75 A R K 0,84 0,85 0,77 D F M 0,96 0,83 D M ¥ 0,87

erozyjne na badanym obszarze przewidywane by³o przy wykorzystaniu algorytmów Rho8, DEMON oraz D8. W pozosta³ych przypadkach by³o ono wyraŸnie mniejsze. Podkreœliæ nale¿y, i¿ w ka¿dym z analizowanych przypadków prognozowany poziom zagro¿enia erozyj-nego uznaæ nale¿y za niski.

Porównuj¹c wp³yw stosowanych algorytmów na wyniki oszacowania erozji aktualnej, poszczególne mapy zagro¿enia erozyjnego zdecydowano siê poddaæ reklasyfikacji, tak aby dla ka¿dego wariantu zidentyfikowaæ obszary zagro¿one w stopniu co najmniej œrednim lub wy¿szym (rys. 2). Taki poziom zagro¿enia oznacza, i¿ w terenach tych nast¹piæ mo¿e trwa³a degradacja profilu glebowego (Józefaciuk i Józefaciuk, 1996). W celu porównania zgodno-œci ocen przeprowadzonych z wykorzystaniem poszczególnych algorytmów, dla uzyska-nych map obszarów zagro¿ouzyska-nych, obliczono wartoœci wspó³czynnika Kappa (Cohen, 1960). Najbardziej podobne do siebie oceny uzyskano stosuj¹c algorytmy: D¥, MD¥ oraz MFD (wsp. Kappa > 0,8). Najni¿sza zgodnoœæ wystêpuje dla par: MFD–DEMON, MFD–Rho8, oraz KRA–Rho8 (wsp. Kappa < 0,6).

Tabela 3. Zagro¿enie erozj¹ wodn¹ powierzchniow¹ na terenie zlewni rzek D³ubnia i Pr¹dnik

m t y r o g l a y t y ¿ U D8 D¥ DEMON Rho8 e n j y z o r e ei n e ¿ o r g a Z [ha] [%] [ha] [%] [ha] [%] [ha] [%] ij z o r e k a r b 36141,26 76,06 36437,76 76,69 35172,77 74,02 35693,57 75,12 a b a ³ s a j z o r e 10697,83 22,51 10529,62 22,16 11622,22 24,46 11068,74 23,30 a n a w o k r ai m u a j z o r e 593,03 1,25 493,07 1,04 626,06 1,32 655,81 1,38 ai n d e r œ a j z o r e 67,50 0,14 47,05 0,10 74,32 0,16 76,77 0,16 a n li s a jz o r e 14,87 0,03 7,83 0,02 19,40 0,04 19,13 0,04 a n li s o z d r a b a j z o r e 0,92 0,002 0,09 0,0002 0,65 0,001 1,40 0,003 m t y r o g l a y t y ¿ U KRA MFD M ¥D e n j y z o r e ei n e ¿ o r g a Z [ha] [%] [ha] [%] [ha] [%] ij z o r e k a r b 37084,79 78,05 36392,45 76,59 36524,90 76,87 a b a ³ s a j z o r e 9953,78 20,95 10599,66 22,31 10438,29 21,97 a n a w o k r ai m u a j z o r e 424,24 0,89 476,28 1,00 496,42 1,04 ai n d e r œ a jz o r e 43,43 0,09 40,75 0,09 47,45 0,10 a n li s a j z o r e 8,82 0,02 6,26 0,01 8,21 0,02 a n li s o z d r a b a j z o r e 0,36 0,0008 0,02 0,00005 0,14 0,0003

Dok³adniejsza analiza map terenów zagro¿onych erozj¹ pozwoli³a równie¿ na stwierdze-nie, i¿ obszary zagro¿enia prognozowane z wykorzystaniem algorytmu KRA prawie w ca³o-œci (w 98,6 % w zlewni Pr¹dnika i 97,6% w zlewni D³ubni) znalaz³y siê równie¿ w prognozie powsta³ej z u¿yciem algorytmu DEMON. Jednak¿e w przypadku zastosowania algorytmu DEMON powierzchnia terenów zagro¿onych erozj¹ by³a prawie dwukrotnie wiêksza.

Dyskusja wyników i wnioski

W przeprowadzonych badaniach porównano mapy czynnika topograficznego oraz mapy klas zagro¿enia gleb aktualn¹ erozj¹ wodn¹ wygenerowane z wykorzystaniem siedmiu algoryt-mów generowania sp³ywu powierzchniowego (D8, Rho8, KRA, D¥, MFD, MD¥ oraz DEMON). Celem badañ by³o poszukiwanie odpowiedzi na pytanie czy dobór algorytmu u¿yte-go na etapie okreœlania dróg sp³ywu powierzchnioweu¿yte-go wp³ywa na wyniki modelowania za-gro¿enia erozyjnego. Uzyskane wyniki pozwoli³y na stwierdzenie, i¿ wp³yw taki istnieje.

Analizuj¹c wartoœci czynnika topograficznego uzyskane w poszczególnych realizacjach stwierdzono ich podobieñstwo dla algorytmów D¥, MD¥ oraz MFD. Zosta³o ono potwier-dzone równie¿ w wynikach oceny zagro¿enia erozyjnego poprzez bardzo wysokie wartoœci wspó³czynnika Kappa, obliczone dla par map przedstawiaj¹cych tereny o zagro¿eniu co najmniej œrednim. Wysokie podobieñstwo uzyskanych rezultatów wynika niew¹tpliwie z po-dobnej zasady dzia³ania tych algorytmów.

W wysokim stopniu skorelowane by³y równie¿ wartoœci czynnika topograficznego uzy-skane z wykorzystaniem algorytmów DEMON i KRA. W tym przypadku wyt³umaczeniem jest zapewne fakt, i¿ w algorytmie DEMON do okreœlenia kierunku sp³ywu stosuje siê podej-œcie z algorytmu KRA. Jednak w przypadku tego pierwszego, w odró¿nieniu od KRA i pozosta³ych algorytmów, przebieg drogi sp³ywu modelowany jest w dwóch wymiarach. Szerokoœæ œcie¿ki sp³ywu pozostaje sta³a na terenach p³askich, zwiêksza siê w obszarach dywergencji i zmniejsza na obszarach konwergencji. W efekcie, w wyniku dzia³ania tego algorytmu uzyskiwano, dla takiego samego jak w przypadku KRA przebiegu œcie¿ki sp³ywu, wiêksze wartoœci jednostkowej powierzchni zasilania. Znalaz³o to swoje odzwierciedlenie w wynikach oceny zagro¿enia erozyjnego. Obszary wskazane jako zagro¿one przy zastosowa-niu algorytmu KRA prawie w ca³oœci zosta³y wskazane równie¿ w przypadku zastosowania algorytmu DEMON. Jednak w tym drugim wypadku zagro¿enie prognozowane by³o rów-nie¿ dla dodatkowych fragmentów modelowanych zlewni.

Wyniki otrzymane przy u¿yciu algorytmów D8 i Rho8 odbiega³y od pozosta³ych. Zauwa¿al-ne jest to zw³aszcza w przypadku zlewni D³ubni i widoczZauwa¿al-ne w niskich wartoœciach wspó³czyn-nika Kappa uzyskanych dla map obszarów zagro¿enia erozyjnego. Prowadzi to do wniosku, i¿ ró¿nice w dzia³aniu tych algorytmów w stosunku do pozosta³ych, ulegaj¹ wzmocnieniu na obszarach o mniejszym nachyleniu, gdzie w przypadku algorytmów umo¿liwiaj¹cych rozp³y-wanie siê wody strumieñ rzadziej kierowany jest do pojedynczej po³o¿onej ni¿ej komórki.

Jeœli wyniki przeprowadzonego modelowania rozpatrywaæ z punktu widzenia zagro¿enia erozyjnego okreœlanego dla ca³oœci zlewni, to obserwowane ró¿nice nie maj¹ istotnego zna-czenia. Znaczenie to roœnie jednak wraz ze zwiêkszaniem skali opracowania. Coraz czêœciej modelowanie erozji przeprowadza siê z rozdzielczoœci¹ przestrzenn¹ kilkunastu (np. Drze-wiecki i in., 2013) lub nawet kilku (Prasuhn i in., 2013) metrów. Wp³yw doboru algorytmu na wyniki oceny zagro¿enia erozyjnego bêdzie w takich przypadkach wiêkszy z dwóch powodów. Z jednej strony, jak wynika z badañ Erskine’a i in. (2006), roœnie on wraz ze wzrostem rozdzielczoœci. Zdaniem cytowanych autorów wynika to z dwóch przyczyn. Wiêksza rozdzielczoœæ pozwala na dok³adniejsze modelowanie sp³ywu powierzchniowego, stwarza-j¹c zarówno mo¿liwoœæ okreœlenia wiêkszej liczby œcie¿ek sp³ywu przed osi¹gniêciem jego koncentracji, jak i zwiêkszaj¹c determinowan¹ przez powierzchnie pojedynczej komórki do-k³adnoœæ oszacowania powierzchni zasilania. (Erskin i in., 2006). Jednoczeœnie obszary od-niesienia (jednostki przestrzenne), dla których (jako ca³oœci) dokonywane jest oszacowanie

zagro¿enia erozyjnego zmniejszaj¹ swoj¹ powierzchnie, nawet do pojedynczych pól (Pra-suhn i in., 2013).

Podsumowuj¹c nale¿y stwierdziæ, i¿ przeprowadzone badania potwierdzaj¹ wnioski sfor-mu³owane w publikacjach Erskine i in. (2006) oraz Mendicino (1999), którzy uznaj¹ algoryt-my umo¿liwiaj¹ce rozp³yw strumienia wody do wielu komórek za bardziej odpowiednie dla okreœlania powierzchni zasilania. Na niekorzyœæ algorytmów kieruj¹cych sp³yw do pojedyn-czej komórki przemawia dodatkowo ich wiêksza czu³oœæ na b³êdy mumerycznego modelu terenu (Desmet i Govers, 1996a). Bior¹c pod uwagê fakt, i¿ algorytmy umo¿liwiaj¹ce roz-p³yw strumienia sroz-p³ywu do wielu (wiêcej ni¿ dwóch) komórek mog¹ zawodziæ w przypadku znacznej wklês³oœci form terenowych oraz linii drena¿owych (Quinn i in., 1991; Freeman 1991), za bardzo interesuj¹cy z punktu widzenia modelowania procesów erozji wodnej gleb uznaæ nale¿y algorytm DEMON. Niestety, ze wzglêdu na swoj¹ z³o¿onoœæ, jest on rzadko implementowany w systemach GIS i stosowany w praktyce modelowania erozji.

Literatura

Biesemans J., Van Meirvenne M., Gabriels D., 2000: Extending the RUSLE with the Monte Carlo error propagation technique to predict long-term average off-site sediment accumulation. Journal of Soil and

Water Conservation 55 (1): 35-42.

Beven K.J., Brazier R.E., 2011: Dealing with Uncertainty in Erosion Model Predictions. [w:] Morgan R.P.C., Nearing M.A. (red.) Handbook of Erosion Modelling, Wiley-Blackwell: 52-79.

Cohen J., 1960: A coefficient of agreement for nominal scales. Educational and Psychological Measurement 20 (1): 37-46.

Costa-Cabral M., Burges S.J., 1994: Digital elevation model networks (DEMON): A model of flow over hillslopes for computation of contributing and dispersal areas. Water Resources Research 30: 1681-1692. Desmet P.J., Govers G., 1996a: Comparison of routing algorithms for digital elevation models and their implications for predicting ephemeral gullies. International Journal of Geographical Information Systems 10: 311-331.

Desmet P.J., Govers G., 1996b: A GIS Procedure for Automatically Calculating the USLE LS Factor on Topographically Complex Landscape Units. Journal of Soil and Water Conservation 51: 427-433. Drzewiecki W., 2010: Badanie zmian przestrzennych struktury u¿ytkowania i funkcji krajobrazu w oparciu

o wieloczasowe obrazy teledetekcyjne jako wsparcie dla planowania krajobrazu. Sprawozdanie meryto-ryczne z projektu badawczego MNiSW nr N526 029 32/2621. Katedra Geoinformacji, Fotogrametrii i Teledetekcji Œrodowiska, Wydzia³ Geodezji Górniczej i In¿ynierii Œrodowiska, AGH w Krakowie. Ma-szynopis.

Drzewiecki W., Mularz S., 2001: Modelowanie erozji wodnej gleb z wykorzystaniem GIS. [w:] Nowoczesne technologie w geodezji i in¿ynierii œrodowiska: konferencja naukowa z okazji jubileuszu 50-lecia Wydzia³u Geodezji Górniczej i In¿ynierii Œrodowiska, AGH, Kraków, 169-186.

Drzewiecki W., Mularz S., 2008: Simulation of water soil erosion effects on sediment delivery to Dobczyce Reservoir. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences XXXVII B8: 787-794.

Drzewiecki W., Mularz S., Twardy S., Kopacz M., 2008: Próba kalibracji modelu RUSLE/SDR dla oceny ³adunku zawiesiny wprowadzanego do zbiornika Dobczyckiego ze zlewni bezpoœredniej. Archiwum

Fotogrametrii, Kartografii i Teledetekcji 18: 83-98.

Drzewiecki W., Wê¿yk P., Pierzchalski M., Szafrañska B., 2013: Quantitative and Qualitative Assessment of Soil Erosion Risk in Ma³opolska (Poland), Supported by an Object-Based Analysis of High-Resolution Satellite Images. Pure and Applied Geophysics, DOI:10.1007/s00024-013-0669-7.

Erskine R.H., Green T.R., Ramirez J.A., MacDonald L.H., 2006: Comparison of grid-based algorithms for computing upslope contributing area. Water Resources Research 42, DOI:10.1029/2005WR004648.

Fairfield J., Leymarie P., 1991: Drainage networks from grid digital elevation model. Water Resources

Rese-arch 27(5): 709-717.

Freeman G.T., 1991: Calculating catchment area with divergent flow based on a regular grid. Computers and

Geosciences 17: 413-422.

Józefaciuk A., Józefaciuk C., 1996: Mechanizm i wskazówki metodyczne badania procesów erozji. Bibliote-ka Monitoringu ŒrodowisBibliote-ka, PIOŒ, Warszawa.

Kondracki J., 1988: Geografia fizyczna Polski. PWN, Warszawa.

Kulas G., 2011: Wp³yw zmian pokrycia i u¿ytkowania terenu w zlewniach Pr¹dnika i D³ubni na zagro¿enie gleb erozj¹ wodn¹. Praca magisterska. Wydzia³ Geodezji Górniczej i In¿ynierii Œrodowiska, AGH w Krakowie.

Lea N.L., 1992: An aspect driven kinematic routing algorithm. [w:] A.J. Parsons i A.D. Abrahams (red.) Overland Flow: Hydraulics and Erosion Mechanics, Chapman and Hall, New York: 147-175.

Licznar P., 2004: Prognozowanie erozyjnoœci deszczy w Polsce na podstawie miesiêcznych sum opadów.

Archiwum Ochrony Œrodowiska 30 (4): 29 -39.

Mendicino G., 1999: Sensitivity Analysis on GIS Procedures for the Estimate of Soil Erosion Risk. Natural

Hazards 20: 231-253.

O'Callaghan J.F., Mark D.M., 1984: The extraction of drainage networks from digital elevation data.

Compu-ter Vision, Graphics and Image Processing 28: 323-344.

Quinn P.F., Beven K.J., Chevallier P., Planchon O., 1991: The prediction of hillslope flow paths for distribu-ted hydrological modelling using digital terrain model. Hydrological Processes 5: 59-79.

Prasuhn V., Liniger H., Gisler S., Herweg K., Candinas A., Clément J.P., 2013: A highresolution soil erosion risk map of Switzerland as strategic policy support system. Land Use Policy 32: 281-291.

Renard K.G., Foster G.R., Weesies G.A., McCool D.K., Yoder D.C. 1997: Predicting Soil Erosion by Water: A Guide to Conservation Planning With the Revised Universal Soil Loss Equation (RUSLE). U.S. Depart-ment of Agriculture. Agriculture Handbook No. 703.

Risse L.M., Nearing M.A., Nicks A.D., Laflen J.M., 1993: Error assessment in the Universal Soil Loss Equation. Soil Science Society of America Journal 57: 825-833.

Seibert J., McGlynn B.L., 2007: A new triangular multiple flow direction algorithm for computing upslope areas from gridded digital elevation models. Water Resources Research 43, DOI:10.1029/2006WR005128. Stuczyñski T., Koza P., £opatka A., Duer I., Jadczyszyn J., 2010: Raport z analizy wskaŸników produktu, rezultatu i oddzia³ywania okreœlonych dla osi 2 PROW 2007-2013 oraz wybranych pytañ oceniaj¹cych zawartych w podrêczniku wspólnych ram monitorowania i oceny. Wytyczne (CMEF) wraz z okreœle-niem Ÿróde³ i dostêpnoœci danych. IUNG-PIB, Pu³awy.

www.minrol.gov.pl/pol/content/download/28450/158380/file/Rap_z_analizy_wskaznikow.pdf

Tarboton D.G., 1997: A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research 33 (2): 309-319.

Tetzlaff B., Friedrich K., Vorderbrügge T., Vereecken H., Wendland F., 2013: Distributed modelling of mean annual soil erosion and sediment delivery rates to surface waters. Catena 102: 13-20.

Wawer R., Nowocieñ E., Podolski B., 2005: Real and Calculated K USLE Erodibility Factor for Selected Polish Soils. Polish Journal of Environmental Studies 14 (5): 655-658.

Wischmeier W.H., Smith D.D., 1978: Predicting Rainfall Erosion Losses – A Guide to Conservation Planning. USDA Handbook 537, Washington, D.C.

Abstract

The objective of the presented research was to evaluate the influence of flow routing algorithm on RUSLE-based soil erosion assessment on catchment scale. Both the amount of soil loss and the pattern of erosion risk classes were compared. Seven flow routing algorithms were tested: D8, Rho8, Kinema-tic Routing, DEMON, D¥, Multiply Flow Direction (MFD) and Triangular Multiply Flow Direction (MD¥).

Based on the results achieved, we can conclude that the choice of flow routing algorithm influences the soil erosion risk assessment. Particular approaches gave similar results when total area of endange-red soils in the catchment was consideendange-red. However, the patterns of erosion differ. Multiple-direction algorithms (especially DEMON) seem to be better suited for water erosion studies.

dr in¿. Wojciech Drzewiecki drzewiec@agh.edu.pl tel. +48 12 617 2288 mgr in¿. Sebastian Ziêtara sebastian.zietara@gmail.com