普通克里金 (OK)插值推导
前言
之前做毕设的时候看文献一知半解,现在找个时间补回来。
1.普通克里金推导
1.1定义
克里金插值公式形如:
z
o
^
=
∑
i
=
1
n
λ
i
z
i
\hat{z_o} = \sum^{n}_{i=1}{\lambda_iz_i}
zo^=i=1∑nλizi
z
o
^
\hat{z_o}
zo^ 为预测点估计值;
λ
i
\lambda_i
λi 为权重系数;
z
i
z_i
zi 已知采样点数据;
同时
λ
i
\lambda_i
λi 为满足估计值
z
o
^
\hat{z_o}
zo^与真实值
z
o
z_o
zo的差最小的一套最优系数
min
λ
i
D
(
z
o
^
−
z
o
)
\min_{\lambda_i} D(\hat{z_o}-z_o)
λiminD(zo^−zo)
z
o
^
\hat{z_o}
zo^满足
E
(
z
o
^
−
z
o
)
=
0
E(\hat{z_o}-z_o)=0
E(zo^−zo)=0
1.2假设条件
普通克里金插值的假设条件为,空间属性
z
z
z是均一的。对于空间任意一点
(
x
,
y
)
(x,y)
(x,y),都有同样的期望
c
c
c与方差
σ
2
σ^2
σ2。 即对任意点
(
x
,
y
)
(x,y)
(x,y)都有
E
[
z
(
x
,
y
)
]
=
E
[
z
]
=
c
E[z(x,y)] = E[z] = c
E[z(x,y)]=E[z]=c
D
[
z
(
x
,
y
)
]
=
σ
2
D[z(x,y)] = \sigma^2
D[z(x,y)]=σ2
换一种说法:任意一点处的值
z
(
x
,
y
)
z(x,y)
z(x,y),都由区域平均值
c
c
c和该点的随机偏差
R
(
x
,
y
)
R(x,y)
R(x,y)组成,即
z
(
x
,
y
)
=
E
[
z
(
x
,
y
)
]
+
R
(
x
,
y
)
]
=
c
+
R
(
x
,
y
)
z(x,y)=E[z(x,y)] + R(x,y)] = c + R(x,y)
z(x,y)=E[z(x,y)]+R(x,y)]=c+R(x,y)
其中
R
(
x
,
y
)
R(x,y)
R(x,y)表示点
(
x
,
y
)
(x,y)
(x,y)处的偏差,其方差均为常数
D [ R ( x , y ) ] = σ 2 D[R(x,y)] = \sigma^2 D[R(x,y)]=σ2
1.3无偏约束条件
先分析无偏估计条件 E ( z o ^ − z o ) = 0 E(\hat{z_o}-z_o)=0 E(zo^−zo)=0,将 z o ^ = ∑ i = 1 n λ i z i \hat{z_o} = \sum^{n}_{i=1}{\lambda_iz_i} zo^=∑i=1nλizi带入则有
E
(
∑
i
=
1
n
λ
i
z
i
−
z
o
)
=
0
E(\sum^{n}_{i=1}{\lambda_iz_i}- z_o)=0
E(i=1∑nλizi−zo)=0
又因为对任意的
z
z
z都有
E
[
z
]
=
c
E[z]=c
E[z]=c,则
c
∑
i
=
1
n
λ
i
−
c
=
0
c \sum^{n}_{i=1}{\lambda_i}- c=0
ci=1∑nλi−c=0
即
∑
i
=
1
n
λ
i
=
1
\sum^{n}_{i=1}{\lambda_i} = 1
i=1∑nλi=1
这是
λ
i
λ_i
λi的约束条件之一。
1.4优化目标/代价函数J
再分析估计误差 D ( z o ^ − z o ) D(\hat{z_o}−z_o) D(zo^−zo)。为方便公式推理,用符号 J J J表示,即
J
=
D
(
z
o
^
−
z
o
)
J = D(\hat{z_o}-z_o)
J=D(zo^−zo)
则有
J = D ( ∑ i = 1 n λ i z i – z o ) = D ( ∑ i = 1 n λ i z i ) – 2 C o v ( ∑ i = 1 n λ i z i , z o ) + C o v ( z o , z o ) = ∑ i = 1 n ∑ j = 0 n λ i λ j C o v ( z i , z j ) – 2 ∑ i = 1 n λ i C o v ( z i , z o ) + C o v ( z o , z o ) \begin{aligned} J &=D(\sum^{n}_{i=1}{\lambda_iz_i} – z_o) \\&=D(\sum^{n}_{i=1}{\lambda_iz_i}) – 2 Cov(\sum^{n}_{i=1}{\lambda_iz_i}, z_o) + Cov(z_o, z_o) \\&= \sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_jCov( z_i, z_j)} – 2 \sum^{n}_{i=1}{\lambda_iCov(z_i, z_o)} + Cov(z_o, z_o) \end{aligned} J=D(i=1∑nλizi–zo)=D(i=1∑nλizi)–2Cov(i=1∑nλizi,zo)+Cov(zo,zo)=i=1∑nj=0∑nλiλjCov(zi,zj)–2i=1∑nλiCov(zi,zo)+Cov(zo,zo)
注解:
D
(
x
)
=
E
(
(
x
−
E
(
x
)
)
2
)
=
E
(
x
2
−
2
x
E
(
x
)
−
(
E
(
x
)
)
2
)
=
E
(
x
2
)
−
2
E
(
x
)
E
(
x
)
+
(
E
(
x
)
)
2
\begin{aligned} D(x)&=E((x-E(x))^2)\\&=E(x^2-2xE(x)-(E(x))^2)\\&=E(x^2)-2E(x)E(x)+(E(x))^2\end{aligned}
D(x)=E((x−E(x))2)=E(x2−2xE(x)−(E(x))2)=E(x2)−2E(x)E(x)+(E(x))2
且有
E
(
∑
i
=
1
n
λ
i
z
i
−
z
o
)
=
0
E(\sum^{n}_{i=1}{\lambda_iz_i}- z_o)=0
E(i=1∑nλizi−zo)=0
于是
J
=
E
(
(
∑
i
=
1
n
λ
i
z
i
–
z
o
)
2
)
=
E
(
(
∑
i
=
1
n
λ
i
z
i
)
2
)
–
2
E
(
∑
i
=
1
n
λ
i
z
i
∗
z
o
)
+
E
(
z
o
2
)
=
∑
i
=
1
n
∑
j
=
0
n
λ
i
λ
j
E
(
z
i
∗
z
j
)
–
2
∑
i
=
1
n
λ
i
E
(
z
i
∗
z
o
)
+
D
(
z
o
)
=
∑
i
=
1
n
∑
j
=
0
n
λ
i
λ
j
C
o
v
(
z
i
,
z
j
)
–
2
∑
i
=
1
n
λ
i
C
o
v
(
z
i
,
z
o
)
+
C
o
v
(
z
o
,
z
o
)
\begin{aligned} J&= E((\sum^{n}_{i=1}{\lambda_iz_i} – z_o)^2) \\&= E((\sum^{n}_{i=1}{\lambda_iz_i})^2) – 2 E(\sum^{n}_{i=1}{\lambda_iz_i}* z_o) + E(z_o^2) \\&= \sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_jE( z_i *z_j)} – 2 \sum^{n}_{i=1}{\lambda_iE(z_i*z_o)} +D(z_o) \\&= \sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_jCov( z_i, z_j)} – 2 \sum^{n}_{i=1}{\lambda_iCov(z_i, z_o)} + Cov(z_o, z_o) \end{aligned}
J=E((i=1∑nλizi–zo)2)=E((i=1∑nλizi)2)–2E(i=1∑nλizi∗zo)+E(zo2)=i=1∑nj=0∑nλiλjE(zi∗zj)–2i=1∑nλiE(zi∗zo)+D(zo)=i=1∑nj=0∑nλiλjCov(zi,zj)–2i=1∑nλiCov(zi,zo)+Cov(zo,zo)
为简化描述,定义符号 C i j = C o v ( z i , z j ) = C o v ( R i , R j ) C_{ij}=Cov(z_i,z_j)=Cov(R_i,R_j) Cij=Cov(zi,zj)=Cov(Ri,Rj),这里 R i = z i – c R_i=z_i–c Ri=zi–c,即点 ( x i , y i ) (x_i,y_i) (xi,yi)处的属性值相对于区域平均属性值的偏差。 则有
J = ∑ i = 1 n ∑ j n λ i λ j C i j – 2 ∑ i = 1 n λ i C i o + C o o J = \sum^{n}_{i=1}\sum^{n}_{j}{\lambda_i\lambda_jC_{ij}} – 2 \sum^{n}_{i=1}{\lambda_iC_{io}} + C_{oo} J=i=1∑nj∑nλiλjCij–2i=1∑nλiCio+Coo
1.5代价函数的最优解
再定义半方差函数 r i j = σ 2 − C i j r_{ij}=σ^2−C_{ij} rij=σ2−Cij,带入 J J J中,有
J
=
∑
i
=
1
n
∑
j
=
0
n
λ
i
λ
j
(
σ
2
–
r
i
j
)
–
2
∑
i
=
1
n
λ
i
(
σ
2
–
r
i
o
)
+
σ
2
–
r
o
o
=
∑
i
=
1
n
∑
j
=
0
n
λ
i
λ
j
(
σ
2
)
−
∑
i
=
1
n
∑
j
=
0
n
λ
i
λ
j
(
r
i
j
)
−
2
∑
i
=
1
n
λ
i
(
σ
2
)
+
2
∑
i
=
1
n
λ
i
(
r
i
o
)
+
σ
2
–
r
o
o
\begin{aligned}J & = \sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_j(\sigma^2 – r_{ij})} – 2 \sum^{n}_{i=1}{\lambda_i(\sigma^2 – r_{io})} + \sigma^2 – r_{oo} \\ &=\sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_j(\sigma^2)} -\sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_j( r_{ij})}-2\sum^{n}_{i=1}{\lambda_i(\sigma^2)}+2 \sum^{n}_{i=1}{\lambda_i(r_{io})}+\sigma^2 – r_{oo} \end{aligned}
J=i=1∑nj=0∑nλiλj(σ2–rij)–2i=1∑nλi(σ2–rio)+σ2–roo=i=1∑nj=0∑nλiλj(σ2)−i=1∑nj=0∑nλiλj(rij)−2i=1∑nλi(σ2)+2i=1∑nλi(rio)+σ2–roo
考虑到
∑
i
=
1
n
λ
i
=
1
\sum^{n}_{i=1}{\lambda_i} = 1
∑i=1nλi=1
J
=
σ
2
−
∑
i
=
1
n
∑
j
n
λ
i
λ
j
(
r
i
j
)
−
2
σ
2
+
2
∑
i
=
1
n
λ
i
(
r
i
o
)
+
σ
2
–
r
o
o
=
2
∑
i
=
1
n
λ
i
(
r
i
o
)
−
∑
i
=
1
n
∑
j
=
0
n
λ
i
λ
j
(
r
i
j
)
–
r
o
o
\begin{aligned}J &= \sigma^2-\sum^{n}_{i=1}\sum^{n}_{j}{\lambda_i\lambda_j(r_{ij})}-2 \sigma^2 +2 \sum^{n}_{i=1}{\lambda_i(r_{io})}+ \sigma^2 – r_{oo}\\&=2 \sum^{n}_{i=1}{\lambda_i(r_{io})} -\sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_j(r_{ij})} – r_{oo} \end{aligned}
J=σ2−i=1∑nj∑nλiλj(rij)−2σ2+2i=1∑nλi(rio)+σ2–roo=2i=1∑nλi(rio)−i=1∑nj=0∑nλiλj(rij)–roo
我们的目标是寻找使J最小的一组
λ
i
λ_i
λi,且
J
J
J是
λ
i
λ_i
λi的函数,因此直接将
J
J
J对
λ
i
λ_i
λi求偏导数令其为
0
0
0即可。即
∂
J
∂
λ
i
=
0
;
i
=
1
,
2
,
⋯
,
n
\frac{\partial J}{\partial \lambda_i}= 0;i=1,2,\cdots,n
∂λi∂J=0;i=1,2,⋯,n
但是要注意的是,我们要保证求解出来的最优
λ
i
λ_i
λi 满足公式
∑
i
=
1
n
λ
i
=
1
\sum^{n}_{i=1}{\lambda_i} = 1
∑i=1nλi=1,这是一个带约束条件的最优化问题。使用拉格朗日乘数法求解,求解方法为构造一个新的目标函数
J
+
2
ϕ
(
∑
i
=
1
n
λ
i
−
1
)
J + 2\phi(\sum^{n}_{i=1}{\lambda_i}-1)
J+2ϕ(i=1∑nλi−1)
其中
ϕ
ϕ
ϕ是拉格朗日乘数。求解使这个代价函数最小的参数集
ϕ
,
λ
1
,
λ
2
,
⋯
,
λ
n
ϕ,λ_1,λ_2,⋯,λ_n
ϕ,λ1,λ2,⋯,λn,则能满足其在
∑
i
=
1
n
λ
i
=
1
\sum^{n}_{i=1}{\lambda_i} = 1
∑i=1nλi=1约束下最小化
J
J
J。即
{
∂
(
J
+
2
ϕ
(
∑
i
=
1
n
λ
i
−
1
)
)
∂
λ
k
=
0
;
k
=
1
,
2
,
⋯
,
n
∂
(
J
+
2
ϕ
(
∑
i
=
1
n
λ
i
−
1
)
)
∂
ϕ
=
0
\left\{\begin{aligned}\frac{\partial(J + 2\phi(\sum^{n}_{i=1}{\lambda_i}-1))}{\partial \lambda_k} &= 0;k=1,2,\cdots,n\\ \frac{\partial(J + 2\phi(\sum^{n}_{i=1}{\lambda_i}-1))}{\partial \phi} &= 0 \end{aligned} \right.
⎩⎪⎪⎪⎨⎪⎪⎪⎧∂λk∂(J+2ϕ(∑i=1nλi−1))∂ϕ∂(J+2ϕ(∑i=1nλi−1))=0;k=1,2,⋯,n=0
{
∂
(
2
∑
i
=
1
n
λ
i
(
r
i
o
)
–
∑
i
=
1
n
∑
j
n
λ
i
λ
j
(
r
i
j
)
–
r
o
o
+
2
ϕ
(
∑
i
=
1
n
λ
i
−
1
)
)
∂
λ
k
=
0
;
k
=
1
,
2
,
⋯
,
n
∂
(
2
∑
i
=
1
n
λ
i
(
r
i
o
)
–
∑
i
=
1
n
∑
j
n
λ
i
λ
j
(
r
i
j
)
–
r
o
o
+
2
ϕ
(
∑
i
=
1
n
λ
i
−
1
)
)
∂
ϕ
=
0
\left\{\begin{aligned} \frac{\partial (2 \sum^{n}_{i=1}{\lambda_i(r_{io})} – \sum^{n}_{i=1}\sum^{n}_{j}{\lambda_i\lambda_j(r_{ij})} – r_{oo}+ 2\phi(\sum^{n}_{i=1}{\lambda_i}-1))}{\partial \lambda_k} & = 0;k=1,2,\cdots,n\\ \frac{\partial ( 2 \sum^{n}_{i=1}{\lambda_i(r_{io})} – \sum^{n}_{i=1}\sum^{n}_{j}{\lambda_i\lambda_j(r_{ij})} – r_{oo}+ 2\phi(\sum^{n}_{i=1}{\lambda_i}-1))}{\partial \phi} &= 0 \end{aligned} \right.
⎩⎪⎪⎪⎨⎪⎪⎪⎧∂λk∂(2∑i=1nλi(rio)–∑i=1n∑jnλiλj(rij)–roo+2ϕ(∑i=1nλi−1))∂ϕ∂(2∑i=1nλi(rio)–∑i=1n∑jnλiλj(rij)–roo+2ϕ(∑i=1nλi−1))=0;k=1,2,⋯,n=0
{
2
r
k
o
–
∑
j
=
1
n
(
r
k
j
+
r
j
k
)
λ
j
+
2
ϕ
=
0
;
k
=
1
,
2
,
⋯
,
n
∑
i
=
1
n
λ
i
=
1
\left\{\begin{aligned} 2r_{ko} – \sum^{n}_{j=1}{(r_{kj}+r_{jk})\lambda_j}+2\phi&=0;k=1,2,\cdots,n\\ \sum^{n}_{i=1}{\lambda_i} &= 1 \end{aligned} \right.
⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧2rko–j=1∑n(rkj+rjk)λj+2ϕi=1∑nλi=0;k=1,2,⋯,n=1
由于
C
i
j
=
C
o
v
(
z
i
,
z
j
)
=
C
j
i
C_{ij}=Cov(z_i,z_j)=C_{ji}
Cij=Cov(zi,zj)=Cji,因此同样地
r
i
j
=
r
j
i
r_{ij}=r_{ji}
rij=rji,那么有
{
r
k
o
–
∑
j
=
1
n
r
k
j
λ
j
+
ϕ
=
0
;
k
=
1
,
2
,
⋯
,
n
∑
i
=
1
n
λ
i
=
1
\left\{\begin{aligned} r_{ko} – \sum^{n}_{j=1}{r_{kj}\lambda_j}+\phi&= 0;k=1,2,\cdots,n\\ \sum^{n}_{i=1}{\lambda_i} &= 1 \end{aligned} \right.
⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧rko–j=1∑nrkjλj+ϕi=1∑nλi=0;k=1,2,⋯,n=1
式子中半方差函数
r
i
j
r_{ij}
rij十分重要,最后会详细解释其计算与定义 在以上计算中我们得到了对于求解权重系数
λ
j
λ_j
λj的方程组。写成线性方程组的形式就是:
{
r
11
λ
1
+
r
12
λ
2
+
⋯
+
r
1
n
λ
n
–
ϕ
=
r
1
o
r
21
λ
1
+
r
22
λ
2
+
⋯
+
r
2
n
λ
n
–
ϕ
=
r
2
o
⋯
r
n
1
λ
1
+
r
n
2
λ
2
+
⋯
+
r
n
n
λ
n
–
ϕ
=
r
n
o
λ
1
+
λ
2
+
⋯
+
λ
n
=
1
{\left\{\begin{aligned} r_{11}\lambda_1 + r_{12}\lambda_2 + \cdots + r_{1n}\lambda_n – \phi&= r_{1o}\\r_{21}\lambda_1 + r_{22}\lambda_2 + \cdots + r_{2n}\lambda_n – \phi&= r_{2o}\\&\cdots\\ r_{n1}\lambda_1 + r_{n2}\lambda_2 + \cdots + r_{nn}\lambda_n – \phi&= r_{no}\\ \lambda_1 + \lambda_2 + \cdots + \lambda_n &= 1\\ \end{aligned} \right.}
⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧r11λ1+r12λ2+⋯+r1nλn–ϕr21λ1+r22λ2+⋯+r2nλn–ϕrn1λ1+rn2λ2+⋯+rnnλn–ϕλ1+λ2+⋯+λn=r1o=r2o⋯=rno=1
写成矩阵形式即为
[
r
11
r
12
⋯
r
1
n
1
r
21
r
22
⋯
r
2
n
1
⋯
⋯
⋯
⋯
⋯
r
n
1
r
n
2
⋯
r
n
n
1
1
1
⋯
1
0
]
[
λ
1
λ
2
⋯
λ
n
−
ϕ
]
=
[
r
1
o
r
2
o
⋯
r
n
o
1
]
\begin{bmatrix}r_{11}&r_{12}&\cdots&r_{1n}&1\\ r_{21}&r_{22}&\cdots&r_{2n}&1\\\cdots&\cdots&\cdots&\cdots&\cdots\\r_{n1}&r_{n2}&\cdots&r_{nn}&1\\1&1&\cdots&1&0\end{bmatrix}\begin{bmatrix} \lambda_1\\ \lambda_2\\\cdots\\\lambda_n\\-\phi\end{bmatrix}=\begin{bmatrix} r_{1o}\\ r_{2o}\\\cdots\\r_{no}\\1\end{bmatrix}
⎣⎢⎢⎢⎢⎡r11r21⋯rn11r12r22⋯rn21⋯⋯⋯⋯⋯r1nr2n⋯rnn111⋯10⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡λ1λ2⋯λn−ϕ⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡r1or2o⋯rno1⎦⎥⎥⎥⎥⎤
对矩阵求逆即可求解。 唯一未知的就是上文中定义的半方差函数
r
i
j
r_{ij}
rij,接下来将详细讨论
注解:关于求解矩阵 A λ = b A\lambda=b Aλ=b,使用cholesky分解比直接求逆快
1.6半方差函数
上文中对半方差函数的定义为
r
i
j
=
σ
2
−
C
i
j
r_{ij} = \sigma^2 -C_{ij}
rij=σ2−Cij
其等价形式为
r
i
j
=
1
2
E
[
(
z
i
−
z
j
)
2
]
r_{ij} = \frac{1}{2}E[(z_i-z_j)^2]
rij=21E[(zi−zj)2]
这也是半方差函数名称的来由,接下来证明这二者是等价的:
根据上文定义
R
i
=
z
i
–
c
R_i=z_i–c
Ri=zi–c,有
z
i
−
z
j
=
R
i
–
R
j
z_i−z_j=R_i–R_j
zi−zj=Ri–Rj,则
r
i
j
=
1
2
E
[
(
R
i
−
R
j
)
2
]
=
1
2
E
[
R
i
2
−
2
R
i
R
j
+
R
j
2
]
=
1
2
E
[
R
i
2
]
+
1
2
E
[
R
j
2
]
−
E
[
R
i
R
j
]
\begin{aligned} r_{ij} &= \frac{1}{2}E[(R_i-R_j)^2]\\&= \frac{1}{2}E[R_i^2-2R_iR_j+R_j^2]\\&= \frac{1}{2}E[R_i^2]+\frac{1}{2}E[R_j^2]-E[R_iR_j] \end{aligned}
rij=21E[(Ri−Rj)2]=21E[Ri2−2RiRj+Rj2]=21E[Ri2]+21E[Rj2]−E[RiRj]
又因为:
E
[
R
i
2
]
=
E
[
R
j
2
]
=
E
[
(
z
i
–
c
)
2
]
=
V
a
r
(
z
i
)
=
σ
2
E[R_i^2] =E[R_j^2] = E[(z_i – c)^2] = Var(z_i) = \sigma^2
E[Ri2]=E[Rj2]=E[(zi–c)2]=Var(zi)=σ2
E
[
R
i
R
j
]
=
E
[
(
z
i
–
c
)
(
z
j
−
c
)
]
=
C
o
v
(
z
i
,
z
j
)
=
C
i
j
E[R_iR_j] = E[(z_i – c)(z_j-c)] = Cov(z_i,z_j) = C_{ij}
E[RiRj]=E[(zi–c)(zj−c)]=Cov(zi,zj)=Cij
于是有
r
i
j
=
1
2
E
[
(
z
i
−
z
j
)
2
]
=
1
2
E
[
R
i
2
]
+
1
2
E
[
R
j
2
]
−
E
[
R
i
R
j
]
=
1
2
σ
2
+
1
2
σ
2
−
C
i
j
=
σ
2
−
C
i
j
\begin{aligned}r_{ij} &= \frac{1}{2}E[(z_i-z_j)^2]\\&= \frac{1}{2}E[R_i^2]+\frac{1}{2}E[R_j^2]-E[R_iR_j]\\&= \frac{1}{2}\sigma^2+\frac{1}{2}\sigma^2- C_{ij}\\&=\sigma^2 -C_{ij}\end{aligned}
rij=21E[(zi−zj)2]=21E[Ri2]+21E[Rj2]−E[RiRj]=21σ2+21σ2−Cij=σ2−Cij
σ
2
−
C
i
j
=
1
2
E
[
(
z
i
−
z
j
)
2
]
\sigma^2 -C_{ij} = \frac{1}{2}E[(z_i-z_j)^2]
σ2−Cij=21E[(zi−zj)2]得证,现在的问题就是如何计算
r
i
j
=
1
2
E
[
(
z
i
−
z
j
)
2
]
r_{ij} = \frac{1}{2}E[(z_i-z_j)^2]
rij=21E[(zi−zj)2]
这时需要用到地理学第一定律,空间上相近的属性相近。
r
i
j
=
1
2
(
z
i
−
z
j
)
2
r_{ij} = \frac{1}{2}(z_i-z_j)^2
rij=21(zi−zj)2表达了属性的相似度;空间的相似度就用距离来表达,定义i与j之间的几何距离
克里金插值假设
r
i
j
r_{ij}
rij与
d
i
j
d_{ij}
dij存在着函数关系,这种函数关系可以是线性、二次函数、指数、对数关系。为了确认这种关系,我们需要首先对观测数据集
{
z
(
x
1
,
y
1
)
,
z
(
x
2
,
y
2
)
,
z
(
x
3
,
y
3
)
,
⋯
,
z
(
x
n
−
1
,
y
n
−
1
)
,
z
(
x
n
,
y
n
)
}
\{z(x_1,y_1),z(x_2,y_2),z(x_3,y_3),\cdots,z(x_{n-1},y_{n-1}),z(x_n,y_n)\}
{z(x1,y1),z(x2,y2),z(x3,y3),⋯,z(xn−1,yn−1),z(xn,yn)}
计算任意两个点的 距离
d
i
j
=
(
x
i
−
x
j
)
2
+
(
y
i
–
y
j
)
2
d_{ij}= \sqrt{(x_i-x_j)^2 + (y_i – y_j)^2}
dij=(xi−xj)2+(yi–yj)2和 半方差
σ
2
−
C
i
j
=
1
2
E
[
(
z
i
−
z
j
)
2
]
\sigma^2 -C_{ij} =\frac{1}{2}E[(z_i-z_j)^2]
σ2−Cij=21E[(zi−zj)2],这时会得到
n
2
n^2
n2个
(
d
i
j
,
r
i
j
)
(d_{ij},r_{ij})
(dij,rij)的数据对。
将所有的
d
d
d和
r
r
r绘制成散点图,寻找一个最优的拟合曲线拟合
d
d
d与
r
r
r的关系,得到函数关系式
r
=
r
(
d
)
r = r(d)
r=r(d)
那么对于任意两点
(
x
i
,
y
i
)
,
(
x
j
,
y
j
)
(x_i,y_i), (x_j,y_j)
(xi,yi),(xj,yj),先计算其距离
d
i
j
d_{ij}
dij,然后根据得到的函数关系就可以得到这两点的半方差.
1.7半方差模型
线性有基台模型:
球状模型:
指数模型:
高斯模型:
不同模型之间可以线性叠加
1.8小结
总的来说,进行克里金插值分为这几个步骤:
- 对于观测数据,两两计算距离与半方差
- 寻找一个拟合曲线拟合距离与半方差的关系,从而能根据任意距离计算出相应的半方差
- 计算出所有已知点之间的半方差 r i j r_{ij} rij
- 对于未知点 z o z_o zo,计算它到所有已知点 z i z_i zi的半方差 r i o r_{io} rio
- 求解第四节中的方程组,得到最优系数 λ i λ_i λi
- 使用最优系数对已知点的属性值进行加权求和,得到未知点 z o z_o zo的估计值
2.Py实现
GeostatsGuy Lectures: https://github.com/GeostatsGuy/PythonNumericalDemos/blob/master/GeostatsPy_kriging.ipynb.
2.1 小结
待补充。。。