别再死记硬背公式了！用Python从零复现Kriging模型（附完整代码与可视化）

用Python实战Kriging模型：从理论盲区到代码落地

第一次接触Kriging模型时，我被那些复杂的数学符号和协方差矩阵吓得不轻。直到有一天，我决定用Python从头实现它——不是调用现成的库，而是真正理解每一步在代码中如何体现。这个过程彻底改变了我对空间插值技术的认知。下面，我将带你用代码"拆解"这个看似高深的地统计学模型。

1. 环境准备与数据生成

工欲善其事，必先利其器。我们先搭建实验环境：

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy.optimize import minimize from sklearn.model_selection import train_test_split

为了演示，我们创建一个具有空间相关性的模拟数据集。这个函数会生成带有噪声的山形曲面数据：

def generate_terrain_data(n_samples=100, noise=0.1): np.random.seed(42) X = np.random.rand(n_samples, 2) * 4 - 2 y = np.exp(-X[:,0]**2 - X[:,1]**2) + np.exp(-(X[:,0]-1)**2 - (X[:,1]-1)**2) y += noise * np.random.randn(n_samples) return X, y X_full, y_full = generate_terrain_data() X_train, X_test, y_train, y_test = train_test_split(X_full, y_full, test_size=0.2)

用3D散点图可视化我们的训练数据：

fig = plt.figure(figsize=(10, 7)) ax = fig.add_subplot(111, projection='3d') ax.scatter(X_train[:,0], X_train[:,1], y_train, c=y_train, cmap='viridis') ax.set_xlabel('X1') ax.set_ylabel('X2') ax.set_zlabel('Y') plt.title('Training Data Spatial Distribution') plt.show()

2. 核心组件实现

2.1 核函数的选择与实现

Kriging的核心在于如何量化空间相关性。我们实现三种常见核函数：

def exponential_kernel(h, theta=1.0): return np.exp(-theta * h) def gaussian_kernel(h, theta=1.0): return np.exp(-(theta * h)**2) def matern_kernel(h, theta=1.0, nu=1.5): h = np.maximum(h, 1e-8) # 避免除以零 if nu == 0.5: return np.exp(-theta * h) elif nu == 1.5: return (1 + np.sqrt(3)*theta*h) * np.exp(-np.sqrt(3)*theta*h) elif nu == 2.5: return (1 + np.sqrt(5)*theta*h + (5/3)*(theta*h)**2) * np.exp(-np.sqrt(5)*theta*h)

可视化不同核函数的形状：

h_vals = np.linspace(0, 3, 100) plt.figure(figsize=(10, 6)) plt.plot(h_vals, exponential_kernel(h_vals), label='Exponential') plt.plot(h_vals, gaussian_kernel(h_vals), label='Gaussian') plt.plot(h_vals, matern_kernel(h_vals, nu=1.5), label='Matern 1.5') plt.xlabel('Distance (h)') plt.ylabel('Correlation') plt.title('Kernel Function Comparison') plt.legend() plt.grid(True) plt.show()

2.2 构建相关矩阵

这是Kriging计算中最关键的步骤之一：

def build_correlation_matrix(X, kernel_func, theta): n = X.shape[0] R = np.zeros((n, n)) for i in range(n): for j in range(i, n): h = np.linalg.norm(X[i] - X[j]) R[i,j] = R[j,i] = kernel_func(h, theta) # 添加小量确保矩阵正定 R += np.eye(n) * 1e-8 return R

提示：在实际应用中，当数据点超过1000个时，建议使用KDTree等空间索引结构加速距离计算

3. 模型训练与超参数优化

3.1 似然函数实现

我们通过最大化对数似然来估计超参数：

def neg_log_likelihood(params, X, y): theta = params[0] mu = params[1] sigma2 = params[2] n = len(y) R = build_correlation_matrix(X, gaussian_kernel, theta) try: L = np.linalg.cholesky(R) Linv = np.linalg.inv(L) Rinv = Linv.T @ Linv except np.linalg.LinAlgError: return np.inf y_mu = y - mu term1 = 0.5 * n * np.log(sigma2) term2 = 0.5 * np.log(np.linalg.det(R)) term3 = (y_mu.T @ Rinv @ y_mu) / (2 * sigma2) return term1 + term2 + term3

3.2 参数优化实战

使用SciPy进行优化：

def train_kriging(X, y): # 初始猜测 initial_theta = 1.0 initial_mu = np.mean(y) initial_sigma2 = np.var(y) bounds = [(1e-3, 100), (None, None), (1e-3, None)] result = minimize( neg_log_likelihood, x0=[initial_theta, initial_mu, initial_sigma2], args=(X, y), bounds=bounds, method='L-BFGS-B' ) return { 'theta': result.x[0], 'mu': result.x[1], 'sigma2': result.x[2], 'X': X, 'y': y }

