7 Modelowanie autokorelacji przestrzennej

Odtworzenie obliczeń z tego rozdziału wymaga załączenia poniższych pakietów oraz wczytania poniższych danych:

7.1 Modelowanie matematycznie autokorelacji przestrzennej

7.1.1 Modelowanie struktury przestrzennej

Semiwariogram empiryczny (wyliczony z danych punktowych) jest:

  • Nieciągły - wartości semiwariancji są średnimi przedziałowymi.
  • Chaotyczny - badana próba jest jedynie przybliżeniem rzeczywistości, dodatkowo obciążonym błędami.

Estymacje i symulacje przestrzenne wymagają modelu struktury przestrzennej analizowanej cechy, a nie tylko wartości empirycznych. Dodatkowo, matematycznie modelowanie wygładza chaotyczne fluktuacje danych empirycznych (rycina 7.1).

Porównanie semiwariogramu empirycznego i modelu semiwariogramu.

Rycina 7.1: Porównanie semiwariogramu empirycznego i modelu semiwariogramu.

7.1.2 Model semiwariogramu

Model semiwariogramu składa się zazwyczaj z trzech podstawowych elementów (rycina 7.2). Są to:

  • Nugget - efekt nuggetowy - pozwala na określenie błędu w danych wejściowych oraz zmienności na dystansie krótszym niż pierwszy odstęp.
  • Sill - semiwariancja progowa - oznacza wariancję badanej zmiennej.
  • Range - zasięg - to odległość do której istnieje przestrzenna korelacja.
Podstawowe elementy modelu semiwariogramu.

Rycina 7.2: Podstawowe elementy modelu semiwariogramu.

7.1.3 Model nuggetowy

Model nuggetowy określa sytuację, w której analizowana zmienna nie wykazuje autokorelacji. Inaczej mówiąc, niepodobieństwo jej wartości nie wzrasta wraz z odległością. Model nuggetowy nie powinien być używany samodzielnie - w większości zastosowań jest on elementem modelu złożonego. Służy on do określania, między innymi, błędu pomiarowego czy zmienności na krótkich odstępach.

7.2 Modele podstawowe

7.2.1 Typy modeli podstawowych

Pakiet gstat zawiera 20 podstawowych modeli geostatystycznych, w tym najczęściej używane takie jak:

  • Nuggetowy (ang. Nugget effect model)
  • Sferyczny (ang. Spherical model)
  • Gaussowski (ang. Gaussian model)
  • Potęgowy (ang. Power model)
  • Wykładniczy (ang. Exponential model)
  • Inne

Do wyświetlenia listy nazw modeli i ich skrótów służy funkcja vgm().

##    short                                      long
## 1    Nug                              Nug (nugget)
## 2    Exp                         Exp (exponential)
## 3    Sph                           Sph (spherical)
## 4    Gau                            Gau (gaussian)
## 5    Exc        Exclass (Exponential class/stable)
## 6    Mat                              Mat (Matern)
## 7    Ste Mat (Matern, M. Stein's parameterization)
## 8    Cir                            Cir (circular)
## 9    Lin                              Lin (linear)
## 10   Bes                              Bes (bessel)
## 11   Pen                      Pen (pentaspherical)
## 12   Per                            Per (periodic)
## 13   Wav                                Wav (wave)
## 14   Hol                                Hol (hole)
## 15   Log                         Log (logarithmic)
## 16   Pow                               Pow (power)
## 17   Spl                              Spl (spline)
## 18   Leg                            Leg (Legendre)
## 19   Err                   Err (Measurement error)
## 20   Int                           Int (Intercept)

Można się również im przyjrzeć używając funkcji show.vgms() (rycina 7.3).

Przykład modeli semiwariogramu dostępnych w pakiecie gstat.

Rycina 7.3: Przykład modeli semiwariogramu dostępnych w pakiecie gstat.

Istnieje możliwość wyświetlenia tylko wybranych modeli podstawowych poprzez argument models (rycina 7.4).

Porównanie modeli sferycznego (Sph), gaussowskiego (Gau), potęgowego (Pow) i wykładniczego (Exp).

Rycina 7.4: Porównanie modeli sferycznego (Sph), gaussowskiego (Gau), potęgowego (Pow) i wykładniczego (Exp).

Dodatkowo, można je porównać na jednym wykresie poprzez argument as.groups = TRUE (rycina 7.5).

Porównanie modeli sferycznego (Sph), gaussowskiego (Gau), potęgowego (Pow) i wykładniczego (Exp) na jednym wykresie.

Rycina 7.5: Porównanie modeli sferycznego (Sph), gaussowskiego (Gau), potęgowego (Pow) i wykładniczego (Exp) na jednym wykresie.

