文章目录

• • @[toc]
• 前言
• 一、克里金插值原理
• • 1.1 概述
  • 1.2 基本公式
  • 1.2 权重$w_i$的确定
  • 1.3 拟合函数的确定
• 二、Python建模与可视化
• • 2.1 Demo
  • • 2.1.1 随机生成已知格网点
    • 2.1.2 拟合
    • 2.1.3 评估内符合精度
    • 2.1.3 内插未知格网点
    • 2.1.4 画图
  • 2.2 结果图
• 参考文献

前言

最近学习了一下克里金插值的原理并做了一点简单尝试，在此文记录！

本文详细介绍克里金(Kriging)插值的原理和python实现。首先介绍克里金插值的原理及理解，再详细介绍基于pykrige库实现的克里金插值Demo，并使用随机生成的数据进行了测试。

---

一、克里金插值原理

1.1 概述

克里金插值是根据空间相关关系通过一系列已知点的属性预测相近的未知点的属性的插值方法。有两点假设。
假设1：一点的属性值与其周围点的属性值有关，并且可以由其周围点的属性值推导出。
假设2：两点属性值差异性(不相关性)与二者间距离在一定距离范围内成正相关（核心思想）。

1.2 基本公式

已知条件：

• 已知点坐标及属性值
• 已知待求点坐标
• 已知点-待求点距离矩阵
• 已知点间距离矩阵

根据假设1，空间内某未知值点估计值$z_0^{\prime }$是部分已知点值的加权和

$z_0^{\prime }={\sum }_{i=1}^n{z}_i{w}_i$

其中n是可以决定该点属性值的点数量(一般取该点周边一定数量的已知点)，
$z_i$代表第i个已知点的属性值，$w_i$代表第i个已知点的权重。

1.2 权重$w_i$的确定

次级假设1：
为得到尽可能准确的结果，$w_i$应使得点处估计值$z_0^{\prime }$与真实值$z_0$之差最小。

$E\left(z_0^{\prime }-z_0\right)=0$

同时我们希望估计值$z_0^{\prime }$与真实值$z_0$之差的方差尽可能小，即

$\min \left(\mathrm{Var}\left(z_0^{\prime }-z_0\right)\right)$

次级假设2：
空间是平稳的，空间任意一点处的值z=z(x,y)由区域平均值c和随
机偏差R(x,y)组成，其中偏差的方差均为常数。

$E(z)=c$
$Var[R(x,y)]={\sigma }^2$

权重计算公式：

$\left[\begin{array}{c}{w}_1\\ ⋮\\ w_n\\ \lambda \end{array}\right]={\left[\begin{array}{ccc}{\gamma }_{11}& \cdots & {\gamma }_{1n}\\ ⋮& \ddots & ⋮\\ {\gamma }_{n1}& \cdots & {\gamma }_{nn}\\ 1& 1& 0\end{array}\right]}^{-1}\left[\begin{array}{c}{\gamma }_{01}\\ ⋮\\ {\gamma }_{0n}\\ 1\end{array}\right]$
其中，i和j代表不同的领域点，0代表待求点。${\gamma }_{ij}$为半方差，可由${\gamma }_{ij}=\frac{{\left(z_{\mathrm{i}}-z_{\mathrm{j}}\right)}^2}{2}$计算。所以上式矩阵的系数已知，但右侧的向量未知，可通过拟合函数$\gamma =f(d)$计算。

1.3 拟合函数的确定

克里金插值常用的拟合函数有：

式中，h就是距离d，C0和C1是带确定的系数，a为假设2存在空间相关性的有限距离。
为了排除偶然性，确定拟合函数的系数时一般先分组计算距离和半方差平均值，而不是输入原始的距离-半方差对计算。

二、Python建模与可视化

2.1 Demo

2.1.1 随机生成已知格网点

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
from mpl_toolkits.basemap import Basemap

