🌡️ ინტერპოლაცია: IDW და Kriging¶

გვაქვს გაზომვები ცალკეულ წერტილებში (ტემპერატურა მეტეოსადგურებზე, სიმაღლე GPS-წერტილებზე, დაბინძურება სინჯის აღების ადგილებზე), მაგრამ გვჭირდება მნიშვნელობა ნებისმიერ სხვა წერტილში ან მთელი რასტრი. ამას ინტერპოლაცია წყვეტს. მისი ყველა მეთოდი ერთ ვარაუდს ეფუძნება: ახლო წერტილები უფრო ჰგვანან ერთმანეთს, ვიდრე შორი (ტობლერის პირველი კანონი). მეთოდები იმით განსხვავდება, რამდენად და როგორ ითვალისწინებენ ამას.

მეთოდი	იდეა	ძლიერი მხარე	სუსტი მხარე
უახლოესი მეზობელი	მნიშვნელობა = უახლოესი ნიმუშის (Voronoi)	უმარტივესი, მონაცემს არ ცვლის	საფეხურებრივი ზედაპირი
IDW	ყველა ნიმუშის შეწონილი საშუალო, წონა მანძილით მცირდება	მარტივი, სწრაფი	ზედაპირი ნიმუშებს შორის ვერ გასცდება მათ დიაპაზონს, ირგვლივ „ხარის თვალები" ჩნდება
Kriging	წონებს მონაცემთა სივრცული კორელაცია (ვარიოგრამა) განსაზღვრავს	იძლევა ცდომილების შეფასებასაც	უფრო რთულია, ვარიოგრამის არჩევა სჭირდება

🎮 სცადეთ თქვენც¶

დააჭირეთ რუკას ნიმუშის დასამატებლად. ფარულ „ნამდვილ" ზედაპირზე ყოველი მეთოდის ცდომილება ქვემოთ ცხრილშია. გადართეთ მეთოდები და ცადეთ სლაიდერები.

დაკვირვება

IDW-ის ხარისხი p: დიდი p-ზე ზედაპირი უახლოესი მეზობლისკენ მიიწევს (საფეხურები), მცირეზე ძალიან გასწორებულია. p = 2 ტრადიციული საწყისი მნიშვნელობაა.
Kriging-ის დისპერსია (ყვითელი რუკა): ნიმუშთან ახლოს ღია (დაბალი), ნიმუშებისგან შორს მუქი (მაღალი). IDW ამას არ გვეუბნება: ის ყველგან თანაბრად თავდაჯერებულია.
ცოტა ნიმუშით და მცირე რადიუსით Kriging შეიძლება IDW-ზე ცუდიც იყოს. ვარიოგრამა მონაცემებს უნდა ერგებოდეს.

⚙️ IDW: მანძილით შეწონილი საშუალო¶

შესაფასებელ წერტილ P-სთვის ყოველ ნიმუშს i ვაძლევთ წონას wᵢ = 1 / dᵢᵖ, სადაც dᵢ არის მანძილი P-დან i-მდე:

z(P) = Σ wᵢ · zᵢ  /  Σ wᵢ

ზუსტად ნიმუშის წერტილზე d = 0, ამიტომ პირდაპირ ნიმუშის მნიშვნელობას ვაბრუნებთ.
წონების ჯამი ნორმირებულია, ამიტომ შედეგი ყოველთვის ნიმუშების მინიმუმსა და მაქსიმუმს შორისაა: IDW ვერ ივარაუდებს მწვერვალს, რომლის ნიმუშიც არ გვაქვს.

⚙️ Kriging: წონები მონაცემთა სტრუქტურიდან¶

IDW-ში წონა მხოლოდ მანძილზეა დამოკიდებული. Kriging ორ რამეს ითვალისწინებს:

რამდენად განსხვავდება მნიშვნელობები მანძილის ზრდასთან ერთად: ამას აღწერს ვარიოგრამა γ(h). ის იზრდება h-თან ერთად და რადიუსზე a (range) გადადის პლატოზე: ამაზე შორ წერტილებს ერთმანეთზე გავლენა აღარ აქვთ.
ნიმუშები ერთმანეთს როგორ ფარავს: ერთმანეთთან ახლოს მყოფი ორი ნიმუში ერთსა და იმავეს გვეუბნება, ამიტომ თითოეულს ნაკლები წონა ეძლევა (დაჯგუფების ეფექტის გასწორება, როგორც Thiessen-ის საშუალოში).