7.3 Metody modelowania

7.3.1 Rodzaje metod modelowania

Istnieją trzy najczęściej spotykane metody modelowania geostatystycznego:

  • Ustawianie “ręczne” parametrów modelu, np. funkcja vgm() z pakietu gstat.
  • Ustawianie “wizualne” parametrów modelu, np. funkcja eyefit() z pakietu geoR.
  • Automatyczny wybór parametrów na podstawie różnych kryteriów statystycznych, np. funkcja fit.variogram() z pakietu gstat(), variofit() z pakietu geoR, autofitVariogram() z pakietu automap.

Odpowiednie określenie modelu matematycznego często nie jest proste. W efekcie automatyczne metody nie zawsze są w stanie dać lepszy wynik od modelowania “ręcznego”. Najlepiej, gdy wybór modelu oparty jest o wiedzę na temat zakładanego procesu przestrzennego.

7.3.2 Funkcja fit.variogram()

Funkcja fit.variogram() z pakietu gstat dopasowuje zasięg oraz semiwariancję progową w oparciu o ustalone “ręcznie” wejściowe parametry modelu3.

7.4 Modelowanie izotropowe

Do zbudowania modelu semiwariogramu należy wykonać szereg kroków:

  1. Stworzyć i wyświetlić semiwariogram empiryczny analizowanej zmiennej z użyciem funkcji variogram() oraz plot().
  2. Zdefiniować wejściowe parametry semiwariogramu. W najprostszej sytuacji wystarczy zdefiniować używany model/e poprzez skróconą nazwę używanej funkcji (model). Możliwe, ale nie wymagane jest także określenie wejściowej semiwariancji cząstkowej (psill) oraz zasięgu modelu (range) w funkcji vgm(). Uzyskany model można przedstawić w funkcji plot() podając nazwę obiektu zawierającego semiwariogram empiryczny oraz obiektu zawierającego model.
  3. Dopasować parametry modelu używając funkcji fit.variogram(). To dopasowanie można również zwizualizować używając funkcji plot().

7.4.1 Model sferyczny

Model sferyczny (Sph) jest jednym z najczęściej stosowanych modeli geostatystycznych. Reprezentuje on cechę, której zmienność wartości ma charakter naprzemiennych płatów niskich i wysokich wartości (rycina 7.6). Średnio te płaty mają średnicę określoną przez zasięg (range) modelu.

Przykład zjawiska reprezentowanego poprzez model sferyczny.

Rycina 7.6: Przykład zjawiska reprezentowanego poprzez model sferyczny.

Stworzenie takiego modelu dla zmiennej temp polega najpierw na zbudowaniu semiwariogramu empirycznego używając funkcji variogram() (rycina 7.7).

Semiwariogram zmiennej temp.

Rycina 7.7: Semiwariogram zmiennej temp.

W kolejnym kroku, przy użyciu funkcji vgm() określa się typ modelu oraz jego podstawowe parametry (rycina 7.8).

##   model psill range
## 1   Sph    10  3000
Model sferyczny zmiennej temp.

Rycina 7.8: Model sferyczny zmiennej temp.

Dodatkowo możliwe jest także automatyczne dopasowanie modelu w oparciu o wstępnie podane parametry używając funkcji fit.variogram() (rycina 7.9).

##   model    psill    range
## 1   Sph 13.20295 4549.691
Automatycznie dopasowany model sferyczny zmiennej temp używając wstępnie podanych parametrów.

Rycina 7.9: Automatycznie dopasowany model sferyczny zmiennej temp używając wstępnie podanych parametrów.

W niektórych przypadkach funkcja fit.variogram() da także zadawalające wyniki, gdy tylko zostanie podany typ modelu (rycina 7.10).

##   model   psill    range
## 1   Sph 13.2018 4549.092
Automatycznie dopasowany model sferyczny zmiennej temp używając jedynie wstępnie podanego typu modelu.

Rycina 7.10: Automatycznie dopasowany model sferyczny zmiennej temp używając jedynie wstępnie podanego typu modelu.

7.4.2 Model wykładniczy

Model wykładniczy (Exp) również jest jednym z najczęściej używanych w geostatystyce. Od modelu sferycznego różni go szczególnie to, że nie ma on skończonego zasięgu. W jego przypadku, zamiast zasięgu podaje się tzw. zasięg praktyczny. Oznacza on odległość na jakiej model osiąga 95% wartości wariancji progowej (rycina 7.11).

Przykład zjawiska reprezentowanego poprzez model wykładniczy.

