7 Modelowanie matematycznie autokorelacji przestrzennej

library('gstat')
library('geostatbook')
data(punkty)

7.1 Modelowanie matematycznie autokorelacji przestrzennej

7.1.1 Modelowanie matematycznie autokorelacji 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.2 Modele podstawowe

7.2.1 Modele podstawowe

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

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

show.vgms()

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

show.vgms(models=c('Nug', 'Sph', 'Gau', 'Pow', 'Exp'), range=1.4, max=2.5)

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

show.vgms(models=c('Nug', 'Sph', 'Gau', 'Pow', 'Exp'), range=1.4, max=2.5, as.groups = TRUE)

7.3 Metody modelowania

7.3.1 Metody 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 Metody modelowania | funkcja eyefit()

v_eye <- eyefit(variog(as.geodata(punkty, 'temp')))
ve_fit <- as.vgm.variomodel(v_eye[[1]])

7.3.3 Metody modelowania | funkcja fit.variogram()

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

7.3.4 Metody modelowania | Liniowy model regionalizacji

W przypadku, gdy analizowane zjawisko jest złożone, odwzorowanie kształtu semiwariogramu empirycznego wymaga połączenia dwóch lub większej liczby modeli podstawowych. W takiej sytuacji konieczne jest spełnienie dwóch warunków:

  • Wszystkie zastosowane modele muszą być dopuszczalne (vgm())
  • Wariancja progowa każdego podstawowego modelu musi być dodatnia

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ć parametry semiwariogramu, tj. semiwariancja cząstkowa (psill), skrócona nazwa używanej funkcji (model) oraz jej zasięg (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 Modelowanie izotropowe | Model nuggetowy

Model nuggetowy (Nug) 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ą. Ten typ modelu najczęściej używany jest w modelach złożonych. W takich wypadkach służy on do określania, między innymi, błędu pomiarowego czy zmienności na krótkich odstępach.

vario <- variogram(temp~1, punkty)
plot(vario)

model_nug <- vgm(psill=10, model='Nug', range=0)
model_nug
##   model psill range
## 1   Nug    10     0
plot(vario, model=model_nug)

fitted_nug <- fit.variogram(vario, model_nug)
fitted_nug
##   model    psill range
## 1   Nug 3.933029     0
plot(vario, model=fitted_nug)

7.4.2 Modelowanie izotropowe | 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.

vario <- variogram(temp~1, punkty)
plot(vario)

model_sph <- vgm(psill=10, model = 'Sph', range=3000)
model_sph
##   model psill range
## 1   Sph    10  3000
plot(vario, model=model_sph)

fitted_sph <- fit.variogram(vario, model_sph)
fitted_sph
##   model   psill    range
## 1   Sph 12.5445 4440.768
plot(vario, model=fitted_sph)

7.4.3 Modelowanie izotropowe | 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.

vario <- variogram(temp~1, punkty)
plot(vario)

model_exp <- vgm(psill=10, model = 'Exp', range=3000)
model_exp
##   model psill range
## 1   Exp    10  3000
plot(vario, model=model_exp)

fitted_exp <- fit.variogram(vario, model_exp)
fitted_exp
##   model    psill    range
## 1   Exp 16.50871 3084.309
plot(vario, model=fitted_exp)

7.4.4 Modelowanie izotropowe | 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.

vario <- variogram(temp~1, punkty)
plot(vario)

model_gau <- vgm(psill=13, model = 'Gau', range=3000)
model_gau
##   model psill range
## 1   Gau    13  3000
plot(vario, model=model_gau)

fitted_gau <- fit.variogram(vario, model_gau)
fitted_gau
##   model    psill    range
## 1   Gau 8.085963 798.7454
plot(vario, model=fitted_gau)

7.4.5 Modelowanie izotropowe | 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.

vario <- variogram(temp~1, punkty)
plot(vario)

model_pow <- vgm(psill=0.03, model = 'Pow', range=0.7)
model_pow
##   model psill range
## 1   Pow  0.03   0.7
plot(vario, model=model_pow)

fitted_pow <- fit.variogram(vario, model_pow)
fitted_pow
##   model      psill     range
## 1   Pow 0.02732273 0.7382382
plot(vario, model=fitted_pow)

7.4.6 Modelowanie izotropowe | 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().

vario <- variogram(temp~1, punkty)
model_zl1 <- vgm(psill=10, model = 'Sph', range = 3000, nugget = 0.5)
model_zl1
##   model psill range
## 1   Nug   0.5     0
## 2   Sph  10.0  3000
plot(vario, model=model_zl1)

fitted_zl1 <- fit.variogram(vario, model_zl1)
fitted_zl1
##   model      psill    range
## 1   Nug  0.7346354    0.000
## 2   Sph 13.2621876 5602.027
plot(vario, model=fitted_zl1)

7.4.7 Modelowanie izotropowe | 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.

vario <- variogram(temp~1, punkty)
model_zl2 <- vgm(10, 'Gau', 3000, add.to = vgm(4, model = 'Gau', range = 500, nugget = 0.5))
model_zl2
##   model psill range
## 1   Nug   0.5     0
## 2   Gau   4.0   500
## 3   Gau  10.0  3000
plot(vario, model=model_zl2)

fitted_zl2 <- fit.variogram(vario, model_zl2)
plot(vario, model=fitted_zl2)

7.5 Modelowanie anizotropowe

7.5.1 Anizotropia

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

  • Zasięg w dominującym kierunku
  • Kąt określający dominujący kierunek
  • Proporcję anizotropii, czyli relację pomiędzy zasięgiem w dominującym kierunku a zasięgiem w przeciwległym 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 głównym kierunku a kierunku prostopadłym. 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).

vario_map <- variogram(temp~1, punkty, cutoff=4000, width=400, map=TRUE)
plot(vario_map)

vario_kier <- variogram(temp~1, punkty, alpha = c(45, 90, 135, 180), cutoff=4000)
plot(vario_kier, plot.numbers=TRUE)

vario_kier_fit <- vgm(psill=8, model='Sph', range=4000, nugget=0.5, anis = c(45, 0.4))
plot(vario_kier, vario_kier_fit, as.table=TRUE)