克里金插值精度问题及变异函数参数差异技术咨询

问题背景与现状 #

我正在对800x500mm墙体结构的法向应变结果做插值，构件中有约150个离散已知采样点，应变存在高、中、低三个明显区域，沿Y方向渐变。

当前遇到的问题：

• 使用Pykrige的泛克里金（Universal Kriging）+高斯变异函数插值，与有限元结果对比的最低误差达37%，数值过高。
• 尝试自定义变异函数模型，得到的参数为：有效变程839.42，基台值、块金效应均为0，插值误差更高。
• 发现Pykrige默认输出的变异函数参数和我用skgstat分析得到的参数差异极大：
  • Pykrige默认：基台值4.94245937e-06，变程1.85944481e+02，块金值1.04982139e-06
  • skgstat自定义：[839.4157043511434, 1.234774627684911e-05, 0]

疑问点：

1. 自定义变异函数模型是否有误？
2. 克里金法是否适用于该场景，还是应该用更简单的插值方法？
3. 是否需要添加趋势模型、用分区克里金，或者引入漂移项？
4. 两种变异函数参数差异的原因是什么？

现有Python代码 #

import numpy as np
import pandas as pd
from pykrige.uk import UniversalKriging
import matplotlib.pyplot as plt
import skgstat as skg

# Load the known strain data and x, y locations from the Excel file using pandas
excel_file = 'Sensor-Results(Exx).xlsx'
data = pd.read_excel(excel_file)

# Extract x, y, and strain values from the DataFrame by using their names on top of the columns
x_known = data['X'].values
y_known = data['Y'].values
strain_known = data['Exx'].values

# Define the grid dimensions
grid_width = 800  # in mm
grid_height = 500  # in mm

# Define the grid division
width1 = 37.5  # in mm
width2 = 50  # in mm
width3 = 282.5  # in mm
width4 = 60  # in mm

# Generate points using the provided divisions (For creating the identical mesh to FEM model)
x1 = np.linspace(0, width1, 1)  # Grid on 0 location
x2 = np.linspace(width1, width1 + width2, 1)  # Grid on 37.5 location
x3 = np.linspace(width1 + width2, 370, 8)  # Grid on 127.857 to 370
x4 = np.linspace(400, width1 + width2 + width3 + width4, 2)
x5 = np.linspace(width1 + width2 + width3 + width4 + width3 / 7, width1 + width2 + width3 * 2 + width4, 7)
x6 = np.linspace(width1 + width2 * 2 + width3 * 2 + width4, grid_width, 2)

# Combine all x points
x_grid = np.unique(np.concatenate((x1, x2, x3, x4, x5, x6)))
y_grid = np.linspace(0, grid_height, 14)

# Create a meshgrid of the grid points
X, Y = np.meshgrid(x_grid, y_grid)

# Flatten the meshgrid arrays
x_flat = X.flatten()
y_flat = Y.flatten()

# Create empirical variogram
coords = np.vstack((x_known, y_known)).T
strain_known = strain_known.flatten()
V = skg.Variogram(coords, strain_known, n_lags=15, normalize=True, model='gaussian', estimator='dowd')
fig = V.plot(show=False)
plt.show()
print(V)

# Extract variogram parameters
variogram_model_parameters = V.parameters
variogram_model = 'gaussian'
print("Variogram Model Parameters:", variogram_model_parameters)

# Perform Universal Kriging interpolation for the full grid with linear drift
uk = UniversalKriging(
    x_known, y_known, strain_known,
    variogram_model=variogram_model,
    variogram_parameters=variogram_model_parameters,
    drift_terms=['regional_linear']
)
strain_interpolated, _ = uk.execute('grid', x_grid, y_grid)

# Replace the interpolated values at known data points with the actual known values
for x, y, strain in zip(x_known, y_known, strain_known):
    xi = np.abs(x_grid - x).argmin()
    yi = np.abs(y_grid - y).argmin()
    strain_interpolated[yi, xi] = strain

