克里金插值法的java实现

项目中要用到克里金插值法，大致了解了一下，今天做个笔记总结一下（有错误请评论指正）

关于克里金插值法，在我看来就是加强版的反距离加权，只不过他的权重系数的确定，复杂一点，是带着你自己的空间模型的分布特性，比如说你要用在气象领域，则权重系数和地质的就是完全不相同的。

我对于克里金方法的理解，认为他的算法可以分成五步（前提是你的模型已确定），第一步是求出每两个已知点之间的距离A，然后带入模型算法中求出对应值，组成已知点的距离矩阵，第二步是求出预测点与已知点的距离，代入模型算法，得到矩阵B，第三步就是通过A的逆阵*B，求出权重系数矩阵，第四步就是通过模型计算每个点的z,第五步就是通过z值和权重系数进行相乘相加，得到最终的预测值。

模型算法的求解，就是根据拟合的量度，克里金算法总共有六个，还是七个模型，可以以距离的半方差当作x,相应的value对应为y，然后查看与模型的拟合程度，确定模型的参数。

我这里是PM2.5的分布预测，在网上找了个关于gri3d的包，但是封装较为严重，想看明白里面的参数配置有点难，就自己写了一个，但是没有考虑一些克里金方面的极值情况，比如h=0等，只是实现了大致的逻辑部分，模型我选择的是幂指数模型，在我看来，选模型较为困难，我自己用python画图求拟合，表示根本拟合不了（可能是样本的原因吧），方程中的参数是算法同事给我的。废话不多说，程序如下：

public double[][] gridData(double[] x, double[] y, double[] z,
                           double[][] X, double[][] Y) {
    int xLength = x.length;
    int xSize = X.length;
    int ySize = X[0].length;
    double[][] data = new double[xLength + 1][xLength + 1];
    double[][] data2 = new double[xSize][ySize];

    for (int i = 0; i <= xLength; i++) {
        for (int j = 0; j <= xLength; j++) {
            if ((i < xLength) && (j < xLength)) {

                if (i != j) {
                    data[i][j] = Math.sqrt((x[i] - x[j]) * (x[i] - x[j])
                            + (y[i] - y[j]) * (y[i] - y[j]));  //算出每两点之间的距离
                } else {
                    data[i][j] = 0.0001;
                }
            } else if ((i == xLength) || (j == xLength)) {
                data[i][j] = 1;
            }
            if (i == xLength && j == xLength) {
                data[i][j] = 0;
            }
        }
    }
    //矩阵的初始化
    Matrix m = new Matrix(data);
    //矩阵求逆  逆矩阵矩阵a*a-1 = 单位矩阵
    Matrix mMatrix = m.inverse();

    for (int i = 0; i < xSize; i++) {
        for (int j = 0; j < ySize; j++) {
            double[] rData = new double[xLength + 1];

            for (int k = 0; k <= xLength; k++) {
                if (k < xLength) {
                    rData[k] = Math.sqrt((X[i][j] - x[k])
                            * (X[i][j] - x[k]) + (Y[i][j] - y[k])
                            * (Y[i][j] - y[k])); //算出预测值的距离
                } else {
                    rData[k] = 1;  //预测值的距离
                }
            }

            double[] rDa = new double[xLength + 1];
            for (int k = 0; k < (xLength + 1); k++) {
                for (int k2 = 0; k2 < (xLength + 1); k2++) {
                    rDa[k] += mMatrix.get(k, k2) * rData[k2];  //返回逆矩阵*rData距离 求权重
                }
            }
            for (int k = 0; k < xLength; k++) {
                data2[i][j] += z[k] * rDa[k];
            }
        }
    }

    return data2;
}

// 变差函数，克里金算法特有的函数
public double function(double h) {
    double r = 0;
    if (h == 0) {
        r = 0;
    } else if (h > 0 && h <= 11) { //此处是他的基程
        Kriging kriging = new Kriging();
        //kriging.variogramModel()
        //r = 0.0 + 1.154 * ((3 * h * h * h * h) / (4 * 8.535 * 8.535 * 8.535 * 8.535));
        r = 0.0 + 0.06 * (1.0D - Math.exp(- h / 11));

    } else {
        r = 3.202;
    }
    return r;

}

幂指数模型也有是-3h/a，感觉就是凭个人意愿吧，无非是a*3就行了。

有问题欢迎评论区指正，谢谢！！

上面的这个基本上就是完整的代码了，在测试时发现了一个问题，就是插值只能一个一个插，如果插太多值的话，插值结果会有问题。

附我的测试代码：

public static void main(String[] args) {
    double[] Lon = {116.376259, 116.376259, 116.376259, 116.376259, 116.376259, 116.376259, 116.376259, 116.375824, 116.375809, 116.375839, 116.375916, 116.375786, 116.375748, 116.375786, 116.375801, 116.375885, 116.37632, 116.37632, 116.37632, 116.374489, 116.374489, 116.375938, 116.376144, 116.375999, 116.376152
            };
    double[] Lat = {39.971764, 39.971764, 39.971764, 39.971764, 39.971764, 39.971764, 39.971764, 39.970943, 39.971077, 39.971268, 39.970734, 39.970695, 39.970432, 39.970318, 39.970619, 39.971119, 39.971336, 39.971336, 39.971336, 39.968445, 39.968445, 39.971176, 39.971252, 39.971279, 39.971516
    };
    double[] PM = {57.576897869, 57.576897869, 57.576897869, 57.576897869, 57.576897869, 57.576897869, 57.576897869, 60.347131145, 61.732247783, 63.80992274, 61.039689464, 63.117364421, 58.962014507, 62.424806102, 61.732247783, 61.039689464, 62.424806102, 60.347131145, 57.576897869, 60.347131145, 60.347131145, 63.117364421, 58.269456188, 61.039689464, 56.191781231//, 54.806664593

    };
    double[] lon = {116.376183};
    double[] lat =  {39.97163};
    double[] m = gridData(Lon, Lat, PM, lon, lat);
    System.out.println(Arrays.toString(m));
}

54是真实值，57是插的值，感觉还行吧，因为设备的值应该是发生了突变