2021-08-20

前言

之前做毕设的时候看文献一知半解，现在找个时间补回来。

1.普通克里金推导

1.1定义

克里金插值公式形如：

$$\widehat{z_o}=\sum _{i=1}^n{\lambda }_i{z}_i$$
$\widehat{z_o}$ 为预测点估计值；
${\lambda }_i$ 为权重系数；
$z_i$ 已知采样点数据；

同时${\lambda }_i$ 为满足估计值 $\widehat{z_o}$与真实值$z_o$的差最小的一套最优系数$$\underset{{\lambda }_i}{\min }D(\widehat{z_o}-z_o)$$
$\widehat{z_o}$满足$$E(\widehat{z_o}-z_o)=0$$

1.2假设条件

普通克里金插值的假设条件为，空间属性$z$是均一的。对于空间任意一点$(x,y)$，都有同样的期望$c$与方差${\sigma }^2$。 即对任意点$(x,y)$都有
$$E[z(x,y)]=E[z]=c$$$$D[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)$处的偏差，其方差均为常数

$$D[R(x,y)]={\sigma }^2$$

1.3无偏约束条件

先分析无偏估计条件$E(\widehat{z_o}-z_o)=0$，将 $\widehat{z_o}={\sum }_{i=1}^n{\lambda }_i{z}_i$带入则有

$$E(\sum _{i=1}^n{\lambda }_i{z}_i-z_o)=0$$
又因为对任意的$z$都有$E[z]=c$，则
$$c\sum _{i=1}^n{\lambda }_i-c=0$$
即
$$\sum _{i=1}^n{\lambda }_i=1$$
这是 ${\lambda }_i$的约束条件之一。

1.4优化目标/代价函数J

再分析估计误差 $D(\widehat{z_o}-z_o)$。为方便公式推理，用符号$J$表示，即

$$J=D(\widehat{z_o}-z_o)$$
则有

$$\begin{array}{rl}J& =D(\sum _{i=1}^n{\lambda }_i{z}_i–z_o)\\ & =D(\sum _{i=1}^n{\lambda }_i{z}_i)–2Cov(\sum _{i=1}^n{\lambda }_i{z}_i,z_o)+Cov(z_o,z_o)\\ & =\sum _{i=1}^n\sum _{j=0}^n{\lambda }_i{\lambda }_{j}Cov(z_i,z_j)–2\sum _{i=1}^n{\lambda }_{i}Cov(z_i,z_o)+Cov(z_o,z_o)\end{array}$$

注解：
$$\begin{array}{rl}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{array}$$
且有
$$E(\sum _{i=1}^n{\lambda }_i{z}_i-z_o)=0$$
于是$$\begin{array}{rl}J& =E((\sum _{i=1}^n{\lambda }_i{z}_i–z_o{)}^2)\\ & =E((\sum _{i=1}^n{\lambda }_i{z}_i{)}^2)–2E(\sum _{i=1}^n{\lambda }_i{z}_i\ast z_o)+E(z_o^2)\\ & =\sum _{i=1}^n\sum _{j=0}^n{\lambda }_i{\lambda }_{j}E(z_i\ast z_j)–2\sum _{i=1}^n{\lambda }_{i}E(z_i\ast z_o)+D(z_o)\\ & =\sum _{i=1}^n\sum _{j=0}^n{\lambda }_i{\lambda }_{j}Cov(z_i,z_j)–2\sum _{i=1}^n{\lambda }_{i}Cov(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 _{i=1}^n\sum _j^n{\lambda }_i{\lambda }_j{C}_{ij}–2\sum _{i=1}^n{\lambda }_i{C}_{io}+C_{oo}$$

1.5代价函数的最优解

再定义半方差函数 $r_{ij}={\sigma }^2-C_{ij}$，带入$J$中，有

$$\begin{array}{rl}J& =\sum _{i=1}^n\sum _{j=0}^n{\lambda }_i{\lambda }_j({\sigma }^2–r_{ij})–2\sum _{i=1}^n{\lambda }_i({\sigma }^2–r_{io})+{\sigma }^2–r_{oo}\\ & =\sum _{i=1}^n\sum _{j=0}^n{\lambda }_i{\lambda }_j({\sigma }^2)-\sum _{i=1}^n\sum _{j=0}^n{\lambda }_i{\lambda }_j(r_{ij})-2\sum _{i=1}^n{\lambda }_i({\sigma }^2)+2\sum _{i=1}^n{\lambda }_i(r_{io})+{\sigma }^2–r_{oo}\end{array}$$
考虑到${\sum }_{i=1}^n{\lambda }_i=1$
$$\begin{array}{rl}J& ={\sigma }^2-\sum _{i=1}^n\sum _j^n{\lambda }_i{\lambda }_j(r_{ij})-2{\sigma }^2+2\sum _{i=1}^n{\lambda }_i(r_{io})+{\sigma }^2–r_{oo}\\ & =2\sum _{i=1}^n{\lambda }_i(r_{io})-\sum _{i=1}^n\sum _{j=0}^n{\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 }_{i=1}^n{\lambda }_i=1$，这是一个带约束条件的最优化问题。使用拉格朗日乘数法求解，求解方法为构造一个新的目标函数
$$J+2\varphi (\sum _{i=1}^n{\lambda }_i-1)$$

其中$\varphi$是拉格朗日乘数。求解使这个代价函数最小的参数集$\varphi ,{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$，则能满足其在 ${\sum }_{i=1}^n{\lambda }_i=1$约束下最小化$J$。即
$$\{\begin{array}{rl}\frac{\partial (J+2\varphi ({\sum }_{i=1}^n{\lambda }_i-1))}{\partial {\lambda }_k}& =0;k=1,2,\cdots ,n\\ \frac{\partial (J+2\varphi ({\sum }_{i=1}^n{\lambda }_i-1))}{\partial \varphi }& =0\end{array}$$
$$\{\begin{array}{rl}\frac{\partial (2{\sum }_{i=1}^n{\lambda }_i(r_{io})–{\sum }_{i=1}^n{\sum }_j^n{\lambda }_i{\lambda }_j(r_{ij})–r_{oo}+2\varphi ({\sum }_{i=1}^n{\lambda }_i-1))}{\partial {\lambda }_k}& =0;k=1,2,\cdots ,n\\ \frac{\partial (2{\sum }_{i=1}^n{\lambda }_i(r_{io})–{\sum }_{i=1}^n{\sum }_j^n{\lambda }_i{\lambda }_j(r_{ij})–r_{oo}+2\varphi ({\sum }_{i=1}^n{\lambda }_i-1))}{\partial \varphi }& =0\end{array}$$
$$\{\begin{array}{rl}2r_{ko}–\sum _{j=1}^n(r_{kj}+r_{jk}){\lambda }_j+2\varphi & =0;k=1,2,\cdots ,n\\ \sum _{i=1}^n{\lambda }_i& =1\end{array}$$
由于 $C_{ij}=Cov(z_i,z_j)=C_{ji}$，因此同样地 $r_{ij}=r_{ji}$，那么有
$$\{\begin{array}{rl}{r}_{ko}–\sum _{j=1}^n{r}_{kj}{\lambda }_j+\varphi & =0;k=1,2,\cdots ,n\\ \sum _{i=1}^n{\lambda }_i& =1\end{array}$$
式子中半方差函数$r_{ij}$十分重要，最后会详细解释其计算与定义 在以上计算中我们得到了对于求解权重系数${\lambda }_j$的方程组。写成线性方程组的形式就是：
$$\{\begin{array}{rl}{r}_{11}{\lambda }_1+r_{12}{\lambda }_2+\cdots +r_{1n}{\lambda }_n–\varphi & =r_{1o}\\ r_{21}{\lambda }_1+r_{22}{\lambda }_2+\cdots +r_{2n}{\lambda }_n–\varphi & =r_{2o}\\ & \cdots \\ r_{n1}{\lambda }_1+r_{n2}{\lambda }_2+\cdots +r_{nn}{\lambda }_n–\varphi & =r_{no}\\ {\lambda }_1+{\lambda }_2+\cdots +{\lambda }_n& =1\end{array}$$

写成矩阵形式即为
$$\left[\begin{array}{ccccc}{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{array}\right]\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ \cdots \\ {\lambda }_n\\ -\varphi \end{array}\right]=\left[\begin{array}{c}{r}_{1o}\\ r_{2o}\\ \cdots \\ r_{no}\\ 1\end{array}\right]$$
对矩阵求逆即可求解。 唯一未知的就是上文中定义的半方差函数 $r_{ij}$，接下来将详细讨论

注解：关于求解矩阵 $A\lambda =b$，使用cholesky分解比直接求逆快

1.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}{rl}{r}_{ij}& =\frac{1}{2}E[(R_i-R_j{)}^2]\\ & =\frac{1}{2}E[R_i^2-2R_i{R}_j+R_j^2]\\ & =\frac{1}{2}E[R_i^2]+\frac{1}{2}E[R_j^2]-E[R_i{R}_j]\end{array}$$

又因为：
$$E[R_i^2]=E[R_j^2]=E[(z_i–c{)}^2]=Var(z_i)={\sigma }^2$$
$$E[R_i{R}_j]=E[(z_i–c)(z_j-c)]=Cov(z_i,z_j)=C_{ij}$$
于是有$$\begin{array}{rl}{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_i{R}_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之间的几何距离

克里金插值假设 $r_{ij}$与 $d_{ij}$存在着函数关系，这种函数关系可以是线性、二次函数、指数、对数关系。为了确认这种关系，我们需要首先对观测数据集
{ 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_{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}$，然后根据得到的函数关系就可以得到这两点的半方差.

1.7半方差模型

线性有基台模型：

球状模型：

指数模型：

高斯模型：

不同模型之间可以线性叠加

1.8小结

总的来说，进行克里金插值分为这几个步骤：

1. 对于观测数据，两两计算距离与半方差
2. 寻找一个拟合曲线拟合距离与半方差的关系，从而能根据任意距离计算出相应的半方差
3. 计算出所有已知点之间的半方差 $r_{ij}$
4. 对于未知点 $z_o$，计算它到所有已知点 $z_i$的半方差$r_{io}$
5. 求解第四节中的方程组，得到最优系数 ${\lambda }_i$
6. 使用最优系数对已知点的属性值进行加权求和，得到未知点 $z_o$的估计值

2.Py实现

GeostatsGuy Lectures： https://github.com/GeostatsGuy/PythonNumericalDemos/blob/master/GeostatsPy_kriging.ipynb.

2.1 小结

待补充。。。