Rycina 7.11: Przykład zjawiska reprezentowanego poprzez model wykładniczy.

Proces tworzenia tego modelu jest bardzo podobny do przedstawionego powyżej - budowany jest semiwariogram empiryczny, a następnie ustalany jest model poprzez podanie kilku jego parametrów (rycina 7.12).

##   model psill range
## 1   Exp    10  3000
Model wykładniczy zmiennej temp.

Rycina 7.12: Model wykładniczy zmiennej temp.

Dalej, funkcja fit.variogram() może pomóc w dopasowaniu tego modelu (rycina 7.13).

##   model    psill    range
## 1   Exp 17.87051 3298.917
Automatycznie dopasowany model wykładniczy zmiennej temp.

Rycina 7.13: Automatycznie dopasowany model wykładniczy zmiennej temp.

7.4.3 Model gaussowski

Model gaussowski (Gau) również posiada zasięg praktyczny definiowany jako 95% wartości wariancji progowej (rycina 7.14). Jego cechą charakterystyczną jest paraboliczny kształt na początkowym odcinku. Jest on najczęściej używany do modelowania cech o regularnej i łagodnej zmienności przestrzennej. Model gaussowski z uwagi na swoje cechy zazwyczaj nie powinien być stosowany samodzielnie, lecz jako element modelu złożonego.

Przykład zjawiska reprezentowanego poprzez model gaussowski.

Rycina 7.14: Przykład zjawiska reprezentowanego poprzez model gaussowski.

Zdefiniowanie modelu gaussowskiego odbywa się używając skrótu "Gau" (rycina 7.15).

##   model psill range
## 1   Gau    13  3000
Model gaussowski zmiennej temp.

Rycina 7.15: Model gaussowski zmiennej temp.

Dopasowanie tego modelu może również nastąpić używając funkcji fit.variogram() (rycina 7.16).

##   model    psill    range
## 1   Gau 8.573835 852.2404
Automatycznie dopasowany model gaussowski zmiennej temp.

Rycina 7.16: Automatycznie dopasowany model gaussowski zmiennej temp.

7.4.4 Model potęgowy

Model potęgowy (Pow) to przykład tzw. modelu nieograniczonego (rycina 7.17). Jego wartość rośnie w nieskończoność, dlatego niemożliwe jest określenie jego zasięgu. W przypadku modelu potęgowego, parametr range oznacza wykładnik potęgowy.

Przykład zjawiska reprezentowanego poprzez model potęgowy.

Rycina 7.17: Przykład zjawiska reprezentowanego poprzez model potęgowy.

Model potęgowy jest określany skrótem "Pow" (ryciny 7.18 i 7.19).

##   model psill range
## 1   Pow  0.03   0.7
Model potęgowy zmiennej temp.

Rycina 7.18: Model potęgowy zmiennej temp.

##   model      psill     range
## 1   Pow 0.02515946 0.7535889
Automatycznie dopasowany model potęgowy zmiennej temp.

Rycina 7.19: Automatycznie dopasowany model potęgowy zmiennej temp.

7.4.5 Porównanie modeli

Z uwagi na swoją charakterystykę, każdy z powyższych modeli ma inny zakres wartości. Aby porównać te modele należy je przedstawić używając tej samej skali kolorystycznej (7.20).

Porównanie powierzchni reprezentowanej przez modele sferyczny, wykładniczy, gaussowski i potęgowy.

Rycina 7.20: Porównanie powierzchni reprezentowanej przez modele sferyczny, wykładniczy, gaussowski i potęgowy.

7.4.6 Modele złożone I

Najczęściej pojedynczy model nie jest w stanie odwzorować dokładnie zmienności przestrzennej analizowanej cechy. W takich sytuacjach konieczne jest połączenie dwóch lub więcej modeli podstawowych. Najbardziej powszechny model złożony składa się z funkcji nuggetowej (dla odległości zero) oraz drugiej funkcji (dla dalszej odległości) (rycina 7.21). Zdefiniowanie takiej funkcji odbywa się poprzez dodanie argumentu nugget w funkcji vgm().

##   model psill range
## 1   Nug   0.5     0
## 2   Sph  10.0  3000
Złożony model zmiennej temp.

Rycina 7.21: Złożony model zmiennej temp.

Dalsze dopasowanie modeli złożonych również można uzyskać używając funkcji fit.variogram() (rycina 7.22).

##   model      psill    range
## 1   Nug  0.6751142    0.000
## 2   Sph 13.7617233 5511.173
Automatycznie dopasowany złożony model zmiennej temp.

Rycina 7.22: Automatycznie dopasowany złożony model zmiennej temp.

7.4.7 Modele złożone II

