Matlab利用球状模型变差函数进行克里金插值求某点的海表盐度

作业记录
题目：设有某一海区，在海面上S1，S2，S3，S4处有四个测量点，其海表盐度分别为30‰，31 ‰，32 ‰ ，33 ‰ 。据此利用克里金方法估计S0点处的海表盐度值Z0

解答：
1.利用球状模型变差函数进行克里金插值思路
利用克里金插值在未有观测值的位点获得海表盐度的值。可以通过4个步骤来获得结果，首先是数据准备，其次是结构分析获得变差函数的关键参数，再次就是求解权重系数，最后获得估计值。
（1）数据准备如题所示
（2）利用结构分析获得变异函数的关键值
设海表盐度S0是二阶平稳的。其在平面上的二维变差函数是一个各向同性的球状模型，如果用h表示距离，则常用的球状模型公式为

（3）求解权系数
通过克里金方程组来求解权系数。首先需要根据变异函数表达式、变异函数和协方差函数的关系式来求取观测点之间以及观测和待估计点之间的协方差矩阵。有了协方差矩阵，即可求解克里金方程组得到权系数。
普通克里金方程组的矩阵形式为


可得
（4）最后一步，加权求和，计算估计值。

2.程序实现
#kriging.m#

clc;clear all
%%设海表盐度S0是二阶平稳的，其在平面上的二维变差函数是一个各向同性的球状模型%%
Z=[30 31 32 33];%S1,S2,S3,S4测量点的已知盐度
a=200;%变程
c=20;%拱高
c0=2;%块金值
z=[50 50;100 0;200 100;0 150];%以左上角顶点为原点，标注S1,S2,S3,S4测量点的位置
h=pdist(z); %pdist函数计算z中各测量点之间的相互距离h
for i=1:6
    C(i)=22-spherical(h(i));%spherical(h(i))是求各距离下的变差函数值，22为基台值；
    %基台值=块金值+拱高；在二阶平稳条件下，变差函数=基台值-协方差
end
%%%%%计算各观测点之间的协方差
C11=c+c0; C22=C11;%C11=C22=C33=C44=基台值=22
C33=C11;  C44=C11;
C12=C(1); C21=C12;%h(1)=sqrt(50*50+50*50)
C13=C(2); C31=C13;%h(2)=sqrt(150*150+50*50)
C14=C(3); C41=C14;%h(3)=sqrt(50*50+100*100)
C23=C(4); C32=C23;%h(4)=sqrt(100*100+100*100)
C24=C(5); C42=C24;%h(5)=sqrt(100*100+150*150)
C34=C(6); C43=C34;%h(6)=sqrt(200*200+50*50)

%%%%%解克里金方程组，求解权系数
K=[C11 C12 C13 C14 1;
    C21 C22 C23 C24 1;
    C31 C32 C33 C34 1;
    C41 C42 C43 C44 1;
    1 1 1 1 0];
%%%%%计算各观测点与待估计点之间的协方差
C01=22-spherical(50);
C02=C14;
C03=22-spherical(150);
C04=C12;

M=[C01;C02;C03;C04;1];

%%%%%加权求和，计算待估计点的盐度值
namta=inv(K)*M;%namta为所求权系数,inv(K)是求K的逆矩阵
Z1=0;%初始化
for i=1:4
Z1(i)=namta(i)*Z(i)
end
Z0=sum(Z1);%%Z0为所求待估计点S0的盐度

#spherical.m#

%%%%%%%%%%%%%%%%%%%%%%球状模型变差函数%%%%%%%%%%%%%%%%%%%%%
%%%假设研究区域满足二阶平稳%%%
function y=spherical(h)
a=200;%变程
c=20;%拱高
c0=2;%块金值

%%%球状模型公式%%%
if h==0  %h表示距离
    y=0;
elseif (h>=0)&(h<=a)
    y=c0+c.*(1.5*h/a-0.5*(h/a).^3);
elseif h>=a
    y=c0+c;
end

end

happy happy study, day day up! over~