克里金(Kriging)插值的原理與公式推導

時間 2019-12-10

標籤克里 kriging 插值原理公式推導简体版

原文原文鏈接

轉自：http://xg1990.com/blog/archives/222ide

學過空間插值的人都知道克里金插值，可是它的變種繁多、公式複雜，還有個半方差函數讓人不知所云
本文講簡單介紹基本克里金插值的原理，及其推理過程。函數

0.引言——從反距離插值(IDW)提及

空間插值問題，就是在已知空間上若干離散點 (xi,yi) 的某一屬性(如氣溫，海拔)的觀測值 zi=z(xi,yi) 的條件下，估計空間上任意一點 (x,y) 的屬性值的問題。優化

直觀來說，根據地理學第必定律，atom

All attribute values on a geographic surface are related to each other, but closer values are more strongly related than are more distant ones.spa

大意就是，地理屬性有空間相關性，相近的事物會更類似。由此人們發明了反距離插值，對於空間上任意一點 (x,y) 的屬性 z=z(x,y) ，定義反距離插值公式估計量code

z^= \sum i = 0 n 1 d α z i

其中 α 一般取1或者2。orm

即，用空間上全部已知點的數據加權求和來估計未知點的值，權重取決於距離的倒數（或者倒數的平方）。那麼，距離近的點，權重就大；距離遠的點，權重就小。blog

反距離插值能夠有效的基於地理學第必定律估計屬性值空間分佈，但仍然存在不少問題：ci

α 的值不肯定
用倒數函數來描述空間關聯程度不夠準確

所以更加準確的克里金插值方法被提出來了數學

1.克里金插值的定義

相比反距離插值，克里金插值公式更加抽象

z o^= \sum i = 0 n λ i z i

其中 zo^ 是點 (xo,yo) 處的估計值，即 zo=z(xo,yo) 。

這裏的 λi 是權重係數。它一樣是用空間上全部已知點的數據加權求和來估計未知點的值。但權重係數並不是距離的倒數，而是可以知足點 (xo,yo) 處的估計值 zo^ 與真實值 zo 的差最小的一套最優係數，即

min λ i V a r (z o^- z o)

同時知足無偏估計的條件

E (z o^- z o) = 0

2.假設條件

不一樣的克里金插值方法的主要差別就是假設條件不一樣。本文僅介紹普通克里金插值的假設條件與應用。

普通克里金插值的假設條件爲，空間屬性 z 是均一的。對於空間任意一點 (x,y) ，都有一樣的指望c與方差 σ2 。

即對任意點 (x,y) 都有

E [z (x, y)] = E [z] = c

V a r [z (x, y)] = σ 2

換一種說法：任意一點處的值 z(x,y) ，都由區域平均值 c 和該點的隨機誤差 R(x,y) 組成，即

z (x, y) = E [z (x, y)] + R (x, y)] = c + R (x, y)

其中 R(x,y) 表示點 (x,y) 處的誤差，其方差均爲常數

V a r [R (x, y)] = σ 2

3.無偏約束條件

先分析無偏估計條件 E(zo^−zo)=0 ，將 zo^=∑ni=0λizi 帶入則有

E (\sum i = 0 n λ i z i - z o) = 0

又由於對任意的z都有 E[z]=c ，則

c \sum i = 0 n λ i - c = 0

即

\sum i = 0 n λ i = 1

這是 λi 的約束條件之一。

4.優化目標/代價函數J

再分析估計偏差 Var(zo^−zo) 。爲方便公式推理，用符號 J 表示，即

J = V a r (z o^- z o)

則有

J = V a r (\sum n i = 0 λ i z i - z o) = V a r (\sum n i = 0 λ i z i) - 2 C o v (\sum n i = 0 λ i z i, z o) + C o v (z o, z o) = \sum n i = 0 \sum n j = 0 λ i λ j C o v (z i, z j) - 2 \sum n i = 0 λ i C o v (z i, z o) + C o v (z o, z o)

爲簡化描述，定義符號 Cij=Cov(zi,zj)=Cov(Ri,Rj) ，這裏 Ri=zi−c ，即點 (xi,yi) 處的屬性值相對於區域平均屬性值的誤差。

則有

J = \sum i = 0 n \sum j n λ i λ j C i j - 2 \sum i = 0 n λ i C i o + C o o

5.代價函數的最優解

再定義半方差函數 rij=σ2−Cij ，帶入J中，有

J = \sum n i = 0 \sum n j = 0 λ i λ j (σ 2 - r i j) - 2 \sum n i = 0 λ i (σ 2 - r i o) + σ 2 - r o o = \sum n i = 0 \sum n j = 0 λ i λ j (σ 2) - \sum n i = 0 \sum n j = 0 λ i λ j (r i j) - 2 \sum n i = 0 λ i (σ 2) + 2 \sum n i = 0 λ i (r i o) + σ 2 - r o o

考慮到 ∑ni=0λi=1

J = σ 2 - \sum n i = 0 \sum n j λ i λ j (r i j) - 2 σ 2 + 2 \sum n i = 0 λ i (r i o) + σ 2 - r o o = 2 \sum n i = 0 λ i (r i o) - \sum n i = 0 \sum n j = 0 λ i λ j (r i j) - r o o

