學過空間插值的人都知道克里金插值，但是它的變種繁多、公式復雜，還有個半方差函數讓人不知所云

本文講簡單介紹基本克里金插值的原理，及其推理過程，全文分為九個部分：

0.引言－從反距離插值說起

1.克里金插值的定義

2.假設條件

3.無偏約束條件

4.優化目標／代價函數

5.代價函數的最優解

6.半方差函數

7.普通克里金與簡單克里金

8.小結

0.引言——從反距離插值(IDW)說起

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

直觀來講，根據地理學第一定律，

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

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

\[\hat{z} = \sum^{n}_{i=1}{\frac{1}{d^\alpha}z_i}\]

其中\(\alpha\)通常取1或者2。

即，用空間上所有已知點的數據加權求和來估計未知點的值，權重取決于距離的倒數(或者倒數的平方)。那么，距離近的點，權重就大；距離遠的點，權重就小。

反距離插值可以有效的基于地理學第一定律估計屬性值空間分布，但仍然存在很多問題：

\(\alpha\)的值不確定

用倒數函數來描述空間關聯程度不夠準確

因此更加準確的克里金插值方法被提出來了

1.克里金插值的定義

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

\[\hat{z_o} = \sum^{n}_{i=1}{\lambda_iz_i}\]

其中\(\hat{z_o}\)是點\((x_o,y_o)\)處的估計值，即\(z_o=z(x_o,y_o)\) 。

這里的\(\lambda_i\)是權重系數。它同樣是用空間上所有已知點的數據加權求和來估計未知點的值。但權重系數并非距離的倒數，而是能夠滿足點\((x_o,y_o)\)處的估計值\(\hat{z_o}\)與真實值\(z_o\)的差最小的一套最優系數，即

\[\min_{\lambda_i} Var(\hat{z_o}-z_o)\]

同時滿足無偏估計的條件

\[E(\hat{z_o}-z_o)=0\]

2.假設條件

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

普通克里金插值的假設條件為，空間屬性\(z\)是均一的。對于空間任意一點\((x,y)\)，都有同樣的期望c與方差\(\sigma^2\)。

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

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

\[Var[z(x,y)] = \sigma^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)\)處的偏差，其方差均為常數

\[Var[R(x,y)] = \sigma^2\]

3.無偏約束條件

先分析無偏估計條件\(E(\hat{z_o}-z_o)=0\)，將\(\hat{z_o} = \sum^{n}_{i=1}{\lambda_iz_i}\)帶入則有

\[E(\sum^{n}_{i=1}{\lambda_iz_i}- z_o)=0\]

又因為對任意的z都有\(E[z] = c\)，則

\[c \sum^{n}_{i=1}{\lambda_i}- c=0\]

即

\[\sum^{n}_{i=1}{\lambda_i} = 1\]

這是\(\lambda_i\)的約束條件之一。

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

再分析估計誤差\(Var(\hat{z_o}-z_o)\)。為方便公式推理，用符號\(J\)表示，即

\[J = Var(\hat{z_o}-z_o)\]

則有

