Alfred Stach Instytut Paleogeografii i Geoekologii

GEOSTATYSTYKAWykłady dla III roku Geografiispecjalność – geoinformacjaKriging wartości kodowanych(Indicator Kriging) Alfred Stach Instytut Paleogeografii i Geoekologii Wydział Nauk Geograficznych i Geologicznych UAM

Lokalizacja stanowisk pomiarowych opadów atmosferycznych na profilu

Maksymalne sumy dobowe opadów zarejestrowane w maju 1980 roku na posterunkach usytuowanych na profilu

Skumulowany rozkład prawdopodobieństwa maksymalnych dobowych sum opadów zarejestrowanych w posterunkach opadowych na terenie Polski w maju 1980 rokuNa wykresie zaznaczono wysokości sum dobowych o prawdopodobieństwie przewyższenia 0,2, 0,4, 0,6, 0,8, 0,9 i 0,95 (percentyle 20, 40, …., 95%)

Wartości progowe maksymalnych sum dobowych opadów w maju 1980 roku o prawdopodobieństwie 0,2, 0,4, 0,6, 0,8, 0,9 i 0,95 (percentyle 20, 40, 60, 80, 90 i 95%) naniesione na dane profilowe

Maksymalne sumy dobowe opadów zarejestrowane w maju 1980 roku na analizowanym profilu przekodowane na wektory danych binarnych w zależności od przekroczenia wartości progowych wyznaczonych z globalnej krzywej skumulowanego rozkładu prawdopodobieństwa (precentyle 20, 40, 60, 80, 90 i 95%)

Semiwariogramy empirycznei ich modele dla wartości kodowanych (percentyle 20, 40, 60, 80, 90 i 95%) maksymalnych sum dobowych opadów na terenie Polski w maju 1980 roku

Estymowany metodą IK profil prawdopodobieństwa maksymalnych opadów dobowych w maju 1980. Pionowymi liniami przerywanymi zaznaczono lokalizację punktów pomiarowych

Wartość oczekiwana (E-mean) maksymalnego opadu dobowego na analizowanym profilu w maju 1980 wyliczona z ccdf estymowanych metodą IK. Zacieniowany pas oznacza zakres odchylenia standardowego estymacji (Conditional Variance). Zaznaczono lokalizację punktów pomiarowych i wysokości rzeczywiście zmierzonych maksymalnych opadów dobowych

Maksymalny opad dobowy na analizowanym profilu w maju 1980 o prawdopodobieństwie wystąpienia 0,9 (A) oraz prawdopodobieństwo wystąpienia opadu dobowego większego lub równego 25 mm (B). Zaznaczono lokalizację punktów pomiarowych i wysokości rzeczywiście zmierzonych maksymalnych opadów dobowych

Estymowane metodą IK warunkowe skumulowane rozkłady prawdopodobieństwa (ccdf) maksymalnych opadów dobowych w maju 1980 roku w trzech lokalizacjach (u1, u2 i u3) na analizowanym profilu. Zaznaczono globalne cdf (V-80) obliczone dla wszystkich danych pomiarowych z całej Polski, a także wartości sum opadów dla cdf odpowiadające prawdopodobieństwu 0,9 Błędy relacji porządkowych ccdf dla lokalizacji u2 i u3 (odpowiednio 0,000418 i 0,0068571)

Błędy relacji porządkowych • Podstawową wadą krigingu wartości kodowanych (IK) jest występowanie błędów relacji porządkowych. W dowolnej lokalizacji u, każde estymowane posteriori prawdopodobieństwo [F(u;zk(n))]* musi należeć do przedziału [0,1], a seria K takich szacunków musi być niemalejącą funkcją wielkości wartości progowej zk:

