最近在研究克里金插值,拜读了@lanainluv的笔记,备受启发, 在这里做一些补充,并分享自己的代码,希望对各位有所帮助,有误的地方请批评指正!(有帮助的话点个赞吧~)
@lanainluv的原文链接(推导过程):(54条消息) (超全面,超基础)Kriging插值推导理论笔记,算法,普通克里金_kridge插值原理_lanainluv的博客-CSDN博客
在编程时,我遇到了一些问题,解决方案如下:
1)下图方程组左边的矩阵不可逆
解决方案:将任意两点的距离代入半变异函数中重新计算矩阵中的半方差,如下所示:
rij=Guassian_func_3(dij);
2)d-r散点图十分混乱,没有规律,拟合困难
解决方案:用k-means聚类对距离进行聚类,聚到一类的距离和半方差均取平均值,再进行拟合,如下图所示。
这里采用的是matlab的拟合工具箱cftool(x,y),使用Gaussian函数进行拟合,也可以用Fourier和多项式函数。
算例结果如下所示:
案例1:
案例2:
这里也是使用Gaussian函数(3项)进行拟合,用sin函数是最好,这里是为了说明算法的有效性。
代码如下:
案例1 主函数:
clear all
close all
clc
S=[0.700000000000000,59.6000000000000;2.10000000000000,82.7000000000000;4.70000000000000,75.1000000000000;4.80000000000000,52.8000000000000;5.90000000000000,67.1000000000000;6,35.7000000000000;6.40000000000000,33.7000000000000;7,46.7000000000000;8.20000000000000,40.1000000000000;13.3000000000000,0.600000000000000;13.3000000000000,68.2000000000000;13.4000000000000,31.3000000000000;17.8000000000000,6.90000000000000;20.1000000000000,66.3000000000000;22.7000000000000,87.6000000000000;23,93.9000000000000;24.3000000000000,73;24.8000000000000,15.1000000000000;24.8000000000000,26.3000000000000;26.4000000000000,58;26.9000000000000,65;27.7000000000000,83.3000000000000;27.9000000000000,90.8000000000000;29.1000000000000,47.9000000000000;29.5000000000000,89.4000000000000;30.1000000000000,6.10000000000000;30.8000000000000,12.1000000000000;32.7000000000000,40.2000000000000;34.8000000000000,8.10000000000000;35.3000000000000,32;37,70.3000000000000;38.2000000000000,77.9000000000000;38.9000000000000,23.3000000000000;39.4000000000000,82.5000000000000;43,4.70000000000000;43.7000000000000,7.60000000000000;46.4000000000000,84.1000000000000;46.7000000000000,10.6000000000000;49.9000000000000,22.1000000000000;51,88.8000000000000;52.8000000000000,68.9000000000000;52.9000000000000,32.7000000000000;55.5000000000000,92.9000000000000;56,1.60000000000000;60.6000000000000,75.2000000000000;62.1000000000000,26.6000000000000;63,12.7000000000000;69,75.6000000000000;70.5000000000000,83.7000000000000;70.9000000000000,11;71.5000000000000,29.5000000000000;78.1000000000000,45.5000000000000;78.2000000000000,9.10000000000000;78.4000000000000,20;80.5000000000000,55.9000000000000;81.1000000000000,51;83.8000000000000,7.90000000000000;84.5000000000000,11;85.2000000000000,67.3000000000000;85.5000000000000,73;86.7000000000000,70.4000000000000;87.2000000000000,55.7000000000000;88.1000000000000,0;88.4000000000000,12.1000000000000;88.4000000000000,99.6000000000000;88.8000000000000,82.9000000000000;88.9000000000000,6.20000000000000;90.6000000000000,7;90.7000000000000,49.6000000000000;91.5000000000000,55.4000000000000;92.9000000000000,46.8000000000000;93.4000000000000,70.9000000000000;94.8000000000000,71.5000000000000;96.2000000000000,84.3000000000000;98.2000000000000,58.2000000000000];
Y=[34.1000000000000;42.2000000000000;39.5000000000000;34.3000000000000;37;35.9000000000000;36.4000000000000;34.6000000000000;35.4000000000000;44.7000000000000;37.8000000000000;37.8000000000000;43.9000000000000;37.7000000000000;42.8000000000000;43.6000000000000;39.3000000000000;42.3000000000000;39.7000000000000;36.9000000000000;37.8000000000000;41.8000000000000;43.3000000000000;36.7000000000000;43;43.6000000000000;42.8000000000000;37.5000000000000;43.3000000000000;38.8000000000000;39.2000000000000;40.7000000000000;40.5000000000000;41.4000000000000;43.3000000000000;43.1000000000000;41.5000000000000;42.6000000000000;40.7000000000000;42;39.3000000000000;39.2000000000000;42.2000000000000;42.7000000000000;40.1000000000000;40.1000000000000;41.8000000000000;40.1000000000000;40.9000000000000;41.7000000000000;40.8000000000000;38.7000000000000;41.7000000000000;40.8000000000000;38.7000000000000;38.6000000000000;41.6000000000000;41.5000000000000;39.4000000000000;39.8000000000000;39.6000000000000;38.8000000000000;41.6000000000000;41.3000000000000;41.2000000000000;40.5000000000000;41.5000000000000;41.5000000000000;38.9000000000000;39;39.1000000000000;39.7000000000000;39.7000000000000;40.3000000000000;39.5000000000000];
n=size(S,1);
rij=zeros(n,n);
dij=zeros(n,n);
d_r=[];
for i=1:n
for j=1:n
rij(i,j)=0.5*(Y(i,1)-Y(j,1)).^2;
dij(i,j)=sqrt(sum((S(i,:)-S(j,:)).^2));%欧式距离
% dij(i,j)=sum(abs(S(i,:)-S(j,:)));%曼哈顿距离(在优化问题中易于线性化)
d_r=[d_r;[dij(i,j),rij(i,j)]];
end
end
K=20;
[idx,~] =kmeans(d_r(:,1),K);
d_r_divided=[];
for i=1:K
d_r_divided=[d_r_divided;mean(d_r(find(idx==i),:),1)];
end
d_r_divided=sortrows(d_r_divided,1);
% d_r_divided=d_r;
figure(1)
scatter(d_r_divided(:,1),d_r_divided(:,2),'LineWidth',2)
hold on
nihe_d_r_divided(:,1)=[min(d_r_divided(:,1)):0.1:1.01*max(d_r_divided(:,1))];
nihe_d_r_divided(:,2)=Guassian_func_3(nihe_d_r_divided(:,1));
plot(nihe_d_r_divided(:,1),nihe_d_r_divided(:,2),'LineWidth',2)
gca = legend('散点','半变异函数','Location','east');
set(gca,'FontSize',20);
% cftool(d_r_divided(:,1),d_r_divided(:,2))
%预测
rij=Guassian_func_3(dij);
rij(logical(eye(size(rij))))=0;
X = gridsamp([0 0;100 100], 40);
[m,~]=size(X);
YX=zeros(m,1);
for i=1:size(X,1)
x=X(i,:);
n=size(S,1);
dix=sqrt(sum((repmat(x,n,1)-S).^2,2));
% dix=sum(abs(repmat(x,n,1)-S),2);
rix=Guassian_func_3(dix);
temp=[[rij,ones(size(rij,1),1)];[ones(1,size(rij,1)),0]];
lambda=temp\[rix;1];
y=sum(lambda(1:n).*Y);
YX(i)=y;
end
X1 = reshape(X(:,1),40,40); X2 = reshape(X(:,2),40,40);
YX = reshape(YX, size(X1));
figure(2), mesh(X1, X2, YX)%绘制预测表面
hold on
plot3(S(:,1),S(:,2),Y,'.k', 'MarkerSize',10)%绘制原始散点数据
案例2 主函数:
clear all
close all
clc
S=[1:0.1:10]';
Y=sin(S);
n=size(S,1);
rij=zeros(n,n);
dij=zeros(n,n);
d_r=[];
for i=1:n
for j=1:n
rij(i,j)=0.5*(Y(i,1)-Y(j,1)).^2;
dij(i,j)=sqrt(sum((S(i,:)-S(j,:)).^2));%欧式距离
d_r=[d_r;[dij(i,j),rij(i,j)]];
end
end
K=20;
[idx,~] =kmeans(d_r(:,1),K);
d_r_divided=[];
for i=1:K
d_r_divided=[d_r_divided;mean(d_r(find(idx==i),:),1)];
end
d_r_divided=sortrows(d_r_divided,1);
% d_r_divided=d_r;
figure(1)
scatter(d_r_divided(:,1),d_r_divided(:,2),'LineWidth',2)
hold on
nihe_d_r_divided(:,1)=[min(d_r_divided(:,1)):0.1:1.01*max(d_r_divided(:,1))];
nihe_d_r_divided(:,2)=Guassian_func_3_test(nihe_d_r_divided(:,1));
plot(nihe_d_r_divided(:,1),nihe_d_r_divided(:,2),'LineWidth',2)
gca = legend('散点','半变异函数','Location','east');
set(gca,'FontSize',20);
% cftool(d_r_divided(:,1),d_r_divided(:,2))
%预测
rij=Guassian_func_3_test(dij);
rij(logical(eye(size(rij))))=0;
X = linspace(0,10,50)';
[m,~]=size(X);
YX=zeros(m,2);
for i=1:size(X,1)
x=X(i,:);
n=size(S,1);
dix=sqrt(sum((repmat(x,n,1)-S).^2,2));
rix=Guassian_func_3_test(dix);
temp=[[rij,ones(size(rij,1),1)];[ones(1,size(rij,1)),0]];
lambda=temp\[rix;1];
y=sum(lambda(1:n).*Y);
YX(i,:)=[y,sin(x)];
end
figure(2)
plot(X,YX(:,1),'r*')
hold on
plot(X,YX(:,2),'bo')
案例1 半变异函数(高斯函数):
function [y] = Guassian_func_3(x)
a1=7.112; b1=99.03; c1=17.06;
a2=6.35; b2=39.31; c2=19.68;
a3=5.351; b3=65.6; c3=18.27;
[m,n]=size(x);
y=zeros(m,n);
for i=1:m
for j=1:n
y(i,j)=a1*exp(-((x(i,j)-b1)/c1)^2) + a2*exp(-((x(i,j)-b2)/c2)^2) + a3*exp(-((x(i,j)-b3)/c3)^2);
end
end
end
案例2 半变异函数(高斯函数):
function [y] = Guassian_func_3_test(x)
a1 = 0.7279;
b1 = 8.405;
c1 = 1.018;
a2 = 0.8317;
b2 = 3.668;
c2 = 1.412;
a3 = 0.5205;
b3 = 2.167;
c3 = 1.183;
[m,n]=size(x);
y=zeros(m,n);
for i=1:m
for j=1:n
y(i,j)=a1*exp(-((x(i,j)-b1)/c1)^2) + a2*exp(-((x(i,j)-b2)/c2)^2) + a3*exp(-((x(i,j)-b3)/c3)^2);
end
end
end
function S = gridsamp(range, q)
%GRIDSAMP n-dimensional grid over given range
%
% Call: S = gridsamp(range, q)
%
% range : 2*n matrix with lower and upper limits
% q : n-vector, q(j) is the number of points
% in the j'th direction.
% If q is a scalar, then all q(j) = q
% S : m*n array with points, m = prod(q)
% hbn@imm.dtu.dk
% Last update June 25, 2002
[mr n] = size(range); dr = diff(range);
if mr ~= 2 | any(dr < 0)
error('range must be an array with two rows and range(1,:) <= range(2,:)')
end
sq = size(q);
if min(sq) > 1 | any(q <= 0)
error('q must be a vector with non-negative elements')
end
p = length(q);
if p == 1, q = repmat(q,1,n);
elseif p ~= n
error(sprintf('length of q must be either 1 or %d',n))
end
% Check for degenerate intervals
i = find(dr == 0);
if ~isempty(i), q(i) = 0*q(i); end
% Recursive computation
if n > 1
A = gridsamp(range(:,2:end), q(2:end)); % Recursive call
[m p] = size(A); q = q(1);
S = [zeros(m*q,1) repmat(A,q,1)];
y = linspace(range(1,1),range(2,1), q);
k = 1:m;
for i = 1 : q
S(k,1) = repmat(y(i),m,1); k = k + m;
end
else
S = linspace(range(1,1),range(2,1), q).';
end