\[\begin{array}{r@{\;=\;}l} J &= Var(\sum^{n}_{i=1}{\lambda_iz_i} – z_o) \\&= Var(\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{array} \]

為簡化描述，定義符號 \(C_{ij} = Cov(z_i,z_j) = Cov(R_i,R_j)\)，這里\(R_i = z_i – c\)，即點\((x_i,y_i)\)處的屬性值相對于區域平均屬性值的偏差。

則有

\[J = \sum^{n}_{i=1}\sum^{n}_{j}{\lambda_i\lambda_jC_{ij}} – 2 \sum^{n}_{i=1}{\lambda_iC_{io}} + C_{oo} \]

5.代價函數的最優解

再定義半方差函數 \(r_{ij} = \sigma^2 -C_{ij}\)，帶入J中，有

\[\begin{array}{r@{\;=\;}l}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{array} \]

考慮到\(\sum^{n}_{i=1}{\lambda_i} = 1\)

\[\begin{array}{r@{\;=\;}l}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{array} \]

我們的目標是尋找使J最小的一組 \(\lambda_i\)，且J是\(\lambda_i\)的函數，因此直接將J對\(\lambda_i\)求偏導數令其為0即可。即

\[\frac{\partial J}{\partial \lambda_i}= 0;i=1,2,\cdots,n\]

但是要注意的是，我們要保證求解出來的最優 \(\lambda_i\) 滿足公式\(\sum^{n}_{i=1}{\lambda_i} = 1\)，這是一個帶約束條件的最優化問題。使用拉格朗日乘數法求解，求解方法為構造一個新的目標函數

\[J + 2\phi(\sum^{n}_{i=1}{\lambda_i}-1)\]

其中\(\phi\)是拉格朗日乘數。求解使這個代價函數最小的參數集\({\phi,\lambda_1,\lambda_2,\cdots,\lambda_n}\)，則能滿足其在\(\sum^{n}_{i=1}{\lambda_i} = 1\)約束下最小化\(J\)。即

\[\left\{\begin{array}{r@{\;=\;}l}\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{array} \right.\]

\[\left\{\begin{array}{r@{\;=\;}l} \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{array} \right.\]

\[\left\{\begin{array}{r@{\;=\;}l} 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{array} \right.\]

由于\(C_{ij}=Cov(z_i,z_j)=C_{ji}\)，因此同樣地\(r_{ij}=r_{ji}\)，那么有

\[\left\{\begin{array}{r@{\;=\;}l} 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{array} \right.\]

式子中半方差函數\(r_{ij}\)十分重要，最后會詳細解釋其計算與定義

在以上計算中我們得到了對于求解權重系數\(\lambda_j\)的方程組。寫成線性方程組的形式就是：

\begin{equation}\left\{\begin{array}{r@{\;=\;}l} 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{array} \right.\end{equation}

寫成矩陣形式即為

\[\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}\]

對矩陣求逆即可求解。

唯一未知的就是上文中定義的半方差函數\(r_{ij}\)，接下來將詳細討論

6.半方差函數

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

\[r_{ij} = \sigma^2 -C_{ij}\]

其等價形式為

\[r_{ij} = \frac{1}{2}E[(z_i-z_j)^2]\]

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

根據上文定義 \(R_i = z_i – c\)，有\(z_i-z_j = R_i – R_j\)，則

\[\begin{array}{r@{\;=\;}l} 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{array} \]

又因為：

\[E[R_i^2] =E[R_j^2] = E[(z_i – c)^2] = Var(z_i) = \sigma^2 \]

\[E[R_iR_j] = E[(z_i – c)(z_j-c)] = Cov(z_i,z_j) = C_{ij}\]

于是有

\[\begin{array}{r@{\;=\;}l} 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{array}\]

\( \sigma^2 -C_{ij} = \frac{1}{2}E[(z_i-z_j)^2]\)得證，現在的問題就是如何計算

\[r_{ij} = \frac{1}{2}E[(z_i-z_j)^2]\]

這時需要用到地理學第一定律，空間上相近的屬性相近。\(r_{ij} = \frac{1}{2}(z_i-z_j)^2\)表達了屬性的相似度；空間的相似度就用距離來表達，定義i與j之間的幾何距離

\[d_{ij} = d(z_i,z_j) = d( (x_i,y_i), (x_j,y_j)) = \sqrt{(x_i-x_j)^2 + (y_i – y_j)^2}\]

克里金插值假設\(r_{ij}\)與\(d_{ij}\)存在著函數關系，這種函數關系可以是線性、二次函數、指數、對數關系。為了確認這種關系，我們需要首先對觀測數據集

\[\{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)\}\]

計算任意兩個點的 距離\(d_{ij}= \sqrt{(x_i-x_j)^2 + (y_i – y_j)^2}\)和 半方差 \(\sigma^2 -C_{ij} =\frac{1}{2}E[(z_i-z_j)^2]\)，這時會得到\(n^2\)個\((d_{ij}, r_{ij})\)的數據對。

將所有的\(d\)和\(r\)繪制成散點圖，尋找一個最優的擬合曲線擬合\(d\)與\(r\)的關系，得到函數關系式

\[r = r(d)\]

那么對于任意兩點\((x_i,y_i), (x_j,y_j)\)，先計算其距離\(d_{ij}\)，然后根據得到的函數關系就可以得到這兩點的半方差\(r_{ij}\)

7. 簡單克里金(simple kriging)與普通克里金(ordinary kriging)的區別

以上介紹的均為普通克里金(ordinary kriging)的公式與推理。

事實上普通克里金插值還有簡化版，即簡單克里金(simple kriging)插值。二者的差異就在于如何定義插值形式：

上文講到，普通克里金插值形式為

\[\hat{z_o} = \sum^{n}_{i=1}{\lambda_iz_i}\]

而簡單克里金的形式則為

\[\hat{z_o} – c= \sum^{n}_{i=1}{\lambda_i(z_i-c)}\]

這里的符號\(c\)在上文介紹過了，是屬性值的數學期望，即\(E[z] = c\)。也就是說，在普通克里金插值中，認為未知點的屬性值是已知點的屬性值的加權求和；而在簡單克里金插值中，假設未知點的屬性值相對于平均值的偏差是已知點的屬性值相對于平均值的偏差的加權求和，用公式表達即為：

\[\hat{R_o} = \sum^{n}_{i=1}{\lambda_iR_i}\]

這里的\(R_i\)在上文定義過了：\(R_i = z_i – c\)。

但是為什么這樣的克里金插值稱為“簡單克里金”呢？由于有假設\(E[z] = c\)，也就是說\(E(R_i + c) = c\)，即\(E(R_i) = 0\)。那么上面的公式\(\hat{R_o} = \sum^{n}_{i=1}{\lambda_iR_i}\)兩邊的期望一定相同，那么在求解未知參數\(\lambda_i\)就不需要有無偏約束條件\(\sum^{n}_{i=1}{\lambda_i} = 1\)。換句話說，這樣的估計公式天生就能滿足無偏條件。因此它被稱為簡單克里金。

從在上文(第4節優化目標/代價函數J)中可以知道，優化目標的推理和求解過程是通過對屬性值相對于期望的偏差量\(R_i\)進行數學計算而進行的。也就是說這兩種克里金插值方法雖然插值形式不一樣，求解方法是一樣的，重要的區別是簡單克里金插值不需要約束條件\(\sum^{n}_{i=1}{\lambda_i} = 1\)，求解方程組為：

\begin{equation}\left\{\begin{array}{r@{\;=\;}l} 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}\\ \end{array} \right.\end{equation}

還有更重要的一點，簡單克里金的插值公式為：

\[\hat{z_o} = \sum^{n}_{i=1}{\lambda_i(z_i-c)}+c\]

換句話說，在計算未知點屬性值\(\hat{z_o}\)前，需要知道該地區的屬性值期望\(c\)。事實上我們在進行插值前很難知道這個地區的真實屬性值期望。有些研究者可能會采用對觀測數據簡單求平均的方法計算期望值\(c\)，而考慮到空間采樣點位置代表性可能有偏差(比如采樣點聚集在某一小片地區，沒有代表性)，簡單平均估計的期望也可能是有偏差的。這是簡單克里金方法的局限性。

8.小結

總的來說，進行克里金插值分為這幾個步驟：

對于觀測數據，兩兩計算距離與半方差

尋找一個擬合曲線擬合距離與半方差的關系，從而能根據任意距離計算出相應的半方差

計算出所有已知點之間的半方差\(r_{ij}\)

對于未知點\(z_o\)，計算它到所有已知點\(z_i\)的半方差\(r_{io}\)

求解第四節中的方程組，得到最優系數\(\lambda_i\)

使用最優系數對已知點的屬性值進行加權求和，得到未知點\(z_o\)的估計值

總結

以上是生活随笔為你收集整理的matlab克里金插值法,克里金(Kriging)插值的原理与公式推导的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。