基于应变场的泛克里金插值精度问题及优化咨询

问题背景 #

对800×500mm的墙体结构进行法向应变插值，构件含约150个离散测点。使用Pykrige的泛克里金（Universal Kriging）+高斯变异函数插值，与有限元结果对比最低误差达37%；尝试自定义变异函数模型时出现异常：有效范围839.42，基台值（sill）和块金效应（nugget）均为0，误差进一步升高。

应变数据来自自建有限元模型，等值线图显示应变在y方向呈渐变的高、中、低分区。

核心疑问 #

1. 变异函数模型构建是否存在错误？
2. 泛克里金是否适用于该场景，还是应采用更简单的插值方法？
3. 是否需要为插值添加趋势模型、采用分区克里金（因应变存在高、中、低三个明显分区）或引入漂移项（drift_term）？
4. 如何判断数据适配何种变异函数？（滞后距分组值增长后趋于平稳，对此存疑）

参数差异问题 #

默认泛克里金输出的变异函数参数：

• 基台值：4.94245937e-06
• 范围：1.85944481e+02
• 块金效应：1.04982139e-06

自定义高斯变异函数参数（来自skgstat）：[839.4157043511434, 1.234774627684911e-05, 0]
两者差异显著，请问原因是什么？

附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. 变异函数模型构建是否有误？

自定义变异函数出现基台值和块金效应为0，核心原因是**normalize=True**参数导致数据被缩放到[0,1]范围，掩盖了真实变异程度，使得拟合参数异常。建议去掉该参数，用原始数据计算变异函数。另外，estimator='dowd'对异常值鲁棒性强，若数据无明显异常，改用默认的'matheron'估计量更贴合应变数据的常规变异特征。

2. 泛克里金是否适用，还是用简单插值？

泛克里金适合存在明显空间趋势的场景（你的数据在y方向有渐变分区），但当前误差过高时，可先对比简单插值方法做基准验证：

• 反距离加权（IDW）：适合趋势明显、测点分布均匀的场景，实现简单，能快速得到基准误差。
• 径向基函数（RBF）：可较好拟合平滑连续场，适配应变这类连续数据。
  如果简单方法误差仍高，说明问题不在插值方法本身，而是数据预处理或变异函数参数设置不当。

3. 趋势模型、分区克里金或漂移项的选择？

• 漂移项：当前用的regional_linear线性漂移可能无法匹配y方向的渐变分区，可尝试'quadratic'二次漂移，或自定义多项式漂移（Pykrige支持自定义漂移函数）。
• 分区克里金：既然应变有明确的高、中、低y向分区，直接按y坐标拆分数据为三个子集，分别做克里金插值再合并结果，能大幅降低跨分区的变异干扰，提升精度。
• 趋势模型验证：先对比普通克里金（不添加漂移）与泛克里金的结果，若普通克里金误差更低，说明当前漂移项选择不当，反而引入了偏差。

4. 如何判断适配的变异函数？

变异函数选择需匹配实验变异曲线形态：

• 高斯模型：曲线平滑趋近基台值，适合连续渐变的场（如你的应变数据）。
• 球状模型：曲线先快速上升后平缓，适合空间相关性随距离快速衰减的场景。
• 指数模型：曲线渐近趋近基台值，相关性衰减速度介于高斯和球状之间。
  你的滞后距分组值趋于平稳，说明数据存在明确空间相关性，且达到一定距离后相关性不再变化，符合高斯/球状模型特征。可通过skgstat的V.plot()查看实验曲线，优先选择拟合R²值最高的模型（skgstat会输出该指标）。

参数差异的原因

Pykrige默认参数与skgstat参数差异源于三点：

1. 数据处理不同：skgstat中normalize=True缩放了数据，导致拟合参数的量级和范围完全不同。
2. 拟合逻辑差异：Pykrige的自动拟合侧重克里金插值的稳定性，而skgstat仅单纯拟合实验变异曲线。
3. 估计量不同：Pykrige默认用Matheron估计量，你在skgstat中用了Dowd估计量，不同估计量计算的实验变异值不同，直接导致拟合参数偏差。

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