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.

7.1.2 Model semiwariogramu

Model semiwariogramu składa się zazwyczaj z trzech podstawowych elementów. 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

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().

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

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

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 kryteró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 modelu. 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.

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. Średnio te płaty mają średnicę określoną przez zasięg (range) modelu.

##   model psill range
## 1   Sph    10  3000

##   model   psill    range
## 1   Sph 12.5445 4440.768

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.

##   model psill range
## 1   Exp    10  3000

##   model    psill    range
## 1   Exp 16.50871 3084.309

7.4.3 Model gaussowski

Model gaussowski (Gau) również posiada zasięg praktyczny definiowany jako 95% wartości wariancji progowej. 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.

##   model psill range
## 1   Gau    13  3000

##   model    psill    range
## 1   Gau 8.085963 798.7454

7.4.4 Model potęgowy

Model potęgowy (Pow) to przykład tzw. modelu nieograniczonego. 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.

##   model psill range
## 1   Pow  0.03   0.7

##   model      psill     range
## 1   Pow 0.02732273 0.7382382

7.4.5 Porównanie modeli

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). 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

##   model      psill    range
## 1   Nug  0.7346354    0.000
## 2   Sph 13.2621876 5602.027

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.

##   model psill range
## 1   Nug   0.5     0
## 2   Gau   4.0   500
## 3   Gau  10.0  3000

##   model    psill    range
## 1   Nug 0.500000    0.000
## 2   Gau 6.131863 1532.805
## 3   Gau 6.131863 1532.805

7.5 Modelowanie anizotropowe

7.5.1 Anizotropia

Uwzględnienie anizotropii wymaga zamiany parametru zasięgu na trzy inne parametry:

  • 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

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).

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?