克里金(克里格)(Corigine)算法.pdf

克里格，或者说克里金插值Kriging。法国krige名字来的。

特点是线性，无偏，方差小，适用于空间分析。所以很适合地质学、气象学、地
理学、制图学等。

相对于其他插值方法。主要缺点：由于他要依次考虑（这也是克里格插值的一般
顺序）计算影响范围，考虑各向异性否，选择变异函数模型，计算变异函数值，
求解权重系数矩阵，拟合待估计点值，所以反映速度很慢。（当然也看你算法设
计和电脑反应速度了呵呵）。而那些趋势面法，样条函数法等。虽然较快，但是
毕竟程度和适合用范围都大受限制。

具体

对比如下：

方法外推能
力逼近程度运算能力适用范围

距离反比加权法分布均匀时
好差快分布均匀

最近邻点插值法不
高强很
快分布均匀

三角网线性插
值高差
慢分布均匀

样条函
数高
强快分布密集时候

克里金插
值高强
慢均可

克里格插值又分为：简单，普通，块，对数，指示性，泛，离析克里金插值等。

克里金插值的变异函数球形模型，指数模型，高斯模型，纯块金模型，幂函数模
型，迪维生模型等。

以下结合我的绘制等值线（等高线）的程序和高斯迭代解矩阵方程方法以及多元
线性回归方法（此两方法实现另补充）说明克里格方法的实现：
注：选择变异函数模型为球形模型，选择插值方法为普通克里金，我为了简化问
题，考虑为各向同性，变差距离为固定。

inti,j,i0,i1,j0,j1,k,l,m,n,p,h;//循环变量
double*r1Matrix;//系数矩阵
double*r0Matrix;//已知向量
double*langtaMatrix;//待求解向量
double*x0;//已知点横坐标
double*y0;//已知点纵坐标
double*densgridz;//存储每次小方格内的已知值。
doubledensgridz0;//待求值
intN1=0;//统计有多少个已知值
doubler[71],r0[71];
intN[70];
for(i=0;i<100;i++)
{
for(j=0;j<100;j++)
{

if(bdataprotected[i*100+j])continue;//原值点不需要插值

//1.遍历所有非保护网格。确定每一个待插值点的r(h)
//每一个网格又从横向和纵向进行搜索，也就是说正方形相关，正方形的边
长以R，格子长度为50；中心距离为25

//首先计算起循环的起始点。
//横向
if(i-25>=0)
i0=i-25;
elsei0=0;
if(i+25<=100)
i1=i+25;
elsei1=100;
//纵向
if(j-25>=0)
j0=j-25;
elsej0=0;
if(j+25<=100)
j1=j+25;
elsej1=100;

//Hmax=int(50*2^.5)=70根据对称性，所有的r(h)除以2即为所得值。
//先待插值点的编程小方格内统计有几个已知点，如果个数小于4，则
不能拟合。
N1=0;
for(l=i0;l<i1;l++)
for(m=j0;m<j1;m++)
{
if(bdataprotected[100*l+m])N1++;
}
if(N1<4)continue;//不符合线性回归条件，本网格不用此方法做插
值
//赋初值
for(l=0;l<=70;l++)
{
r0[l]=0.0;
r[l]=0.0;
}
//计算此插值点方格内的变异函数值
for(intl=i0;l<i1;l++)
for(m=j0;m<j1;m++)
if(i!=l&&j!=m)//不计待估计值本身
{
//自循环统计算网格内部的互相之间的h，N和z。
for(n=i0;n<i1;n++)
for(p=j0;p<j1;p++)
if(l!=n&&p!=m)//不对h=0计算
{
if(bdataprotected[l*100+m]&&bdataprotected[n*100+p])//保证用
样本值而非估计值来进行估计
{
//计算分离距离
h=int(sqrt(l-n)*(l-n)+(m-p)*(m-p));
//计算不同分离距离的样本差的平方和
r0[h]=r0[h]+(densgrid[l*100+m]-densgrid[
n*100+p])*(densgrid[l*100+m]-densgrid[n*100+p]);
//不同的分离距离的样本对数
N[h]++;
}
}
}

//不同分离距离的变异函数值
for(h=1;h<=70;h++)
r[h]=r0[h]/(2*N[h])/2;
//2.通过所有的r(h)拟合计算球形模型的参数。方法为多元线性回归
//参数初始化
doublex1[70],y1[70];intn1=70;doublea1[4]={0,0,0