Bardziej złożone modele można tworzyć z pomocą argumentu add.to. Przyjmuje on kolejny obiekt funkcji vgm() i poprzez połączenie tych dwóch obiektów otrzymuje model złożony. Na poniższym przykładzie stworzony został model złożony składający się z modelu nuggetowego oraz dwóch modeli gaussowskich (ryciny 7.23 i 7.24).

##   model psill range
## 1   Nug   0.5     0
## 2   Gau   4.0   500
## 3   Gau  10.0  3000
Złożony model zmiennej temp składający się z modelu nuggetowego oraz dwóch modeli gaussowskich.

Rycina 7.23: Złożony model zmiennej temp składający się z modelu nuggetowego oraz dwóch modeli gaussowskich.

Automatycznie dopasowany złożony model zmiennej temp składający się z modelu nuggetowego oraz dwóch modeli gaussowskich.

Rycina 7.24: Automatycznie dopasowany złożony model zmiennej temp składający się z modelu nuggetowego oraz dwóch modeli gaussowskich.

7.5 Modelowanie anizotropowe

7.5.1 Anizotropia

Uwzględnienie anizotropii wymaga zamiany parametru zasięgu na trzy inne parametry (rycina 7.25):

  • Kąt określający dominujący kierunek.
  • Zasięg w dominującym kierunku.
  • Proporcję anizotropii, czyli relację pomiędzy zasięgiem w przeciwległym kierunku do zasięgu w dominującym kierunku.
Podstawowe parametry mapy semiwariogramu.

Rycina 7.25: Podstawowe parametry mapy semiwariogramu.

W pakiecie gstat odbywa się to poprzez dodanie argumentu alpha do funkcji variogram(). Należy w niej zdefiniować analizowane kierunki, które zostały określone na podstawie mapy semiwariogramu. Następnie w funkcji vgm() należy podać nowy argument anis. Przyjmuje on dwie wartości. Pierwsza z nich (45 w przykładzie poniżej) oznacza dominujący kierunek anizotropii, druga zaś (0.4) mówi o tzw. proporcji anizotropii. Proporcja anizotropii jest to relacja pomiędzy zmiennością na kierunku prostopadłym a głównym kierunku. Na poniższym przykładzie zasięg ustalony dla głównego kierunku wynosi 4000 metrów. Wartość proporcji anizotropii, 0.4, w tym wypadku oznacza że dla prostopadłego kierunku zasięg będzie wynosił 1600 metrów (4000 metrów x 0.4) (rycina 7.26).

Mapa semiwariogramu zmiennej temp.

Rycina 7.26: Mapa semiwariogramu zmiennej temp.

Anizotropia może być także reprezentowana używając semiwariogramów kierunkowych (rycina 7.27)

Semiwariogramy kierunkowe zmiennej temp.

Rycina 7.27: Semiwariogramy kierunkowe zmiennej temp.

Następnie semiwariogramy kierunkowe mogą być modelowane ręcznie (rycina 7.28) lub też używając automatycznego dopasowania (rycina 7.29).

Modele kierunkowe zmiennej temp.

Rycina 7.28: Modele kierunkowe zmiennej temp.

Automatycznie dopasowane modele kierunkowe zmiennej temp.

Rycina 7.29: Automatycznie dopasowane modele kierunkowe zmiennej temp.

7.6 Zadania

Przyjrzyj się danym z obiektu punkty_pref. Możesz go wczytać używając poniższego kodu:

  1. Zbuduj modele semiwariogramu zmiennej srtm używając modelu sferycznego używając zarówno ręcznie ustalonych parametrów oraz funkcji fit.variogram(). Porównaj graficznie uzyskane modele.
  2. Zbuduj modele semiwariogramu zmiennej srtm używając modelu nuggetowego, sferycznego, wykładniczego, gausowskiego i potęgowego. Porównaj graficznie uzyskane modele.
  3. Stwórz złożony model semiwariogramu zmiennej srtm używając modelu nuggetowego i sferycznego.
  4. W oparciu o mapę semiwariogramu, zbuduj semiwariogramy kierunkowe zmiennej srtm dla kierunków wykazujących anizotropię przestrzenną. Następnie zbuduj modele semiwariogramu dla uzyskanych semiwariogramów kierunkowych.
  5. (Dodatkowe) Spróbuj użyć jednego z modeli podstawowych, który nie był opisywany w tym rozdziale. Czym ten wybrany model się charakteryzuje?

  1. Więcej na temat działania tej funkcji można przeczytać we wpisie na stronie https://www.r-spatial.org/r/2016/02/14/gstat-variogram-fitting.html.