# Create a DataFrame with all grid points and interpolated strain values
interpolated_df = pd.DataFrame({
    'X': x_flat,
    'Y': y_flat,
    'Strain': strain_interpolated.flatten()
})

# Sort the DataFrame to achieve the desired zigzag order
interpolated_sorted = pd.DataFrame(columns=['X', 'Y', 'Strain'])

for x_val in x_grid:
    temp_df = interpolated_df[interpolated_df['X'] == x_val]
    idx = np.where(x_grid == x_val)[0][0]
    if idx % 2 == 0:  # Even index
        interpolated_sorted = pd.concat([interpolated_sorted, temp_df])
    else:  # Odd index
        interpolated_sorted = pd.concat([interpolated_sorted, temp_df.iloc[::-1]])

# Save the interpolated data
interpolated_sorted.to_csv('24_06_21_Kriging-Variogram(Exx).txt', sep='\t', index=False, header=True)

# Plot the heat map
strain_grid = strain_interpolated.reshape(X.shape)
plt.figure(figsize=(10, 6))
plt.contourf(X, Y, strain_grid, cmap='rainbow')
plt.colorbar()

# Scatter plot with connecting lines
plt.scatter(x_known, y_known, color='black', s=10)
plt.plot(x_known, y_known, color='red', linewidth=1, linestyle='-')

plt.xlabel('X (mm)')
plt.ylabel('Y (mm)')
plt.title('Strain Heatmap Exx')
plt.show()

问题解答 #

一、变异函数参数差异的核心原因 #

1. 参数顺序不匹配：Pykrige与skgstat对高斯变异函数的参数存储顺序完全不同：
   • Pykrige的顺序是：[块金值, 基台值, 变程]
   • skgstat的顺序是：[变程, 基台值, 块金值]
     你看到的数值差异很大程度上是因为参数顺序搞反了。
2. 归一化设置干扰：你在skgstat中设置了normalize=True，会将数据归一化后再计算变异函数，导致参数量级和Pykrige用原始数据计算的结果完全不匹配。
3. 拟合逻辑差异：两个工具包的变异函数拟合算法、初始值设置不同，当数据存在全局渐变趋势时，拟合结果的偏差会被放大。

二、自定义变异函数模型的问题 #

你的自定义模型出现基台值和块金值为0，主要是两个因素导致：

• 归一化后的数据方差被压缩，拟合算法无法有效区分基台值（数据总方差）和块金值（随机噪声）；
• 应变沿Y方向的渐变趋势属于全局趋势，会干扰变异函数对局部空间相关性的捕捉，导致模型拟合失效。

三、克里金法的适用性与优化方向 #

克里金法完全适用于你的场景，但需要针对性调整：

1. 移除归一化设置：在skgstat创建变异函数时去掉normalize=True，用原始数据计算，保证参数量级与Pykrige匹配。
2. 优化漂移项/趋势模型：
   • 尝试drift_terms=['quadratic']或自定义多项式漂移，适配Y方向的渐变趋势；
   • 先对数据做趋势拟合（比如用线性回归移除Y方向趋势），再用普通克里金插值残差，最后将趋势加回插值结果。
3. 采用分区克里金：将数据按高、中、低应变区域拆分，分别做克里金插值后拼接结果，每个区域内的空间相关性更强，误差会显著降低。
4. 调整变异函数拟合参数：
   • 尝试球形、指数等其他变异函数模型；
   • 减少n_lags至8-10（避免滞后太多导致数据稀疏）；
   • 换用默认的estimator='matheron'代替'dowd'，提升拟合稳定性。

四、替代插值方法建议 #

如果克里金调整后效果仍不理想，可以尝试：

• 径向基函数插值（RBF）：使用scipy的Rbf类，选择'multiquadric'或'gaussian'核函数，适配渐变趋势数据；
• 反距离权重插值（IDW）：实现简单，对局部区域的拟合效果较好，适合有明显分区特征的数据；
• 双样条插值：使用scipy的bisplrep和bisplev，适合平滑连续的应变数据。

内容的提问来源于stack exchange，提问作者Shalva Esakia