Ordinary kriging-ში წონები წრფივი სისტემიდან მოდის:

| γ(xᵢ,xⱼ)  1 |   | w |   | γ(xᵢ,P) |
|   1ᵀ      0 | · | μ | = |    1    |

ბოლო მწკრივი ამბობს Σ w = 1, μ არის ლაგრანჟის მამრავლი. შედეგი: z(P) = Σ wᵢ zᵢ და დისპერსია σ²(P) = Σ wᵢ γ(xᵢ,P) + μ.

ჩვენი დემო იყენებს ექსპონენციალურ ვარიოგრამას γ(h) = 1 − e^(−3h/a) ნაგეტის გარეშე. ასეთ შემთხვევაში ზღურბლის (sill) სიდიდე წონებზე გავლენას არ ახდენს, ამიტომ მხოლოდ a რჩება პარამეტრად. რეალურ მონაცემებზე ვარიოგრამას ემპირიული მონაცემებიდან ირჩევენ: ითვლიან γ-ს წერტილთა წყვილებისთვის მანძილის ინტერვალებით და მასზე მოდელს ადებენ (spherical, exponential, gaussian).

სირთულე

ერთ წერტილში Kriging-ს n+1 განტოლებათა სისტემის ამოხსნა სჭირდება, O(n³). მატრიცის ინვერსიას ერთხელ ვთვლით და ყოველი უჯრისთვის მხოლოდ ვექტორზე ვამრავლებთ, O(n²). ამიტომ რასტრის დიდი ზომისთვის ხშირად მხოლოდ უახლოეს k ნიმუშს იღებენ ("moving neighborhood"), რისთვისაც k-d tree გამოიყენება.

🐍 რეალიზაცია სუფთა Python-ით¶

interp.py
import math


def nearest_value(samples, x, y):
    """უახლოესი მეზობელი = Voronoi უჯრის მნიშვნელობა."""
    return min(samples, key=lambda s: (s[0] - x) ** 2 + (s[1] - y) ** 2)[2]


def idw(samples, x, y, power=2.0):
    """Inverse Distance Weighting. samples = [(x, y, მნიშვნელობა), ...]."""
    num = den = 0.0
    for sx, sy, v in samples:
        d2 = (sx - x) ** 2 + (sy - y) ** 2
        if d2 == 0:
            return v                           # ზუსტად ნიმუშის წერტილზე
        w = d2 ** (-power / 2)                 # 1 / d^power
        num += w * v
        den += w
    return num / den


def exp_variogram(h, rng):
    """ექსპონენციალური ვარიოგრამა (ნაგეტი 0, ზღურბლი 1)."""
    return 1 - math.exp(-3 * h / rng)


def solve(A, b):
    """გაუსის მეთოდი მთავარი ელემენტის არჩევით."""
    n = len(A)
    M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for c in range(n):
        p = max(range(c, n), key=lambda r: abs(M[r][c]))
        M[c], M[p] = M[p], M[c]
        for r in range(c + 1, n):
            f = M[r][c] / M[c][c]
            for k in range(c, n + 1):
                M[r][k] -= f * M[c][k]
    x = [0.0] * n
    for r in range(n - 1, -1, -1):
        x[r] = (M[r][n] - sum(M[r][k] * x[k] for k in range(r + 1, n))) / M[r][r]
    return x


def ordinary_kriging(samples, x, y, rng):
    """Ordinary kriging: აბრუნებს (შეფასება, დისპერსია)."""
    n = len(samples)
    A = [[exp_variogram(math.hypot(a[0] - b[0], a[1] - b[1]), rng) for b in samples] + [1.0] for a in samples]
    A.append([1.0] * n + [0.0])
    g0 = [exp_variogram(math.hypot(s[0] - x, s[1] - y), rng) for s in samples] + [1.0]
    sol = solve(A, g0)
    w, mu = sol[:n], sol[n]
    z = sum(wi * s[2] for wi, s in zip(w, samples))
    var = sum(wi * gi for wi, gi in zip(w, g0[:n])) + mu
    return z, var