Błędy relacji porządkowych I • Występowanie błędów relacji porządkowych pierwszego rodzaju wynika z samej natury algorytmu krigingu, który jest liniową, nie wypukłą, kombinacją danych pomiarowych. • Pociąga to za sobą możliwość obliczenia ujemnych wag dla poszczególnych danych pomiarowych znajdujących się w zasięgu sąsiedztwa szukania. • Sytuacja taka ma miejsce jeśli zachodzi zjawisko ekranowania, tj. zlokalizowany bliżej punktu estymacji u0 punkt danych u2 częściowo „niweluje” wpływ leżącego dalej na tym samym kierunku punktu u1. • Ta cecha algorytmu ma zarówno zalety, jak i wady. Z jednej strony umożliwia uzyskanie estymacji, które wykraczają poza zakres danych pomiarowych, z drugiej mogą być to czasami wyniki nierealistyczne, takie jak ujemne stężenia, czy proporcje większe od 1. • Błędy tego rodzaju występują częściej, i ich rozmiary są większe, w zwykłym krigingu (OK) niż w prostym krigingu (SK), oraz w wielozmiennym kokrigingu niż w krigingu. Jest to efektem występujących w owych algorytmach (OK, SCK, OCK) ograniczeń wielkości wag (wymuszających ich sumowanie do 1 lub do 0)

Ekranowanie danych w krigingu Ilustracja sytuacji występowania ujemnych wag w algorytmie zwykłego krigingu (OK). A – wykres i wzór sferycznego modelu semiwariogramu użytego w obliczeniach: wariancja nuggetowa (C0) = 10, wariancja progowa (C1) = 90, zasięg (a) = 100 jednostek. B – układ przestrzenny estymowanej lokalizacji (0) i punktów danych (1-5) oraz wartości wyliczonych dla tej konfiguracji wag OK. Sytuacja bez ekranowania. C – układ przestrzenny estymowanej lokalizacji (0) i punktów danych (1-5) oraz wartości wyliczonych dla tej konfiguracji wag OK. Sytuacja z ekranowaniem punktu 1 przez punkt 2. Cieniowany okrąg (B i C) wskazuje na zasięg autokorelacji.

Błędy relacji porządkowych II • Występowanie błędów drugiego rodzaju wynika dodatkowo z faktu, że każde z K prawdopodobieństw jest estymowane osobno, oraz że często w konkretnych klasach z (przedziałach wartości analizowanej cechy) w lokalnym sąsiedztwie brak jest danych pomiarowych. Jakie ma to konsekwencje zaprezentowano w poniższym przykładzie. • Zakładamy że w klasie (z7, z8] nie ma danych pomiarowych. Obie estymacje IK dla wartości progowych z7 i z8 oparte są zatem na tym samym zbiorze danych kodowanych ponieważ: • Różnice między tymi dwoma estymacjami IK są wówczas jedynie efektem liniowej kombinacji różnic pomiędzy wagami IK dla obu wartości progowych z7 i z8:

Błędy relacji porządkowych II • Wartość ujemna różnicy pociąga za sobą naruszenie relacji porządkowej. W sytuacji kiedy oba modele semiwariogramów I(h, z7) i I(h, z8) są identyczne, także oba zbiory wag IK będą takie same, ponieważ dla obu wartości progowych w obliczeniach zostaną wykorzystane te same lokalizacje danych pomiarowych: • Różnica wynosi wówczas zero, stąd nie ma naruszenia relacji porządkowej. W przeciwnym wypadku, istotnej różnicy dwóch kolejnych modeli semiwariogramów wartości kodowanych, w tym przypadku między progami z7 i z8, powstają dwa odmienne zbiory wag IK pociągając za sobą ryzyko wystąpienia błędów relacji porządkowych.

Sposoby eliminacji błędów relacji porządkowych • Błędy relacji porządkowych są w estymacjach IK stosunkowo częste, ale ich rozmiar jest zazwyczaj niewielki – około 0,01. Aby ograniczyć ich ilość i rozmiar stosuje się dwie strategie: • Błędy drugiego rodzaju w zasadzie są łatwe do wyeliminowania jeśli dla wszystkich wartości progowych użyje się tego samego modelu struktury przestrzennej – semiwariogramu. • To dość radykalne podejście jest często stosowane pod nazwą median Indicator Kriging (mIK). Nazwa sugeruje, i rzeczywiście tak bywa najczęściej, że w algorytmie tym stosuje się model struktury przestrzennej danych kodowanych w stosunku do wartości mediany (50 percentyla). Nie jest jednakże jakaś ścisła reguła. • Zalety mIK związane są nie tylko z eliminacją większości naruszeń relacji porządkowych. Jest to przede wszystkim metoda mniej pracochłonna – modelowanie jednego semiwariogramu zamiast kilku, czy kilkunastu, ale przede wszystkim znacznie szybsza w obliczeniach. Dla każdej lokalizacji (węzła siatki interpolacyjnej) obliczany jest bowiem tylko jeden układ równań krigingu. • Popularność mIK wynika również z faktu, że mimo tak znacznego uproszczenia procedury, uzyskiwane wyniki są zazwyczaj tylko nieznacznie gorsze od uzyskanych za pomocą „pełnego” krigingu wartości kodowanych.

