用Python从零构建克里金插值模型告别公式恐惧的实战指南当第一次接触克里金插值时大多数人都被那堆晦涩的协方差矩阵和最大似然估计吓退。但今天我们要用Python和NumPy这把瑞士军刀把看似高深的理论拆解成可运行的代码块——你会发现真正理解克里金的最好方式就是亲手实现它。1. 环境准备与数据生成工欲善其事必先利其器。我们从一个干净的Python环境开始import numpy as np import matplotlib.pyplot as plt from scipy.spatial.distance import cdist np.random.seed(42) # 固定随机种子确保可复现1.1 生成模拟观测数据克里金插值需要已知点数据作为基础。我们创建包含20个随机观测点的数据集# 生成二维空间中的随机观测点 obs_points np.random.rand(20, 2) * 10 # 为每个点生成模拟值这里用正弦函数制造空间相关性 obs_values np.sin(obs_points[:,0]/3) * np.cos(obs_points[:,1]/3) np.random.normal(0, 0.1, 20)用Matplotlib可视化这些观测点plt.scatter(obs_points[:,0], obs_points[:,1], cobs_values, cmapviridis, s100) plt.colorbar(label观测值) plt.title(空间观测点分布) plt.show()2. 变异函数空间相关性的量化工具变异函数(variogram)是克里金的核心它量化了空间两点间的差异随距离变化的规律。2.1 计算经验变异函数定义计算经验变异函数的函数def empirical_variogram(points, values): 计算经验变异函数 :param points: 观测点坐标数组 (n,2) :param values: 观测值数组 (n,) :return: 距离数组, 变异函数值数组 # 计算所有点对间的距离 dist_matrix cdist(points, points) # 计算所有点对的值差异平方 value_diff (values[:, None] - values[None, :])**2 # 提取上三角矩阵避免重复计算 upper_tri np.triu_indices_from(dist_matrix, k1) distances dist_matrix[upper_tri] gamma 0.5 * value_diff[upper_tri] return distances, gamma2.2 拟合理论变异函数模型常见理论模型有指数型、高斯型和球型。我们实现指数模型def exponential_variogram(h, nugget, sill, range_): 指数变异函数模型 return nugget (sill - nugget) * (1 - np.exp(-3 * h / range_))通过非线性最小二乘法拟合参数from scipy.optimize import curve_fit # 计算经验值 distances, gamma empirical_variogram(obs_points, obs_values) # 拟合理论模型 popt, _ curve_fit(exponential_variogram, distances, gamma, p0[0.1, 1.0, 5.0], # 初始猜测块金值、基台值、变程 bounds(0, [10, 10, 20])) nugget, sill, range_ popt print(f拟合参数块金值{nugget:.3f}, 基台值{sill:.3f}, 变程{range_:.3f})3. 构建克里金矩阵系统克里金插值的本质是求解一个权重系统使预测值最优。3.1 构建克里金矩阵实现普通克里金(Ordinary Kriging)的矩阵构建def build_kriging_matrix(points, variogram_func): 构建克里金矩阵系统 :param points: 观测点坐标 (n,2) :param variogram_func: 变异函数 :return: 左侧矩阵A, 右侧向量b n len(points) A np.zeros((n1, n1)) # 填充变异函数部分 dist_matrix cdist(points, points) A[:n, :n] variogram_func(dist_matrix) # 添加无偏约束条件 A[:n, n] 1 A[n, :n] 1 # 最后一格设为0拉格朗日乘数 A[n, n] 0 return A3.2 预测未知点实现预测函数def kriging_predict(train_points, train_values, test_point, variogram_func): 在单个测试点进行克里金预测 :param train_points: 训练点坐标 (n,2) :param train_values: 训练值 (n,) :param test_point: 测试点坐标 (2,) :param variogram_func: 变异函数 :return: 预测值, 预测方差 n len(train_points) A build_kriging_matrix(train_points, variogram_func) # 构建右侧向量b b np.zeros(n1) dists cdist([test_point], train_points).ravel() b[:n] variogram_func(dists) b[n] 1 # 无偏约束 # 求解权重 weights np.linalg.solve(A, b) # 计算预测值和方差 prediction np.dot(weights[:n], train_values) variance np.dot(weights, b) return prediction, variance4. 网格预测与可视化现在我们可以对整个区域进行预测了。4.1 生成预测网格创建覆盖整个区域的网格点def create_prediction_grid(points, resolution50): 创建覆盖观测点的预测网格 x_min, y_min np.min(points, axis0) x_max, y_max np.max(points, axis0) # 扩展10%的边界 x_pad (x_max - x_min) * 0.1 y_pad (y_max - y_min) * 0.1 x np.linspace(x_min - x_pad, x_max x_pad, resolution) y np.linspace(y_min - y_pad, y_max y_pad, resolution) xx, yy np.meshgrid(x, y) grid_points np.column_stack([xx.ravel(), yy.ravel()]) return xx, yy, grid_points4.2 执行网格预测定义包装好的变异函数def variogram(h): return exponential_variogram(h, nugget, sill, range_) xx, yy, grid_points create_prediction_grid(obs_points) predictions np.zeros(len(grid_points)) variances np.zeros(len(grid_points)) for i, point in enumerate(grid_points): pred, var kriging_predict(obs_points, obs_values, point, variogram) predictions[i] pred variances[i] var4.3 结果可视化绘制预测表面和方差图fig, (ax1, ax2) plt.subplots(1, 2, figsize(14, 6)) # 预测值图 contour ax1.contourf(xx, yy, predictions.reshape(xx.shape), levels20, cmapviridis) ax1.scatter(obs_points[:,0], obs_points[:,1], cred, s50, label观测点) ax1.set_title(克里金预测表面) plt.colorbar(contour, axax1) # 方差图 var_contour ax2.contourf(xx, yy, variances.reshape(xx.shape), levels20, cmapplasma) ax2.scatter(obs_points[:,0], obs_points[:,1], cwhite, s50) ax2.set_title(预测方差) plt.colorbar(var_contour, axax2) plt.tight_layout() plt.show()5. 性能优化与实用技巧当数据量增大时基础实现会遇到性能瓶颈。以下是几个关键优化点5.1 矩阵求解优化对于大型数据集直接求解线性系统效率低下。可以利用from scipy.sparse.linalg import spsolve from scipy.sparse import csr_matrix # 将矩阵转换为稀疏格式 A_sparse csr_matrix(A) weights spsolve(A_sparse, b)5.2 局部邻域选择不必使用全部观测点只需选择预测点周围的邻居def select_neighbors(pred_point, train_points, radius): 选择预测点半径范围内的邻居点 distances cdist([pred_point], train_points).ravel() mask distances radius return train_points[mask], distances[mask]5.3 并行计算利用多核加速网格预测from joblib import Parallel, delayed def parallel_predict(grid_points): results Parallel(n_jobs-1)( delayed(kriging_predict)(obs_points, obs_values, point, variogram) for point in grid_points ) return np.array(results) predictions, variances zip(*parallel_predict(grid_points))6. 不同变异函数模型对比除了指数模型还有其他常用变异函数形式6.1 高斯模型def gaussian_variogram(h, nugget, sill, range_): return nugget (sill - nugget) * (1 - np.exp(-3 * (h/range_)**2))6.2 球型模型def spherical_variogram(h, nugget, sill, range_): condition h range_ return np.where( condition, nugget (sill - nugget) * (1.5*h/range_ - 0.5*(h/range_)**3), nugget (sill - nugget) )6.3 模型选择建议模型类型适用场景平滑性指数型中等连续性中等高斯型高度连续非常平滑球型明确变程中等选择模型时可以计算不同模型的均方根误差(RMSE)from sklearn.model_selection import KFold def cross_validate(points, values, variogram_func, k5): k折交叉验证评估变异函数模型 kf KFold(n_splitsk) errors [] for train_idx, test_idx in kf.split(points): train_pts, test_pts points[train_idx], points[test_idx] train_vals, test_vals values[train_idx], values[test_idx] preds [] for pt in test_pts: pred, _ kriging_predict(train_pts, train_vals, pt, variogram_func) preds.append(pred) rmse np.sqrt(np.mean((np.array(preds) - test_vals)**2)) errors.append(rmse) return np.mean(errors)
别再死记硬背公式了！用Python从零实现一个简易克里金(Kriging)插值模型
用Python从零构建克里金插值模型告别公式恐惧的实战指南当第一次接触克里金插值时大多数人都被那堆晦涩的协方差矩阵和最大似然估计吓退。但今天我们要用Python和NumPy这把瑞士军刀把看似高深的理论拆解成可运行的代码块——你会发现真正理解克里金的最好方式就是亲手实现它。1. 环境准备与数据生成工欲善其事必先利其器。我们从一个干净的Python环境开始import numpy as np import matplotlib.pyplot as plt from scipy.spatial.distance import cdist np.random.seed(42) # 固定随机种子确保可复现1.1 生成模拟观测数据克里金插值需要已知点数据作为基础。我们创建包含20个随机观测点的数据集# 生成二维空间中的随机观测点 obs_points np.random.rand(20, 2) * 10 # 为每个点生成模拟值这里用正弦函数制造空间相关性 obs_values np.sin(obs_points[:,0]/3) * np.cos(obs_points[:,1]/3) np.random.normal(0, 0.1, 20)用Matplotlib可视化这些观测点plt.scatter(obs_points[:,0], obs_points[:,1], cobs_values, cmapviridis, s100) plt.colorbar(label观测值) plt.title(空间观测点分布) plt.show()2. 变异函数空间相关性的量化工具变异函数(variogram)是克里金的核心它量化了空间两点间的差异随距离变化的规律。2.1 计算经验变异函数定义计算经验变异函数的函数def empirical_variogram(points, values): 计算经验变异函数 :param points: 观测点坐标数组 (n,2) :param values: 观测值数组 (n,) :return: 距离数组, 变异函数值数组 # 计算所有点对间的距离 dist_matrix cdist(points, points) # 计算所有点对的值差异平方 value_diff (values[:, None] - values[None, :])**2 # 提取上三角矩阵避免重复计算 upper_tri np.triu_indices_from(dist_matrix, k1) distances dist_matrix[upper_tri] gamma 0.5 * value_diff[upper_tri] return distances, gamma2.2 拟合理论变异函数模型常见理论模型有指数型、高斯型和球型。我们实现指数模型def exponential_variogram(h, nugget, sill, range_): 指数变异函数模型 return nugget (sill - nugget) * (1 - np.exp(-3 * h / range_))通过非线性最小二乘法拟合参数from scipy.optimize import curve_fit # 计算经验值 distances, gamma empirical_variogram(obs_points, obs_values) # 拟合理论模型 popt, _ curve_fit(exponential_variogram, distances, gamma, p0[0.1, 1.0, 5.0], # 初始猜测块金值、基台值、变程 bounds(0, [10, 10, 20])) nugget, sill, range_ popt print(f拟合参数块金值{nugget:.3f}, 基台值{sill:.3f}, 变程{range_:.3f})3. 构建克里金矩阵系统克里金插值的本质是求解一个权重系统使预测值最优。3.1 构建克里金矩阵实现普通克里金(Ordinary Kriging)的矩阵构建def build_kriging_matrix(points, variogram_func): 构建克里金矩阵系统 :param points: 观测点坐标 (n,2) :param variogram_func: 变异函数 :return: 左侧矩阵A, 右侧向量b n len(points) A np.zeros((n1, n1)) # 填充变异函数部分 dist_matrix cdist(points, points) A[:n, :n] variogram_func(dist_matrix) # 添加无偏约束条件 A[:n, n] 1 A[n, :n] 1 # 最后一格设为0拉格朗日乘数 A[n, n] 0 return A3.2 预测未知点实现预测函数def kriging_predict(train_points, train_values, test_point, variogram_func): 在单个测试点进行克里金预测 :param train_points: 训练点坐标 (n,2) :param train_values: 训练值 (n,) :param test_point: 测试点坐标 (2,) :param variogram_func: 变异函数 :return: 预测值, 预测方差 n len(train_points) A build_kriging_matrix(train_points, variogram_func) # 构建右侧向量b b np.zeros(n1) dists cdist([test_point], train_points).ravel() b[:n] variogram_func(dists) b[n] 1 # 无偏约束 # 求解权重 weights np.linalg.solve(A, b) # 计算预测值和方差 prediction np.dot(weights[:n], train_values) variance np.dot(weights, b) return prediction, variance4. 网格预测与可视化现在我们可以对整个区域进行预测了。4.1 生成预测网格创建覆盖整个区域的网格点def create_prediction_grid(points, resolution50): 创建覆盖观测点的预测网格 x_min, y_min np.min(points, axis0) x_max, y_max np.max(points, axis0) # 扩展10%的边界 x_pad (x_max - x_min) * 0.1 y_pad (y_max - y_min) * 0.1 x np.linspace(x_min - x_pad, x_max x_pad, resolution) y np.linspace(y_min - y_pad, y_max y_pad, resolution) xx, yy np.meshgrid(x, y) grid_points np.column_stack([xx.ravel(), yy.ravel()]) return xx, yy, grid_points4.2 执行网格预测定义包装好的变异函数def variogram(h): return exponential_variogram(h, nugget, sill, range_) xx, yy, grid_points create_prediction_grid(obs_points) predictions np.zeros(len(grid_points)) variances np.zeros(len(grid_points)) for i, point in enumerate(grid_points): pred, var kriging_predict(obs_points, obs_values, point, variogram) predictions[i] pred variances[i] var4.3 结果可视化绘制预测表面和方差图fig, (ax1, ax2) plt.subplots(1, 2, figsize(14, 6)) # 预测值图 contour ax1.contourf(xx, yy, predictions.reshape(xx.shape), levels20, cmapviridis) ax1.scatter(obs_points[:,0], obs_points[:,1], cred, s50, label观测点) ax1.set_title(克里金预测表面) plt.colorbar(contour, axax1) # 方差图 var_contour ax2.contourf(xx, yy, variances.reshape(xx.shape), levels20, cmapplasma) ax2.scatter(obs_points[:,0], obs_points[:,1], cwhite, s50) ax2.set_title(预测方差) plt.colorbar(var_contour, axax2) plt.tight_layout() plt.show()5. 性能优化与实用技巧当数据量增大时基础实现会遇到性能瓶颈。以下是几个关键优化点5.1 矩阵求解优化对于大型数据集直接求解线性系统效率低下。可以利用from scipy.sparse.linalg import spsolve from scipy.sparse import csr_matrix # 将矩阵转换为稀疏格式 A_sparse csr_matrix(A) weights spsolve(A_sparse, b)5.2 局部邻域选择不必使用全部观测点只需选择预测点周围的邻居def select_neighbors(pred_point, train_points, radius): 选择预测点半径范围内的邻居点 distances cdist([pred_point], train_points).ravel() mask distances radius return train_points[mask], distances[mask]5.3 并行计算利用多核加速网格预测from joblib import Parallel, delayed def parallel_predict(grid_points): results Parallel(n_jobs-1)( delayed(kriging_predict)(obs_points, obs_values, point, variogram) for point in grid_points ) return np.array(results) predictions, variances zip(*parallel_predict(grid_points))6. 不同变异函数模型对比除了指数模型还有其他常用变异函数形式6.1 高斯模型def gaussian_variogram(h, nugget, sill, range_): return nugget (sill - nugget) * (1 - np.exp(-3 * (h/range_)**2))6.2 球型模型def spherical_variogram(h, nugget, sill, range_): condition h range_ return np.where( condition, nugget (sill - nugget) * (1.5*h/range_ - 0.5*(h/range_)**3), nugget (sill - nugget) )6.3 模型选择建议模型类型适用场景平滑性指数型中等连续性中等高斯型高度连续非常平滑球型明确变程中等选择模型时可以计算不同模型的均方根误差(RMSE)from sklearn.model_selection import KFold def cross_validate(points, values, variogram_func, k5): k折交叉验证评估变异函数模型 kf KFold(n_splitsk) errors [] for train_idx, test_idx in kf.split(points): train_pts, test_pts points[train_idx], points[test_idx] train_vals, test_vals values[train_idx], values[test_idx] preds [] for pt in test_pts: pred, _ kriging_predict(train_pts, train_vals, pt, variogram_func) preds.append(pred) rmse np.sqrt(np.mean((np.array(preds) - test_vals)**2)) errors.append(rmse) return np.mean(errors)