გამოყენება¶

import math
import random
from interp import *

random.seed(5)


def truth(x, y):
    return 20 + 16 * math.sin(x / 17) * math.cos(y / 11) + 10 * math.sin((x + y) / 26)


samples = [(x, y, truth(x, y)) for x, y in ((random.uniform(0, 100), random.uniform(0, 50)) for _ in range(25))]

q = (40, 20)
print('ნამდვილი: ', round(truth(*q), 2))
print('უახლოესი: ', round(nearest_value(samples, *q), 2))
print('IDW p=2:  ', round(idw(samples, *q, power=2), 2))
z, var = ordinary_kriging(samples, *q, rng=30)
print('Kriging:  ', round(z, 2), f'(დისპერსია {var:.2f})')

ნამდვილი:  24.63
უახლოესი:  30.2
IDW p=2:   25.08
Kriging:   25.66 (დისპერსია 0.79)

შემოწმება

ჩვენი Kriging-ის წონები იგივეა, რაც numpy.linalg.solve-ით იგივე სისტემის ამოხსნით მიიღება (განსხვავება ~10⁻¹⁵), წონების ჯამი 1-ია და ნიმუშის წერტილზე შედეგი ზუსტად ნიმუშის მნიშვნელობაა. JavaScript-ის ვერსია, რომელიც დემოშია, იგივე რიცხვებს იძლევა.

📚 იგივე ბიბლიოთეკით¶

SciPyPyKrigeQGIS

import numpy as np
from scipy.interpolate import griddata

xy = np.array([(s[0], s[1]) for s in samples])
v = np.array([s[2] for s in samples])

nearest = griddata(xy, v, [(40, 20)], method='nearest')
linear = griddata(xy, v, [(40, 20)], method='linear')   # TIN + ბარიცენტრული

method='nearest' ჩვენი nearest_value-ის იდენტურია (შემოწმებულია). 'linear' ააგებს Delaunay-ს და სამკუთხედებში ბარიცენტრულად აინტერპოლირებს, თუმცა ამოზნექილი გარსის გარეთ აბრუნებს nan-ს.

from pykrige.ok import OrdinaryKriging

ok = OrdinaryKriging(xs, ys, vs, variogram_model='exponential')
z, variance = ok.execute('grid', grid_x, grid_y)

PyKrige ვარიოგრამის პარამეტრებს (nugget, sill, range) თავად არჩევს მონაცემებიდან.

import processing

processing.run('gdal:gridinversedistance', {
    'INPUT': points_layer, 'Z_FIELD': 'value', 'POWER': 2,
    'OUTPUT': 'idw.tif'})

რა არის შემოწმებული

SciPy-ს nearest და linear მეთოდებით მუშაობა გავუშვი. PyKrige-ისა და QGIS-ის ფრაგმენტები არ გამიშვია ამ გარემოში, ამიტომ პარამეტრების სახელები დოკუმენტაციასთან შეამოწმეთ.

✍️ დავალებები¶

დავალება 1: ვინ არის უკეთესი?

გამოიყენეთ leave-one-out შემოწმება (cross-validation): თითოეული ნიმუში გამოაკელით, შეაფასეთ იგი დანარჩენებით და დათვალეთ საშუალო კვადრატული ცდომილება (RMSE). შეადარეთ უახლოესი მეზობელი, IDW (p = 1…4) და Kriging (a = 10, 20, 40, 80).

ამოხსნა

def loo(fn):
    err = []
    for i, s in enumerate(samples):
        rest = samples[:i] + samples[i + 1:]
        err.append(fn(rest, s[0], s[1]) - s[2])
    return math.sqrt(sum(e * e for e in err) / len(err))

print('უახლოესი', round(loo(nearest_value), 2))
for p in (1, 2, 3, 4):
    print('IDW p =', p, round(loo(lambda T, x, y: idw(T, x, y, p)), 2))
for a in (10, 20, 40, 80):
    print('Kriging a =', a, round(loo(lambda T, x, y: ordinary_kriging(T, x, y, a)[0]), 2))

