最近因为项目需要,研究了下克里金插值算法。在地质学中,克里金插值算法是一种使用的空间属性估计技术,克里金插值说到底是个回归问题,且依据的因素只有两个位置之间的距离。
克里金插值算法又分为很多中,比如普通克里金插值,简单克里金插值等,不同的克里金插值算法只是假设条件不同。下面以普通克里金为例来说明其原理。
普通克里金插值算法的假设条件是:空间属性是均匀的,对于空间任意一点(x,y)都有同样的数学期望和方差。普通克里金插值的目的是对于空间一点(x,y)找到一组系数
,估计其属性值为:
![](https://i-blog.csdnimg.cn/blog_migrate/e938833f2c3b0f6b7b8fbf9ed6787a98.png)
![](https://i-blog.csdnimg.cn/blog_migrate/e7d579f337ae8a8bdef9576e0a6f18a2.png)
![](https://i-blog.csdnimg.cn/blog_migrate/c7c55220fd2ba1edc294439a3d131329.png)
![](https://i-blog.csdnimg.cn/blog_migrate/79113c8be244d1a2d35aae49d6c19b35.png)
无偏性:
![](https://i-blog.csdnimg.cn/blog_migrate/92175e6fb99c82346ba48fc91c3eb171.png)
最优性:
![](https://i-blog.csdnimg.cn/blog_migrate/80c9206e44b709405e9ac5e93c511a71.png)
参考:
https://blog.csdn.net/sinat_23619409/article/details/79974990,易知,无偏性可以推出
,最终问题转化为一个带约束条件的最优化问题:
![](https://i-blog.csdnimg.cn/blog_migrate/e41327df1ab7e75c7546175a1f1217e5.png)
最终,可以写成:
![](https://i-blog.csdnimg.cn/blog_migrate/d4b13c671cc123f5dba430086256e639.png)
定义
为半方差函数,在地质学上,一般
,所以,上面的等号左边的
都是已知的,因为属性值是已知的。所以逆也是可以求得的。只有能够求得右边的
,那么系数
便可以求得,普通克里金插值任务也就完成了。
![](https://i-blog.csdnimg.cn/blog_migrate/7c387f3a352b7412f05dba6569d160ed.png)
![](https://i-blog.csdnimg.cn/blog_migrate/9ff4045fec1c19e454d429a2130d17c9.png)
![](https://i-blog.csdnimg.cn/blog_migrate/9b4a794a264f330061519534a6691d75.png)
![](https://i-blog.csdnimg.cn/blog_migrate/1d9ba11603230dc1469e33bc73ce032d.png)
![](https://i-blog.csdnimg.cn/blog_migrate/db42bed9ea5249e16cc429571254e145.png)
如何求
?对所有已知属性的空间点,求两两的距离以及半方差函数值。用一个函数式模拟二者的关系,可以是多项式关系,可以是线性关系,可以是指数,对数关系。
![](https://i-blog.csdnimg.cn/blog_migrate/e6045a36d7e0160fff8c402a5e21b79e.png)
然后用高斯牛顿迭代法,进行迭代,求出拟合系数。
用拟合的关系式,求出
,问题就可以解决。理论上看上去不难。我在算法实现过程中发现几个问题:
![](https://i-blog.csdnimg.cn/blog_migrate/4dbe9f83b02fab13b653513c450c3f38.png)
1)我所用到的已知点的属性值都很接近,
有很多都是0或者接近0的,所以
![](https://i-blog.csdnimg.cn/blog_migrate/ea889ac05d6aec6d4adfc99fbce119ac.png)
![](https://i-blog.csdnimg.cn/blog_migrate/7e3f8745c9f8642c7af3aaa8ca062b75.png)
是一个稀疏矩阵,逆不存在,即使用伪逆最后的结果也很怪异。最后我采取的办法是求
也是用拟合的关系式求取,虽然和实际可能有些偏差,但是总比无法求解强。
![](https://i-blog.csdnimg.cn/blog_migrate/7e3f8745c9f8642c7af3aaa8ca062b75.png)
2)距离和半方差值得关系式对结果至关重要,既要选取和已知的观察点比较吻合的关系,表达式也不能太复杂,避免高斯牛顿迭代法计算复杂或者无法收敛。
参考: