簡介：
Kriging是一組統計技術，用來將隨機場的值(例如，地形的高程，z，作為地理位置的函數)從其在附近位置的觀測值中插值到一個未觀測到的位置。
令\((x,y)\)為某一空間數據點的坐標，\(Z(x,y)\)表示其值，對於某個點\((x_0,y_0)\)的值記為\(z_0=Z(x_0, y_0)\)，則克里金插值公式如下：
\(\hat{z}_0=\sum_{i=1}^n\omega_iz_i\)
\(\hat{z}_0\)表示點\((x_0,y_0)\)處的估計值，它由該點附近n個點的觀測值\(z_i\)(1<=i<=n)與相對該點\((x_0,y_0)\)的權重系數\(\omega_i\)的線性組合決定。
一、數學原理
1、假設條件
（1）數學期望未知但是一個恆定的值\(E[Z(x,y)]={\mu}\);
（2）\(Z(x,y)\)的變異函數是已知的\({\gamma}(x, y)_{ij} = E[(z_i - z_j)^2]\).
解釋：普通克里金（Ordinary Kriging)的假設條件為：空間屬性z是均一的，對空間任意一點（x,y)都有相同的期望
因此可以據此得出任意一點\((x,y)\)的都由區域均值\(\mu\)和該點的隨機偏差\(R(x,y)\)組成，即\(Z(x,y)=E[Z(x,y)]+R(x,y)]=\mu+R(x,y)\)
\(R(x,y)=Z(x,y)-\mu\)，而方差\(\sigma^2=Var(Z(x,y))=E[(Z(x,y)-\mu)^2]=E[R(x,y)^2]\)
所以隨機偏差的方差也是常數\(\sigma^2\)，\(Var[R(x,y)]=\sigma^2\)
下文中絕大多數為方便起見都將\(Z(x,y)\)記作\(z_i\)，\(R(x,y)\)記作\(r_i\)
2、無偏估計條件
普通克里金的權重\(\omega\)需要滿足無偏估計，即
\(\sum_{i=1}^{n}{\omega_i}=1\)
證明如下：
無偏估計條件是 \(E(\hat{z_{0}}-z_{0})=0\), 代入\(\hat{z_{0}}=\sum_{i=1}^{n}{\omega_iz_i}\)得：
\(E(\sum_{i=1}^{n}{\omega_iz_i}-z_0)=0\)
因為對任意的z都存在E[z] = μ, 則有
\(μ\sum_{i=1}^{n}{\omega_i}-μ=0\)
約去μ得證
\(\sum_{i=1}^{n}{\omega_i}=1\)
3、半方差
先復習下數學期望、方差、協方差公式之間的關系：
設X、Y是隨機變量，C、a、b均為常數,記\(D(X)\)和\(D(Y)\)分別表示各自的方差，\(E(X)\)和\(E(Y)\)分別表示各自的數學期望，\(Cov(X,Y)\)表示X與Y的協方差，則有如下關系：
\(E(C)=C\)
\(E(CX)=CE(X)\)
\(E(CX)=CE(X)\)
\(E(X+Y)=E(X)+E(Y)\)
當X和Y相互獨立時，\(E(XY)=E(X)E(Y)\)
\(D(X+Y)=D(X)+D(Y)+2Cov(X,Y)\)
\(D(X-Y)=D(X)+D(Y)-2Cov(X,Y)\)
\(Cov(X,Y)=E(XY)-E(X)E(Y)\)
\(Cov(X,Y)=Cov(Y,X)\)
\(Cov(aX,bY)=ab{\cdot}Cov(X,Y)\)
\(Cov(X_1+X_2,Y)=Cov(X_1,Y)+Cov(X_2,Y)\)
\(Cov(X+a,Y+b)=Cov(X,Y)\)
\(Cov(X,X)=D(X)\)
\(Cov(Y,Y)=D(Y)\)
半方差定義：\(r_{ij}=\sigma^2-Cov(z_i,z_j)\)，\(\sigma^2\)表示方差；其等價形式為\(r_{ij}=\frac{1}{2}E[(z_i-z_j)^2]\)，證明如下：
由上文知\(r_i=z_i-\mu\)，則\(z_i-z_j=r_i-r_j\)，從而
\(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]-E[r_ir_j]+\frac{1}{2}E[r_j^2]\)
又\(E[r_i^2]=E[r_j^2]=E[(z_i-\mu)^2]=Var(z_i)=\sigma^2\)
\(E[r_ir_j]=E[(z_i-\mu)(z_j-\mu)]=Cov(z_i,z_j)\)
於是有
\(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-Cov(z_i,z_j)\)
\(=\sigma^2-Cov(z_i,z_j)\)
4、優化目標函數
克里金誤差（kriging variance or kriging error)給出如下：
\(\sigma_k^2(z_0):=Var(\hat{z_0}-z_0)\)
推導\(\sigma_k^2(z_0):=Var(\hat{z_0}-z_0)\)
=\(Var(\hat{z_0})+Var(z_0)-2{\cdot}Cov(\hat{z_0}-z_0)\)
又\(Var(\hat{z_0})=Var(\sum_{i=1}^n\omega_iz_i)\)
\(=Cov(\sum_{i=1}^n\omega_iz_i,\sum_{i=1}^n\omega_iz_i)\)
\(=Cov(\sum_{i=1}^n\omega_iz_i,\sum_{j=1}^n\omega_jz_j)\)(第二個參數中的i用j替換)
\(=\sum_{i=1}^n\omega_i\sum_{j=1}^n\omega_j{\cdot}Cov(z_i,z_j)\)
\(=\sum_{i=1}^n\sum_{j=1}^n\omega_i\omega_j{\cdot}Cov(z_i,z_j)\)
則有：\(\sigma_k^2(z_0):=Var(\hat{z_0}-z_0)\)
=\(\sum_{i=1}^n\sum_{j=1}^n\omega_i\omega_j{\cdot}Cov(z_i,z_j)+Cov(z_0,z_0)-2\sum_{i=1}^n\omega_i{\cdot}Cov(z_i,z_0)\)
由上文知:\(Cov(z_i,z_j)=\sigma^2-r_{ij}\)，\(\sum_{i=1}^n=1\)
從而\(\sigma_k^2(z_0):=\sum_{i=1}^n\sum_{j=1}^n\omega_i\omega_j(\sigma^2-r_{ij})+(\sigma^2-r_{00})-2\sum_{i=1}^nw_i(\sigma^2-r_{i0})\)
=\(\sum_{i=1}^n\sum_{j=1}^nw_iw_j(\sigma^2)-\sum_{i=1}^n\sum_{j=1}^nw_iw_j(r_{ij})+(\sigma^2-r_{00})-2\sum_{i=1}^nw_i(\sigma^2)+2\sum_{i=1}^nw_i(r_{i0})\)
=\(\sigma^2-\sum_{i=1}^n\sum_{j=1}^nw_iw_j(r_{ij})+\sigma^2-r_{00}-2\sigma^2+2\sum_{i=1}^nw_i(r_{i0})\)
=\(2\sum_{i=1}^nw_i(r_{i0})-\sum_{i=1}^n\sum_{j=1}^nw_iw_j(r_{ij})-r_{00}\)
接下來尋找使\(\sigma_k^2(z_0)\)最小的一組\(w_i\)，應用拉格朗日乘數法求解。
回顧拉格朗日乘數法知識：
設給定二元函數z=ƒ(x,y)和附加條件φ(x,y)=0，為尋找z=ƒ(x,y)在附加條件下的極值點，先做拉格朗日函數 ，其中λ為參數。
令F(x,y,λ)對x和y和λ的一階偏導數等於零，即
F'x=ƒ'x(x,y)+λφ'x(x,y)=0 [1]
F'y=ƒ'y(x,y)+λφ'y(x,y)=0
F'λ=φ(x,y)=0
由上述方程組解出x,y及λ，如此求得的(x,y)，就是函數z=ƒ(x,y)在附加條件φ(x,y)=0下的可能極值點。
若這樣的點只有一個，由實際問題可直接確定此即所求的點。
為使\(\sigma_k^2(z_0)\)最小，構造拉格朗日函數\(J=F(w_0,w_1,...,w_n,\lambda)=f(w_0,w_1,...,w_n)+\lambda\phi(\sum_{i=1}^n-1)=\sigma_k^2(z_0)+\lambda\phi(\sum_{i=1}^n-1)\)
其中\(\lambda\phi(\sum_{i=1}^n-1)\)是約束條件，\(\lambda\)是拉格朗日乘數。
將\(J\)對\(\omega_i\)求偏導並令其值為0，即\(\frac{\partial J}{\partial \omega_i}= 0;i=1,2,\cdots,n\)
求解過程暫略，最終會得到如下的矩陣方程：
\(\begin{bmatrix}w_1\\ w_2\\ ...\\w_n\\ k\lambda\\\end{bmatrix}=\begin{bmatrix} r_{11}&r_{12}&...&r_{1n}&1\\ r_{21}&r_{22}&...&r_{2n}&1\\ ...&...&...&...&...\\ r_{n1}&r_{n2}&...&r_{nn}&1\\ 1& 1& ...& 1 & 0 \end{bmatrix}^T\begin{bmatrix} r_{1o}\\ r_{2o}\\ ...\\ r_{no}\\ 1\\ \end{bmatrix}\)
向量\(W=\begin{bmatrix}w_1\\ w_2\\ ...\\w_n\\ k\lambda\\\end{bmatrix}\)中的\(w\)就是用於最后計算預測值的權重。
5、半變異函數模型擬合
對於4中的求解需要先擬合半方差函數以獲得半變異模型，然后再根據半變異模型獲取滿足空間自相關的一定范圍（range)內的觀察值，運用這些觀測值結合4來計算最終的預測值。
半方差的計算方式如下：計算預測點附近N個觀測值兩兩之間的空間位置歐幾里德距離h和值誤差（均方誤差或半方差等）r，然后繪制r關於h的散點圖，一般情況下，如果地理屬性之間存在空間自相關，那么
r會隨着h的增大逐漸增大，h大到一定程度（range)后r會收斂於某個值（sill)，表示超過h這個距離后其他的觀測值對預測值的自相關性已經可以忽略不計了。此時就應該用range范圍內的觀測數據去計算預測值。
幾種常見的函數模型如下：
線性模型用的應該比較少，線性模型可以用線性最小二乘法擬合以獲得函數參數，對於球體、指數、高斯等其他曲線模型可以用非線性最小二乘法擬合以獲得函數參數。
關於擬合函數模型這塊暫略，以后再補充筆記。
個人推測實際應用克里金插值的散點數量如果比較大應該不是把所有的點都擬合一次來進行預測，觀察ArcMap的克里金插值是設置的橢圓半徑范圍內的最大最小相鄰元素數目，所以應該就是用預測值附近選取的N個點擬合半變異模型，獲得nugget、range、sill參數，然后用range范圍內的點估計預測值。如下圖所示操作：
我用線性最小二乘擬合模型計算存在些問題，就不貼代碼了。
注釋：個人水平所限，筆記中難免存在些錯誤，如果讀者發現了可留言指出，一起學習成長進步哦:)。
References：
http://wiki.gis.com/wiki/index.php/Kriging#Ordinary_kriging_equation
https://desktop.arcgis.com/zh-cn/arcmap/10.3/tools/3d-analyst-toolbox/how-kriging-works.htm
https://xg1990.com/blog/en/archives/222