咱們的目標是尋找使J最小的一組 λi ，且J是 λi 的函數，所以直接將J對 λi 求偏導數令其爲0便可。即

\partial J \partial λ i = 0; i = 1, 2, \dots, n

可是要注意的是，咱們要保證求解出來的最優 λi 知足公式 ∑ni=0λi=1 ，這是一個帶約束條件的最優化問題。使用拉格朗日乘數法求解，求解方法爲構造一個新的目標函數

J + ϕ (\sum i = 0 n λ i - 1)

其中 ϕ 是拉格朗日乘數。求解使這個代價函數最小的參數集 ϕ,λ1,λ2,⋯,λn ，則能知足其在 ∑ni=0λi=1 約束下最小化 J 。即

⎧ ⎩ ⎨ ⎪ ⎪ \partial ( J + ϕ ( \sum n i = 0 λ i - 1 ) ) \partial λ i \partial ( J + ϕ ( \sum n i = 0 λ i - 1 ) ) \partial ϕ = 0; i = 1, 2, \dots, n = 0

⎧ ⎩ ⎨ ⎪ ⎪ \partial ( 2 \sum n i = 0 λ i ( r i o ) - \sum n i = 0 \sum n j λ i λ j ( r i j ) - r o o ) \partial λ i \partial ( \partial 2 \sum n i = 0 λ i ( r i o ) - \sum n i = 0 \sum n j λ i λ j ( r i j ) - r o o ) \partial ϕ = 0; i = 1, 2, \dots, n = 0

{2 r i o - \sum n j = 1 (r i j + r j i) λ j \sum n i = 0 λ i = 0; i = 1, 2, \dots, n = 1

因爲 Cij=Cov(zi,zj)=Cji ，所以一樣地 rij=rji ，那麼有

{r i o - \sum n j = 1 r i j λ j \sum n i = 0 λ i = 0; i = 1, 2, \dots, n = 1

式子中半方差函數 rij 十分重要，最後會詳細解釋其計算與定義

在以上計算中咱們獲得了對於求解權重係數 λj 的方程組。寫成線性方程組的形式就是：

⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ r 11 λ 1 + r 12 λ 2 + \dots + r 1 n λ n + ϕ r 21 λ 1 + r 22 λ 2 + \dots + r 2 n λ n + ϕ r n 1 λ 1 + r n 2 λ 2 + \dots + r n n λ n + ϕ λ 1 + λ 2 + \dots + λ n = r 1 o = r 2 o \dots = r n o = 1 (1)

寫成矩陣形式即爲

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ r 11 r 21 \dots r n 1 1 r 12 r 22 \dots r n 2 1 \dots \dots \dots \dots \dots r 1 n r 2 n \dots r n n 1 11 \dots 10 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ λ 1 λ 2 \dots λ n 0 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ r 1 o r 2 o \dots r n o 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

對矩陣求逆便可求解。

惟一未知的就是上文中定義的半方差函數 rij ，接下來將詳細討論

6.半方差函數

上文中對半方差函數的定義爲

r i j = σ 2 - C i j

其等價形式爲

r i j = 1 2 E [(z i - z j) 2]

這也是半方差函數名稱的來由，接下來證實這兩者是等價的：

根據上文定義 Ri=zi−c ，有 zi−zj=Ri−Rj ，則

r i j = 1 2 E [(R i - R j) 2] = 1 2 E [R 2 i - 2 R i R j + R 2 j] = 1 2 E [R 2 i] + 1 2 E [R 2 j] - E [R i R j]

又由於：

E [R 2 i] = E [R 2 j] = E [(z i - c) 2] = V a r (z i) = σ 2

E [R i R j] = E [(z i - c) (z j - c)] = C o v (z i, z j) = C i j

因而有

r i j = 1 2 E [(z i - z j) 2] = 1 2 E [R 2 i] + 1 2 E [R 2 j] - E [R i R j] = 1 2 σ 2 + 1 2 σ 2 - C i j = σ 2 - C i j

σ2−Cij=12E[(zi−zj)2] 得證，如今的問題就是如何計算

r i j = 1 2 E [(z i - z j) 2]

這時須要用到地理學第必定律，空間上相近的屬性相近。 rij=12(zi−zj)2 表達了屬性的類似度；空間的類似度就用距離來表達，定義i與j之間的幾何距離

d i j = d (z i, z j) = d ((x i, y i), (x j, y j)) = (x i - x j) 2 + (y i - y j) 2 - - - - - - - - - - - - - - - - - - \sqrt

克里金插值假設 rij 與 dij 存在着函數關係，這種函數關係能夠是線性、二次函數、指數、對數關係。爲了確認這種關係，咱們須要首先對觀測數據集

{z (x 1, y 1), z (x 2, y 2), z (x 3, y 3), \dots, z (x n - 1, y n - 1), z (x n, y n)}

計算任意兩個點的距離 dij=(xi−xj)2+(yi−yj)2−−−−−−−−−−−−−−−−−−√ 和半方差 σ2−Cij=12E[(zi−zj)2] ，這時會獲得 n2 個

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。