უახლოესი 4.46
IDW p = 1 7.41
IDW p = 2 5.16
IDW p = 3 4.15
IDW p = 4 3.75
Kriging a = 10 8.77
Kriging a = 20 6.59
Kriging a = 40 4.83
Kriging a = 80 4.27

დასკვნა. „ყველაზე სწორი მეთოდი" არ არსებობს. ჩვენი ზედაპირი ნელა იცვლება (ტალღის სიგრძე ~100 ერთეული), ამიტომ დიდი რადიუსი (a = 80) პატარაზე გაცილებით უკეთესია, ხოლო IDW p = 4-ზე აქ Kriging-საც კი სჯობს. ვარიოგრამის რადიუსი ზედაპირის ბუნებას უნდა ერგებოდეს: a = 10 ძალიან ცუდია, რადგან ნიმუშებს ერთმანეთზე გავლენას თითქმის არ აძლევს. ამიტომ რეალურ მონაცემებზე a-ს ემპირიული ვარიოგრამიდან ირჩევენ და არა გამოცნობით.

დავალება 2: ექსტრაპოლაცია

რას აბრუნებს თითოეული მეთოდი ნიმუშებიდან ძალიან შორს, წერტილში (300, 300)? შეადარეთ ნიმუშების საშუალოს.

ამოხსნა

q = (300, 300)
print(round(sum(s[2] for s in samples) / len(samples), 2),
      round(nearest_value(samples, *q), 2),
      round(idw(samples, *q, power=2), 2),
      round(ordinary_kriging(samples, *q, rng=30)[0], 2))

20.43 17.89 19.7 21.21

(საშუალო, უახლოესი, IDW, Kriging.) IDW და Kriging შორს მნიშვნელობას ერთგვარ „საშუალოსკენ" წევენ, უახლოესი კი უბრალოდ ყველაზე ახლო ნიმუშს იმეორებს. Kriging-ის შედეგი დაჯგუფების გასწორებული საშუალოა, ამიტომ უბრალო საშუალოს ზუსტად არ ემთხვევა. მთავარი ის არის, რომ ნიმუშებისგან შორს ნებისმიერი მეთოდი მონაცემზე კი არა, დაშვებაზე დგას, ამიტომ ასეთი წერტილების ნდობა არ შეიძლება. შემდეგი დავალება სწორედ ამას აჩვენებს.

დავალება 3: დისპერსია როგორც ნდობა

გამოთვალეთ Kriging-ის დისპერსია: (ა) ნიმუშიდან ძალიან ახლოს, (ბ) რუკის შუაგულში (50, 25), (გ) შორს (300, 300).

ამოხსნა

s0 = samples[0]
for q in [(s0[0] + 1, s0[1]), (50, 25), (300, 300)]:
    print(q, round(ordinary_kriging(samples, *q, rng=30)[1], 3))

(63.29016948897019, 37.08934946303647) 0.179
(50, 25) 0.644
(300, 300) 1.092

დისპერსია ნიმუშთან ახლოს მცირეა, ნიმუშებს შორის საშუალო, შორს კი ზღურბლს (≈1) უახლოვდება: აქ ნიმუშებისგან ინფორმაცია აღარ გვრჩება. ამ რუკის გამოყენებით შეიძლება გადაწყვეტა, სად ღირს დამატებითი გაზომვა.

დავალება 4: SciPy-ს linear და ჩვენი მეთოდები

ბადეზე x = 0, 4, …, 100, y = 0, 4, …, 50 გამოიყენეთ griddata(method='linear') და შეამოწმეთ, რამდენ წერტილში აბრუნებს nan-ს. რატომ?

ამოხსნა

import numpy as np
from scipy.interpolate import griddata

grid = [(x, y) for x in range(0, 101, 4) for y in range(0, 51, 4)]
lin = griddata(xy, v, np.array(grid), method='linear')
print(int(np.isnan(lin).sum()), 'of', len(grid))

92 of 338

27% წერტილი ნიმუშების ამოზნექილი გარსის (convex hull) გარეთაა და სამკუთხედი მათ ვერ ფარავს. IDW და Kriging ასეთ წერტილებშიც აბრუნებენ მნიშვნელობას, მაგრამ ის ექსტრაპოლაციაა (იხ. დავალება 2). ამიტომ რუკის კიდეები სიფრთხილით უნდა წავიკითხოთ.