kringing.txt

数值分析算法学习笔记之克里格插值及其实现2009-04-21 22:25克里格

数值分析算法学习笔记之克里格插值及其实现2009-04-21 22:25克里格，或者说克里金插值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 2
        c0=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//左上三角
     {
      //求解已知点的h
      h=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;
     }
    }
   //已知向量矩阵赋值ri0
   for(ii=0;ii<N1;ii++)
   {
    //计算h 
    h=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;
   
}