Sposoby eliminacji błędów relacji porządkowych • W sytuacji kiedy nie można zastosować metody mIK zaleca się takie modelowanie struktury przestrzennej dla kolejnych wartości progowych, aby unikać gwałtownych zmian parametrów modeli. Można to osiągnąć na przykład poprzez użycie dla wszystkich wartości progowych różnych kombinacji liniowych tych samych elementarnych struktur. • Parametry modeli semiwariogramów danych kodowanych (wariancja progowa, zasięg, kierunek i proporcja anizotropii) powinny zmieniać się stopniowo od jednej wartości progowej do następnej. Nie jest to zazwyczaj żadne istotne ograniczenie, ponieważ w „naturze” zmiany struktury przestrzennej dla różnych klas wielkości analizowanego parametru zazwyczaj zachodzą w sposób stopniowy – płynny. • Zupełnie inne podejście do problemu redukcji błędów relacji porządkowych zakłada nie „sztywne” ustalenie jednej serii wartości progowych zk, ale ich dynamiczną modyfikację osobno dla każdego sąsiedztwa szukania w zależności od zakresu wartości tam występujących. Unika się w ten sposób, często w tradycyjnym IK występującej sytuacji, że w pewnych klasach wielkości nie ma danych pomiarowych. Potrzebne odpowiednie modele semiwariogramów dla zmiennych wartości progowych są interpolowane z podanych wcześniej przez „operatora”.

Usuwanie błędów relacji porządkowych • Wymienione procedury redukują, ale całkowicie nie eliminują problemu naruszeń relacji porządkowych. Dlatego też konieczna jest dodatkowa, finalna operacja korekty uzyskanych za pomocą algorytmu IK wartości ccdf. • Najczęściej stosuje się prostą procedurę uśredniana korekt wartości rosnących i malejących: Błędy relacji porządkowych uzyskanych z obliczeń IK wartości ccdf i ich korekta. Objaśnienia: a – „niezależne” wartości ccdf wyliczone algorytmem IK, b – korekta wartości rosnących (upward correction), c – korekta wartości malejących (downward correction), d – wynikowe ccdf uzyskane z uśrednienia obu wartości skorygowanych.

Przykład korekty relacji porządkowych Przykłady raportów dotyczące ilości i rozmiarów korekt relacji porządkowych warunkowych kumulacyjnych funkcji rozkładu maksymalnych sum dobowych opadów: A – jednowymiarowy przykład (profil) z maja 1980, B – maksymalne sumy dobowe opadów w roku 1974 na całym obszarze Polski.

Przykład korekty relacji porządkowych

Interpolacja i ekstrapolacja wynikowej ccdf • Działanie algorytmu IK można porównać do korekty, czy też modyfikacji, na podstawie informacji lokalnych, globalnego dyskretnego cdf. • Otrzymujemy w efekcie punktową, dyskretną, warunkową funkcję rozkładu prawdopodobieństwa (ccdf). • Aby móc ją w pełni wykorzystać do różnorodnych zastosowań, musimy w ostatnim etapie obliczeń dokonać operacji odwrotnej do tej która rozpoczynała całą procedurę – z dyskretnej, nieciągłej ccdf uzyskać z powrotem rozkład ciągły. • Praktycznie rzecz biorąc pociąga to za sobą konieczność ustalenia sposobu za pomocą którego można oszacować dowolną wartość ccdf, a nie tylko dla K wybranych progów.