# 定义边界
lonLeft = 100
lonRight = 139
latDown = 20
latUp = 60

# 随机生成数据，即已知格网点
np.random.seed(42)
num_points = 50
lons = np.random.uniform(lonLeft, lonRight, num_points)
lats = np.random.uniform(latDown, latUp, num_points)
data = np.sin(lons) + np.cos(lats) + np.random.normal(0, 0.1, num_points)

2.1.2 拟合

variogram_model即拟合函数，可选参数有linear, power, gaussian, spherical, exponential, hole-effect.

from pykrige.ok import OrdinaryKriging
import matplotlib
from mpl_toolkits.basemap import Basemap

OK = OrdinaryKriging(lons, lats, data, variogram_model = 'gaussian', nlags=12)

2.1.3 评估内符合精度

## 评估拟合精度
dataPrd, ss1 = OK.execute('points', lons, lats)
compressed_array = dataPrd.compressed()
result_list_from_compressed = compressed_array.tolist()
squared_diffs = [(result_list_from_compressed[iLoop] - data[iLoop]) ** 2 for iLoop in range(len(data))]

# 求和
sum_squared_diffs = sum(squared_diffs)
 
# 取平均
mean_squared_diffs = sum_squared_diffs / len(data)
 
# 开方
rmse = math.sqrt(mean_squared_diffs)
print(rmse)

1.83099255103024e-15

2.1.3 内插未知格网点

# 生成未知格网点
grid_lon = np.linspace(100,139,400)
grid_lat = np.linspace(20,60,400)

## 内插未知点
z1, ss1 = OK.execute('grid', grid_lon, grid_lat)

2.1.4 画图

#转换成网格
xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)
#将插值网格数据整理
df_grid =pd.DataFrame(dict(long=xgrid.flatten(),lat=ygrid.flatten()))
#添加插值结果
df_grid["Krig_gaussian"] = z1.flatten()


fig,ax = plt.subplots(figsize=(6,4.5),dpi=130,facecolor="white")
map_base = Basemap(lonLeft,latDown, lonRight, latUp,
                  projection="cyl",lon_0 = 119,lat_0 = 33,ax = ax)
# lat = np.arange(30,36,1)
# lon = np.arange(116,122,1)
map_base.drawparallels(np.arange(30,36,1), labels=[1,0,0,0],fontsize=12,zorder=0) #画纬度线
map_base.drawmeridians(np.arange(116,122,1), labels=[0,0,0,1],fontsize=12,zorder=0) #画经度线
# map_base.readshapefile( name = "Js", default_encoding="ISO-8859-1",
#                        drawbounds=True)
cp=map_base.pcolormesh(xgrid, ygrid, data=z1.data,cmap='Spectral_r')  
colorbar = map_base.colorbar(cp,size='3%',pad="5%",label="Kriging_inter")
#设置colorbar
colorbar.outline.set_edgecolor('none')
 
for spine in ['top','left','right','bottom']:
    ax.spines[spine].set_visible(None) #隐去轴脊
 
ax.text(.5,1.1,"Map Charts in Python Exercise 02:Map Kriging Grid line",transform = ax.transAxes,ha='center', 
        va='center',fontweight="bold",fontsize=14)
ax.text(.5,1.03, "processed map charts with Basemap",
        transform = ax.transAxes,ha='center', va='center',fontsize = 10,color='black')
ax.text(.83,-.06,'\nVisualization by DataCharm',transform = ax.transAxes,
        ha='center', va='center',fontsize = 8,color='black')

2.2 结果图

50个已知点：

100个已知点：

参考文献

【GIS算法】克里金插值原理详解
https://www.bilibili.com/video/BV1bT4y1C7z6

空间绘图 | Python-pykrige包-克里金(Kriging)插值计算及可视化绘制tps://blog.csdn.net/qq_40483688/article/details/135821334

python克里金插值建模
https://blog.51cto.com/u_16175499/12763524