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

克里格，或者说克里金插值Kriging。

法国krige名字来的。

特点是线性，无偏，方差小，适用于空间分析。

所以很适合地质学、气象学、地理学、制图学等。

相对于其他插值方法。

主要缺点：由于他要依次考虑（这也是克里格插值的一般顺序）计算影响范围，考虑各向异性否，选择变异函数模型，计算变异函数值，求解权重系数矩阵，拟合待估计点值，所以反映速度很慢。

（当然也看你算法设计和电脑反应速度了呵呵）。

而那些趋势面法，样条函数法等。

虽然较快，但是毕竟程度和适合用范围都大受限制。

具体对比如下：方法外推能力逼近程度运算能力适用范围距离反比加权法分布均匀时好差快分布均匀最近邻点插值法不高强很快分布均匀三角网线性插值高差慢分布均匀样条函数高强快分布密集时候克里金插值高强慢均可克里格插值又分为：简单，普通，块，对数，指示性，泛，离析克里金插值等。

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

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

int i,j,i0,i1,j0,j1,k,l,m,n,p,h;//循环变量double *r1Matrix;//系数矩阵double *r0Matrix;//已知向量double *langtaMatrix;//待求解向量double *x0;//已知点横坐标double *y0;//已知点纵坐标double * densgridz;//存储每次小方格内的已知值。

double densgridz0;//待求值int N1=0;//统计有多少个已知值double r[71],r0[71];int N[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;else i0=0;if(i+25<=100)i1=i+25;else i1=100;//纵向if(j-25>=0)j0=j-25;else j0=0;if(j+25<=100)j1=j+25;else j1=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(int l=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)拟合计算球形模型的参数。

方法为多元线性回归//参数初始化double x1[70],y1[70];int n1=70; double a1[4]={0,0,0,0};//a1数组存储拟合结果int m1=4;double dt1[4];//dt1数组时拟合精度double b0,b1,b2,b3;double c0,c1,a;for(h=0;h<70;h++){x1[h]=double(h+1.0);y1[h]=r[h+1];//从h=1开始有数据}//带入多元线性回归拟合模型MLR(x1,y1,n1,a1,m1,dt1);//求得球形模型的参数如下。

注意x也即h的平均值为1到70中间数。

35.5； b0=a1[0]-a1[1]*35.5+a1[2]*35.5*35.5-a1[3]*35.5*35.5*35.5;b1=a1[1]-2*35.5*a1[2]+3*a1[3]*35.5*35.5;//b2==0 忽略b3=a1[3];//化为球形模型的标准形式参数参见空间变异理论及应用 p16 及空间插值技术开发及实现 //2.3.1 2c0=b0;a=sqrt(-b1/(3*b3));c1=3/(2*a*b1);//3.求解权重系数矩阵方程的解//开辟空间r1Matrix=new double[(N1+1)*(N1+1)];r0Matrix=new double[N1+1];langtaMatrix=new double[N1+1];x0=new int[N1];y0=new int[N1];densgridz=new double[N1];//求取已知点的横纵坐标及其值int i0=0;for(l=i0;l<i1;l++)for(m=j0;m<j1;m++){if(bdataprotected[100*l+m]){x0[i0]=gridx[100*l+m];y0[i0]=gridy[100*l+m];densgridz[i0]= densgridz[100*l+m];i0++;}}//给系数矩阵赋值for(int ii=0;ii<N1+1;ii++)for(int jj=0;jj<N1+1;jj++){//最后一行最后一列分别赋值if(ii==N1) r1Matrix[ii*(N1+1)+jj]=1;else if(jj=N1) r1Matrix[ii*(N1+1)+jj]=1;else if(ii==N1 && JJ==N1) r1Matrix[ii*(N1+1)+jj];else//左上三角{//求解已知点的hh=int(sqrt((x0[ii]-x0[jj])* (x0[ii]-x0[jj])+ (y0[ii]-y0[jj])* (y0[ii]-y0[jj])));//计算rij=r(h)r1Matrix[ii*(N1+1)+jj]=b0+b1*h+b3*h*h*h;}}//已知向量矩阵赋值ri0for(ii=0;ii<N1;ii++){//计算hh=int(sqrt((x0[ii]-gridx[i*100+j])* (x0[ii]-gridx[i*100+j])+ (y0[ii]-gridy[100*i+j])* (y0[ii]-gridy[100*i+j])));r0Matrix[ii]= b0+b1*h+b3*h*h*h;}//最后一个单独赋值r0Matrix[N1]=1;double eps=0.0001;N1=N1+1；//用求解矩阵,用高斯迭代法if (gaos(N1, r1Matrix,r0Matrix, langtaMatrix,eps)>0){//得到待估计点数值for(ii=0;ii<N1;ii++)densgridz0= densgridz0+densgridz[ii]*langtaMatrix[ii]; //赋值给对应网格densgrid[i*100+j]= densgridz0;}else continue;//如果无解，则继续//释放内存空间if( r1Matrix!=NULL) delete r1Matirx;if(r0Matrix!=NULL) delete r0Matirx;if(langtaMatrix!=NULL) delete langtaMatrix; if(x0!=NULL) delete x0;if(y0!=NULL) delete y0;if(densgridz!=NULL) delete densgirdz;}。