Interpolacja i ekstrapolacja wynikowej ccdf Problem ten zazwyczaj „rozbija się” na dwa cząstkowe: (1) interpolację ccdf w obrębie klas wyznaczonych przez kolejne wartości progowe, (2) ekstrapolację poza progami skrajnymi, tj. minimalnym i maksymalnym Budowa ciągłego ccdf dla lokalizacji u1 z jednowymiarowego przykładu obliczeń krigingu wartości kodowanych. Objaśnienia: A – ciągły cdf dla całego zbioru danych (a) i dyskretny ccdf uzyskany z obliczeń IK dla lokalizacji u1, B – to samo co w A plus: d – ekstrapolacja potęgowa dolnego ogona rozkładu ( = 4,0), c – interpolacja liniowa pomiędzy granicami klas i, b – ekstrapolacja hiperboliczna górnego ogona rozkładu ( = 2,5).

Interpolacja i ekstrapolacja wynikowej ccdf • Do interpolacji ccdf pomiędzy wartościami progowymi (zk-1, zk) wykorzystywany jest zazwyczaj model liniowy. Używając tego modelu zakładamy istnienie w klasach rozkładu równomiernego • Do ekstrapolacji dolnego ogona rozkładu używany jest najczęściej model potęgowy • Do ekstrapolacji górnego ogona używany jest model potęgowy lub hiperboliczny Potęgowa (a) i hiperboliczna (b) interpolacja / ekstrapolacja skumulowanego rozkładu zmiennej

Zestawienie wad i zalet IK WADY: • utrata części informacji ze względu dyskredytyzację ciągłej dystrybuanty empirycznej, • pracochłonność – konieczność czasochłonnego budowania modelu semiwariancji dla każdej wartości progowej; często występujące trudności w określeniu modeli dla wartości bardzo niskich i bardzo wysokich zmuszają do subiektywnych decyzji, a te rodzą wątpliwości co do optymalności uzyskanych estymacji, • wykraczanie estymowanych prawdopodobieństw poza dopuszczalny zakres (0, 1), oraz błędy w ich relacjach porządkowych, • arbitralnie przyjmowana metoda interpolacji/ekstrapolacji uzyskanej warunkowej dystrybuanty. ZALETY: • potwierdzona w dziesiątkach zastosowań i testów metodycznych skuteczność, • brak trudnych do weryfikacji założeń dotyczących rozkładu statystycznego populacji (metoda nieparametryczna), • żadna z alternatywnych metod nie jest wyraźnie lepsza, • alternatywne metody są bardziej skomplikowane = bardziej „podatne” na błędy metodyczne, • łatwa możliwość uwzględnienia danych uzupełniających („twardych” i „miękkich”). • powszechna dostępność oprogramowania

Analizowane dane Fragment doliny lodowca Ebba na Spitsbergenie Zachodnim ok. 7843’N i 1644’E

Analizowane dane Dolina Ebby – analizowany obszar:

Analizowane dane Zdjęcie satelitarne Aster – Terra z 13 lipca 2002 roku Światło widzialne i bliska podczerwień. Rozdzielczość – 15 m

Analizowane dane Fragment mapy geomorfologicznej otoczenia fiordu Petuniabukta (Karczewski 1990).

Zdjęcie z 13.VII. 2002 r. Fragment mapy geomorfologicznej Zdjęcie z 5.VIII. 2002 r. Oryginał i klasyfikacja obszar 1500  1950 m = 2,925 km2(100  130 pikseli = 13000 danych)kanał 3n – bliska podczerwień250 podstawowych losowych próbek100 dodatkowych losowych próbek Analizowane dane:

1 2 3 4 Obraz rzeczywisty Interpolacja danych jakościowych– semiwariogramy kategorii

1 2 Statystyki klasyfikacji przestrzennej 3 4 Obraz rzeczywisty Interpolacja IK Prawdo- podobieństwo przynależności do regionu (klasy) Interpolacja danych jakościowych– kriging kategorii (IK)

SCHEMAT OPRÓBOWANIA STOKU.POBÓR RDZENI GLEBOWYCH.

POBÓR RDZENIGLEBOWYCHI OPISBARWY GLEB

Empiryczne semiwariogramy wskaźnikowe (indicator semivariogram)dla poszczególnych klas barw poziomu akumulacyjno-próchnicznego na stoku A

Prawdopodobieństwo przynależności do klas barwpoziomu akumulacyjno-próchnicznego na stoku A

Alfred Stach Instytut Paleogeografii i Geoekologii