训练模型并检查结果：

model = train_kriging(X_train, y_train) print(f"Optimized parameters: theta={model['theta']:.3f}, mu={model['mu']:.3f}, sigma2={model['sigma2']:.3f}")

4. 预测与可视化

4.1 预测函数实现

def predict(x_new, model, kernel_func=gaussian_kernel): X = model['X'] y = model['y'] theta = model['theta'] mu = model['mu'] n = len(y) R = build_correlation_matrix(X, kernel_func, theta) y_mu = y - mu # 计算新点与训练点的相关性向量 r = np.zeros(n) for i in range(n): h = np.linalg.norm(x_new - X[i]) r[i] = kernel_func(h, theta) try: L = np.linalg.cholesky(R) Linv = np.linalg.inv(L) Rinv = Linv.T @ Linv except np.linalg.LinAlgError: return mu, 0 # 预测均值 y_pred = mu + r.T @ Rinv @ y_mu # 预测方差 mse = model['sigma2'] * (1 - r.T @ Rinv @ r) return y_pred, mse

4.2 空间预测可视化

生成网格预测：

def plot_prediction_surface(model): x1 = np.linspace(-2, 2, 50) x2 = np.linspace(-2, 2, 50) X1, X2 = np.meshgrid(x1, x2) Z = np.zeros_like(X1) for i in range(X1.shape[0]): for j in range(X1.shape[1]): Z[i,j], _ = predict(np.array([X1[i,j], X2[i,j]]), model) fig = plt.figure(figsize=(14, 6)) ax1 = fig.add_subplot(121, projection='3d') ax1.plot_surface(X1, X2, Z, cmap='viridis', alpha=0.8) ax1.scatter(X_train[:,0], X_train[:,1], y_train, c='r', s=20) ax1.set_title('Kriging Prediction Surface') ax2 = fig.add_subplot(122) contour = ax2.contourf(X1, X2, Z, levels=20, cmap='viridis') ax2.scatter(X_train[:,0], X_train[:,1], c='r', s=20) plt.colorbar(contour) ax2.set_title('Contour Plot with Training Points') plt.tight_layout() plt.show() plot_prediction_surface(model)

4.3 与scikit-learn的对比

from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import RBF sklearn_kernel = RBF(length_scale=1/model['theta']) gp = GaussianProcessRegressor(kernel=sklearn_kernel, alpha=1e-8) gp.fit(X_train, y_train) # 比较预测结果 test_pred_custom = np.array([predict(x, model)[0] for x in X_test]) test_pred_sklearn = gp.predict(X_test) plt.figure(figsize=(10, 6)) plt.scatter(y_test, test_pred_custom, label='Custom Implementation') plt.scatter(y_test, test_pred_sklearn, label='scikit-learn') plt.plot([min(y_test), max(y_test)], [min(y_test), max(y_test)], 'k--') plt.xlabel('True Values') plt.ylabel('Predictions') plt.title('Comparison with scikit-learn Implementation') plt.legend() plt.grid(True) plt.show()

5. 高级话题与性能优化

5.1 处理大规模数据

当数据量增大时，我们需要优化实现：

from scipy.spatial import KDTree def fast_build_correlation_matrix(X, kernel_func, theta, leaf_size=40): tree = KDTree(X, leafsize=leaf_size) n = X.shape[0] R = np.eye(n) # 对角线初始化为1 for i in range(n): # 只计算与邻近点的相关性 indices = tree.query_ball_point(X[i], r=3/theta) for j in indices: if j > i: # 避免重复计算 h = np.linalg.norm(X[i] - X[j]) R[i,j] = R[j,i] = kernel_func(h, theta) return R

5.2 多维度参数优化

对于高维输入空间，我们可以为每个维度设置不同的θ：

def anisotropic_kernel(X1, X2, thetas): sq_dist = 0 for i in range(len(thetas)): sq_dist += (thetas[i] * (X1[i] - X2[i]))**2 return np.exp(-sq_dist)

5.3 交叉验证实现

def cross_validate_kriging(X, y, n_folds=5, theta_range=(0.1, 10)): fold_size = len(X) // n_folds indices = np.random.permutation(len(X)) results = [] for i in range(n_folds): test_indices = indices[i*fold_size : (i+1)*fold_size] train_indices = np.setdiff1d(indices, test_indices) X_train, X_test = X[train_indices], X[test_indices] y_train, y_test = y[train_indices], y[test_indices] model = train_kriging(X_train, y_train) test_pred = np.array([predict(x, model)[0] for x in X_test]) mse = np.mean((test_pred - y_test)**2) results.append(mse) return np.mean(results)