通过C#代码实现空间离散点的克里金(kriging)插值

1. 求取已知点的距离以及点对的半方差
2. 筛选第一步求取的结果，计算出几个均值点，用于拟合
3. 选定拟合模型，为了方便代码实现，我选择了指数模型
4. 用指数模型去拟合第二步得出的均值点，得出偏基台值 c” role=”presentation” style=”position: relative;”>cc
5. 根据拟合得到的模型，按照公式通过已知点高程计算位置点高程

点对的半方差计算公式：
Semivariogram(distanceh)=12∗(valuei–valuej)2” role=”presentation” style=”position: relative;”>Semivariogram(distanceh)=12∗(valuei–valuej)2Semivariogram(distanceh)=12∗(valuei–valuej)2
计算结果为一堆横坐标为distanceh” role=”presentation” style=”position: relative;”>distancehdistanceh的点。
将这些点按照横坐标分为约10份，（此处可按照实际情况修改），计算每个区间里的点的横坐标纵坐标均值，这样就能得到10个点代表刚刚的整个计算结果。
如果点数过多，拟合的效果以及效率会受到影响。

以下提到的数学参数请对照《计算原理》一文中的定义。

定义方法：

private double CalCij(double x1, double y1, double x2, double y2)
{
    double distance = Math.Sqrt(
        Math.Pow(x1 - x2, 2) + Math.Pow(y1 - y2, 2));
    if (distance == 0)
    {
        return formula_c;
    }
    else
    {
        return formula_c* Math.Exp(-distance / formula_r);
    }
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13

计算两点之间的Cij” role=”presentation” style=”position: relative;”>CijCij 就是模型结果中的一个值。
计算矩阵 K” role=”presentation” style=”position: relative;”>KK：

//size为已知点的个数
var K = new DenseMatrix(size, size);
for (int m = 0; m < size; m++)
        for (int n = 0; n < size; n++)
            K[m, n] = CalCij(
             pointList[m].X,
             pointList[m].Y, 
             pointList[n].X,
             pointList[n].Y);

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9

计算矩阵 K” role=”presentation” style=”position: relative;”>KK：

Kn = K.Inverse();

• 1

假设求坐标为 (m,n)” role=”presentation” style=”position: relative;”>(m,n)(m,n) 的未知点的高程：
1. 计算向量 D” role=”presentation” style=”position: relative;”>DD：

var D = new DenseVector(size);
for (int p = 0; p < size; p++)
    D[p] = CalCij(randomPointList[p].point.X,
        randomPointList[p].point.Y, m, n);

• 1
• 2
• 3
• 4

1. 计算λ(i)” role=”presentation” style=”position: relative;”>λ(i)λ(i) 个已知点对当前未知点的影响权重：

var namuta = Kn.LeftMultiply(D);

• 1

1. 计算 Z(xi)” role=”presentation” style=”position: relative;”>Z(xi)Z(xi) 个点的高程值：

for (int q = 0; q < size; q++)
    interpolationDEMData[m, n] += 
        namuta[q] * randomPointList[q].altitudeValue;

• 1
• 2
• 3

至此，坐标为 (m,n)” role=”presentation” style=”position: relative;”>(m,n)(m,n) 数组中。
下图为插值结果：

结果与原始数据相比较，还是比较准确的。