% 克里金插值示例

clc;

clear;

% 生成模拟数据

x = linspace(0, 10, 11);

y = linspace(0, 10, 11);

[X, Y] = meshgrid(x, y);

Z = sin(sqrt(X.^2 + Y.^2)) + 0.1 * randn(size(X));

% 设置克里金参数

nugget = 0.1; % 块金值

range = 1; % 范围

sill = 1; % 基台值

azimuth = 0; % 方位角

% 计算距离矩阵

D = pdist2(cat(3, X(:), Y(:)), cat(3, X(:), Y(:)));

D = D / range;

% 计算半方差矩阵

V = nugget + (sill - nugget) * exp(-D);

% 计算权重矩阵

W = V .* ones(size(D));

for i = 1:size(D, 1)

    W(i, i) = sum(W(i, :)) - W(i, i);

end

% 克里金插值

[~, I] = max(W);

Z_pred = Z(I) + (Z(I) - Z(I(I))) * (1 - exp(-D(I)));

% 绘制结果

figure;

surf(X, Y, Z, 'FaceColor', 'b');

hold on;

surf(X, Y, Z_pred, 'FaceColor', 'r');

title('克里金插值示例');

xlabel('X');

ylabel('Y');

zlabel('Z');

legend('原始数据', '插值结果');