1
Geostatistické metody interpolace ­ STRUKTURNÍ ANALÝZA A KRIGING
Kriging patří do tzv. stochastických metod. Doposud uvedené metody se naopak
označují jako deterministické a žádná z nich dosud neřešila následující problémy:
* počet bodů nutných k výpočtu lokálního průměru (interpolované hodnoty)
* velikost, orientaci a tvar okolí
* zda neexistuje jiná cesta k definování vah než funkce vzdálenosti bodů
* jaké jsou chyby a nejistoty spojené s interpolovanými hodnotami
Odpovědi na tyto otázky poskytují geostatistické postupy založené na tzv. strukturní
analýze a krigingu.
V případě krigingu je interpolovaný povrch tvořen ze tří složek (obr. 1). První složku
představuje tzv. strukturální komponenta s konstantním průměrem či trendem obecný
trend (tzv. drift). Druhou složku povrchu představují kolísání (drobné
sníženiny či vyvýšeniny), jejichž podstatu lze vyjádřit určitou matematickou funkcí
jako v případě trendu, ale která vyjadřují určitou prostorovou korelaci (tzv.
regionalizovaná proměnná) ­ právě tato složka je geostatistickými metodami
odhadována. Třetí složku potom představuje náhodný šum.
Obr. 1. Základní komponenty spojitého povrchu (i ­ trendová složka ­ drift; ii ­ tzv.
regionalizovaná proměnná; iii ­ náhodná složka)
Hodnota náhodné proměnné Z v bodě x bude
( ) '',
)()(  ++= xxxZ
x - pozice v 1, 2 či 3 rozměrném prostoru
Z - interpolovaná proměnná
Z(x) - hodnota proměnné v bodě x
(x) - deterministická složka (trend)
'(x) - stochastická složka (regionalizovaná proměnná) - lokálně proměnné, ale
prostorově závislé reziduum od (x)
'' - náhodná, prostorově nezávislá složka, gaussovský šum s nulovým průměrem a
s rozptylem 2
. Velké písmeno Z značí, že se jedná o náhodnou funkci a ne o
měřenou hodnotu proměnné z.
Regionalizovaná proměnná se odhaduje v tzv. strukturní analýze pomocí hodnot
strukturní funkce - tzv. semivariance (h) ­ je to míra nepodobnosti hodnot dvou
bodů vzdálených od sebe hodnotu h.
( ) ''
)()(  ++= hxxZ
2
1. Strukturní analýza
Je to procedura zahrnující výpočet strukturálních funkcí a je výchozím krokem
geostatistického modelování.
Obr. 2 Příklad výpočtu měr prostorové variability pro 1D (řadu hodnot)
Ke kvantifikaci prostorové autokorelace, která vyjadřuje skutečnost, že objekty blízké
si jsou více podobné než objekty vzdálenější slouží strukturální funkce. Strukturální
funkce měří sílu korelačního vztahu jako funkci vzdálenosti. Na obr. 2 je pro
jednoduchost odvozena strukturální funkce pro řadu pravidelně rozmístěných bodů
na linii (tedy pro 1D). Charakteristiky, které popisují prostorovou variabilitu lze odvodit
z měr úrovně a variability následovně:
průměr = (1+3+6+5+3+1+2+3)/8=3,0
rozptyl=[(1-3)2
+(3-3)2
+(6-3)2
+(5-3)2
+(3-3)2
+(1-3)2
+(2-3)2
+(2-3)2
]/8=2,75
kovariance(1)=[(1-3)*(3-3)+(3-3)*(6-3)+(6-3)*(5-3)+(5-3)*(3-3)+(3-3)*(1-3)+(1-3)*(2-3)+(2-3)*(3-
3)]/7=1,14
semivariance(1)=[(1-3)2
+(3-6)2
+(6-5)2
+(5-3)2
+(3-1)2
+(1-2)2
+(2-3)2
]/7=3,43
semivariance(2)=[(1-6)2
+(3-5)2
+(6-3)2
+(5-1)2
+(3-2)2
+(1-3)2
]/6=9,83
semivariance(3)=[(1-5)2
+(3-3)2
+(6-1)2
+(5-2)2
+(3-3)2
]/5=12,50
...
Abychom popsali závislost hodnot studovaného jevu v prostoru, vyneseme hodnoty
semivariancí pro všechny dvojice bodů do semivariogramu obdobně jako ve výše
uvedeném případě. Semivariogram je strukturální funkce, která popisuje závislost
průměrné kvadratické diference hodnot prostorové proměnné veličiny Z na
vzdálenosti h. Semivarianci lze odhadnout z naměřených dat podle následujícího
vztahu:
( ) ( ) ( )( )=
+-=
n
i
ii hxzxz
n
h
1
2
2
1
^ (1)
kde:
n - počet dvojic bodů pozorování proměnné s atributem z vzdálených o hodnotu h
h - tzv. lag - vzdálenost dané dvojice bodů.
Obr. 4 Experimentální semivariogram (+) s charakteristickými hodnotami pro
vzdálenosti h (ˇ) a proložený teoretický model semivariogramu (plná čára)
3
Graf hodnot ( )h^ a h se označuje jako experimentální semivariogram a je prvním
krokem ke kvantitativní deskripci regionalizované proměnné. I když je semivariogram
vektorová funkce, sestavuje se často jako izotropní (tj. bez ohledu na směr) pro
celkové charakterizování hodnoceného náhodného pole nebo tehdy, je-li k dispozici
omezený počet pozorování. Experimentálním semivariogramem se v následujícím
kroku prokládá teoretický model.
Prvky semivariogramu
Na obr. 5 uveden často používaný tzv. sférický model s vysvětlením používané
terminologie.
Obr.5 Příklad teoretického semivariogramu ­ sférický model. Parametry
semivariogramu: a - dosah (range), d ­ rozpětí, c0 - zbytkový rozptyl (nugget), c=c0 +
c1 - práh (sill), h ­ lag (krok vzdálenosti)
Na horizontální ose je vynášena vzdálenost h mezi jednotlivými vstupními body
interpolovaného povrchu (tzv. lag), na vertikální ose potom rozptyl zkoumané
proměnné jako funkce vzájemných vzdáleností jednotlivých měřených bodů ­ tedy
semivariance. Je mírou, která vyjadřuje, jak velké je okolí, daného bodu, ve kterém
se nacházejí body sousední, jejichž hodnota interpolovaného atributu závisí
(koreluje) s hodnotou v tomto bodě. Takto vynesenými body je proložena křivka
mající charakteristický tvar. Je-li vzdálenost mezi dvěma body malá, jejich hodnoty
jsou podobné a hodnota semivariance je také malá. Se zvětšující se vzdáleností
hodnota semivariance roste. Při určité vzdálenosti dvou bodů je možno říci, že jejich
hodnoty (např. výšky) spolu již nekorelují a hodnota semivariance i se zvyšováním
vzdálenosti již neroste, ale zůstává konstantní. Plochá část semivariogramu určuje
tzv. práh (sill) a je rovna rozptylu zpracovávaných dat. Existence prahu značí, že se
zvětšující se vzdáleností se hodnoty semivariance nemění. Kritická hodnota
vzdálenosti, na níž se křivka semivariogramu stává rovnoběžnou s vodorovnou osou
se označuje jako dosah (range). Dosah definuje pro daný bod velikost okolí, které je
nutné uvažovat při interpolaci hodnoty v daném bodě.
Z výše uvedeného by také mělo platit, že proložená křivka semivariogramu by měla
procházet počátkem souřadné soustavy (nulová vzdálenost mezi body znamená
zákonitě také nulovou hodnotu rozptylu).Velmi často nenabývají experimentální
semivariogramy v počátku nulové hodnoty; protínají osu y v nenulové hodnotě, která
je nazývána zbytkový rozptyl (nugget variance). To může ukazovat na rozptyl
menší než je "vzorkovací" vzdálenost, nebo na malou přesnost měření, kdy např.
jsou v datech obsaženy dva vzorky ze stejného místa, pokaždé s jinou hodnotou.
Zbytkový rozptyl je tak vyjádřením náhodného šumu '' .V případě nulového
zbytkového rozptylu je krigování exaktním interpolátorem.
Efekt anizotropie - rozdíly v hodnotách semivariance v závislosti na směru
4
Teoretický semivariogram - model, který nejlépe aproximuje průběh
experimentálního semivariogramu v okolí počátku a prahu, proces hledání
teoretického semivariagramu se někdy označuje jako strukturní analýza.
Příklad teoretického semivariogramu: Sférický model
Obr. 9 Sférický model semivariogramu














-+=
3
10
2
1
2
3
)(
a
h
a
h
cch ........... pro ah 
10)( cch += ........... pro ah >
Lineární model s prahem - jednoduchý a poměrně často využívaný zvláště
programy provádějícími interpolaci pomocí krigování na základě automaticky
vypočítaného a vyhodnoceného semivariogramu. Při provádění strukturální analýzy
se využívá raději jiných přechodových modelů.
Obr. 13 Lineární model semivariogramus prahem
bhch += 0)( ........... pro ah 
10)( cch += ........... pro ah >
5
2. KRIGING ­ geostatistické metody interpolace
Interpolovaná odhadovaná hodnota je vypočtena jako lineární kombinace vstupních
hodnot:
=
=
n
i
ii xzxz
1
0 )()(^  (2)
kde pro váhy platí
=
=
n
i
i
1
1
K určení vah i potřebných pro interpolaci se využívá strukturní analýzy semivariogramu.
Váhy i jsou zvoleny tak, aby odhad )(^ 0xz byl nestranný a odhad
rozptylu 2
e byl menší, než jakákoliv jiná lineární kombinace pozorovaných hodnot
(minimální). Pro minimální rozptyl hodnot )(^ 0xz platí výraz :
=
+=
n
i
iie xx
1
0
2
),(^  (3)
kde:
=
=+
n
i
jjii xxxx
1
0 ),(),(  pro všechna j. (4)
Hodnota ),( ji xx je semivariance proměnné z mezi body xi a xj. Hodnota ),( 0xxi je
semivariance proměnné z mezi bodem xi a bodem x0, pro který je hodnota proměnné
z zjišťována. Obě hodnoty lze získat z vhodného teoretického modelu
semivariogramu. Hodnota  je tzv. Lagrangeův multiplikátor, který zajišťuje
požadavek minimalizace odchylek a zároveň podmínku, že suma vah je rovna jedné.
Uvedená metoda se označuje jako základní (ordinary) kriging a je možné ji použít
pro interpolaci v pravidelné mřížce hodnot, ke konstrukci map (např. izolinií).
6
PŘÍKLAD: Výpočet neznámé hodnoty v bodě metodou základního krigingu.
Na základě změřených hodnot veličiny Z v pěti bodech (i=1,..., 5) (viz. obrázek) máme za
úkol odhadnout hodnotu Z bodě (i=0) o souřadnicích (x=5, y=5) metodou krigingu.
Obr. 1 Vstupní data pro lokální odhad metodou základního krigování (podtržená čísla značí
hodnotu atributu v bodě)
Na základě předem provedené strukturní analýzy použijeme sférický semivariogram.














-+=
3
10
2
1
2
3
)(
a
h
a
h
cch ........... pro ah 
10)( cch += ........... pro ah >
Parametry použitého teoretického semivariogramu odhadnuté z předchozí strukturní analýzy
jsou následující:
c0 = 2,5
c1 = 7,5
a = 10,0 (dosah)
Data v pěti měřených bodech mají následující souřadnice
i x y z
1 2 2 3
2 3 7 4
3 9 9 2
4 6 5 4
5 5 3 6
Pokud budeme dále značít:
A ­ matice semivariancí mezi všemi dvojicemi bodů
b ­ vektor semivariancí mezi všemi body a bodem predikovaným
 ­ vektor vah jednotlivých bodů
 ­ tzv. Lagrangeův člen
potom základní vztah pro odhad metodou korigování (rovnice 4) lze psát jako:
7
bA =
Pro vlastní řešení je nutné vypočítat váhy , které musí splňovat podmínku 1=
Uvedený základní vztah lze vyjádřit jako soustavu rovnic:






















=












































1
.
.
.
.
.
.
*.
01...11
1...
...
...
...
1...
1...
0
20
10
2
1
21
22221
11211
nnnnnn
n
n










V tomto zápisu poslední řádek a poslední sloupec v první matici a hodnota Lagrangeova
členu  jsou použity pro zajištění podmínky sumy vah 1= . Hodnota Lagrangeova
multiplikátoru  také slouží pro výpočet rozptylu odhadnuté hodnoty. Uvedená soustava
rovnic nám poskytne hodnoty všech vah  a hodnotu . V maticovém zápisu lze tedy psát:






=-


bA 1
Aby bylo možné vyčíslit hodnoty semivariancí, je v prvním kroku zapotřebí vytvořit matici
vzdáleností mezi datovými body:
i 1 2 3 4 5
1 0,000 5,099 9,899 5,000 3,162
2 5,099 0,000 6,325 3,606 4,472
3 9,899 6,325 0,000 5,000 7,211
4 5,000 3,606 5,000 0,000 2,236
5 3,162 4,472 7,211 2,236 0,000
Vektor vzdáleností mezi měřenými body a bodem predikovaným:
i 0
1 4,234
2 2,828
3 5,657
4 1,000
5 2,000
Těchto vzdáleností využijeme k výpočtu semivariancí pro sférický model semivariogramu
s výše uvedenými parametry c0 , c1 , a ­ tedy k sestavení matice A a vektoru b:
Matice A:
i 1 2 3 4 5 6
1 2,500 7,739 9,999 7,656 5,939 1,000
2 7,739 2,500 8,667 6,381 7,196 1,000
3 9,999 8,667 2,500 7,656 9,206 1,000
4 7,656 6,381 7,656 2,500 4,936 1,000
5 5,939 7,196 9,206 4,936 2,500 1,000
6 1,000 1,000 1,000 1,000 1,000 1,000
Ve výše uvedené matici má řádek navíc (i=6) zajistit podmínku, že váhy budou mít sumu
rovnu jedné.
8
Vektor b:
i 0
1 7,151
2 5,597
3 8,815
4 3,621
5 4,720
6 1,000
Inverzní matce A
-1
:
i 1 2 3 4 5 6
1 -,172 ,050 ,022 -,026 ,126 ,273
2 ,050 -,167 ,032 ,077 ,007 ,207
3 ,022 ,032 -,111 ,066 -,010 ,357
4 -,026 ,077 ,066 -,307 ,190 ,030
5 ,126 ,007 -,010 ,190 -,313 ,134
6 ,273 ,207 ,357 ,003 ,134 -6,873
Řešením výše uvedené soustavy rovnic lze pro jednotlivé vzdálenosti získat hodnoty vah :
i 
1 0,0175
2 0,2281
3 -0,0891 vypočtené hodnoty vah
4 0,6437
5 0,1998
6 0,1182 vypočtená hodnota 
Pro váhy i=1,...5 platí, že jejich suma se rovná jedné, v posledním řádku je hodnota
Lagrangeova členu .
Vzdálenosti měřených bodů od bodu predikovaného, již přísluší výše určené váhy:
i 0
1 4,234
2 2,828
3 5,657
4 1,000
5 2,000
Podle rovnice (2) bude odhad hodnoty Z v bodě (i=0) o souřadnicích (x=5, y=5):
Z(xi=0) = 0,0175*3+0,2281*4-0,0891*2+0,6437*4+0,1998*6 =
Z(xi=0) =4,560
Rozptyl odhadu (rovnice 3):
e
2
= [0,0175*7,151+0,2281*5,597-0,0891*8,815+0,6437*3621+0,1998*4,720]+  =
e
2
= 3,890 + 0,1182 =
e
2
= 4,008
9
Typy krigování
Na rozdíl od deterministických metod interpolace nabízí metody krigingu vedle
odhadů vlastní interpolované hodnoty také odhady pravděpodobnosti výskytu
těchto hodnot a dále odhady chyb predikce. Existuje celá řada typů krigování, pro
jejich charakterizování předpokládejme jednoduchý model:
( ) )()( iii xxxZ  +=
kde Z(xi) je proměnná v bodě xi, která se skládá z deterministické hodnoty trendu
(xi) a autokorelované náhodné proměnné (xi).
Základní krigování (ordinary kriging) -  je neznámá hodnota, kterou odhadujeme
z měřených dat
Obr. Princip základního krigování
Univerzální krigování (Universal kriging) - (x) je deterministická funkce - např.
polynom druhého stupně jako na přiloženém obrázku.
Obr. Princip univerzálního krigování
Co-kriging
Máme dvě proměnné z1 a z2, které vykazují prostorovou korelaci. Pak lze využít
hodnot proměnné z2, k interpolaci hodnot proměnné z1. Tento koncept je vhodný
zvláště v případech, kdy je proměnná z2 snáze získatelná a rozšiřitelný i na více než
dvě proměnné. Přitom pro přesnější odhady se používá jak autokorelace jednotlivých
proměnných, tak vzájemné (cross) korelace všech použitých proměnných. Základní
co-kriging využívá následujících modelů:
( ) )(111 xxZ  +=
( ) )(222 xxZ  +=
10
kde 1 a 2 jsou neznámé konstanty. Dále dostáváme dvě náhodné chyby 1(x) a
2(x). Základní co-kriging odhaduje hodnotu proměnné Z1(x0) stejně jako základní
krigování, ovšem navíc využívá kovariance s hodnotu Z2(x).
Obr. Princip co-krigingu
Hodnocení a verifikace modelů
Krigování jako interpolační metoda umožňuje pro každý interpolovaný bod
odhadnout potenciální velikost chyby odhadu. Vedle map predikovaných hodnot tak
lze především konstruovat mapy hodnot 2
e (rozptyl krigingu), které vypovídají o
spolehlivosti interpolovaných hodnot. Tyto hodnoty se obvykle prezentují v podobě
map druhé odmocniny 2
e - tzv. směrodatné chyby (odchylky) krigingu (Standard
error map), protože tyto mají stejné jednotky jako predikované hodnoty.
Vyjdeme-li z výše uvedeného příkladu, kdy rozptyl odhadu je e
2
= 4,008. Potom
směrodatná chyba krigingu bude e = 2,002. Budeme-li předpokládat, že chyby
predikce mají normální rozdělení, potom 95% interval spolehlivosti predikovaných
hodnot lze určit z následujícího vztahu:
2
0 96,1)( exZ 
kde Z(x0) je odhad hodnoty proměnné z v bodě x0 a e
2
je rozptyl odhadu. V našem
případě tedy při opakovaném použití stejného modelu padne 95 % odhadovaných
hodnot do intervalu (4,560  1,96*2,002) tj. (0,64;8,48)
Křížová validace modelu - k vytvoření spojitého povrchu jsou použita všechna
vstupní data v měřených bodech. Poté jsou jednotlivé body měření (červené) po
jednom postupně vynechány ze vstupní množiny dat a ze zbývajících (modrých) je
vypočtena hodnota v místě vynechaného bodu.
Obr. Princip křížové validace modelu
11
Statistické zhodnocení
Pozorované a vypočtené hodnoty mohou být porovnávány jednoduchými
charakteristikami (měrami):Je-li )(^
ixZ predikovaná hodnota pro daný bod xi, kterou
obdržíme v procesu křížové validace, potom:
MPE ­ mean prediction error - průměr rozdílů měřených a předikovaných hodnot hodnoty
chyb odhadů by měly být nestranné ­ tedy jejich průměr by se měl rovnat
nule.
n
xzxZ
MPE
n
i
ii=
=
1
))()(^(
RMSPE (root mean square prediction error) ­ druhá odmocnina průměrného
čtverce vzdálenosti vypočtených hodnot (červené body) od teoretických (zelená
přímka v grafech). Tato hodnota slouží k porovnání několika různých modelů. Čím
menší je RMSPE, tím vhodnější je model (tím bližší jsou vypočtené hodnoty
hodnotám měřeným).
n
xzxZ
RMSPE
n
i
ii=
=
1
2
))()(^(
Validace modelu - vstupní soubor měřených hodnot rozdělí na dvě části ­ data
trénovací a testovací. Trénovací množina dat se použije pro odhad trendu a
autokorelačního modelu. Pokud sestavený model vyhovuje trénovacím datům, je
ověřen na datech testovacích.
Obecnou vlastností krigingu jako interpolační metody je podhodnocení vysokých
hodnot a naopak nadhodnocení hodnot nízkých. Tato vlastnost se projeví menší
hodnotou směrnice přímky proložené korelačním polem.
Obr. Korelační pole měřených a predikovaných hodnot
V případě značného podílu šumové složky (např. v důsledku chyb v měření) či v
případě značně komplexního povrchu nedává kriging lepší výsledky než jiné
interpolátory. Na rozdíl o jiných metod kriging nabízí objektivní, a priori metodu
odhadu vhodného okolí pro vlastní interpolaci. Řeší tedy otázku počtu bodů v okolí
daného bodu, otázku velikosti a tvaru tohoto okolí. V případě existence bariér
(náhlých skoků v hodnotách interpolovaného povrchu) nedává kriging dobré výsledky
a je nutné jej rozdělit na elementární části neobsahující bariéry.
Doplňující informace k prostorové analýze a metodám interpolace:
http://www.spatialanalysisonline.com/output/
http://uncert.mines.edu/tutor/
http://www.u.arizona.edu/~donaldm/homepage/glossary.html