克里金插值调用matlab工具箱

绝地灬酷狼 2021-11-02 01:56 362阅读 0赞

克里金插值

克里金插值是依据协方差函数对随机过程或随机场进行空间建模和插值的回归算法。

克里金插值法的公式为：

式中为待插入的各点的重金属污染值，为已知点的重金属污染值，为每个点的权重值。

用BLUP理论求解克里金权重：

将随机场中变量的估计表示为包含随机误差的线性系统，则BLUP可表示为选择线性系统参数使估计值和真实值方差最小：

式中为未知点，\{为随机场的样本，为权重系数，通常被称为克里金权重。由方差定义可知，当估计值和真实值的数学期望相同时，两者的方差最小

使用上述BLUP条件求解的权重系数包含样本点与未知点间的协方差函数。

克里金法是一种在许多领域都很有用的地址统计格网化方法，很符合本题分析污染物的浓度分布，而且克里金插值产生的结果更自然，能够有效的避免异常值的产生，也能给出标准误差，这得益于克里金插值算法考虑了被估计点的位置与已知点位置的相互之间的关系，也考虑了已知点位置之间的关系。所以，更能客观的反应污染物的分布规律，估值的精度也就相对较高。

代码如下

%采样地形图绘制
    clc,clear
    data=xlsread('zz.xls','附件1','B4:D322');%读取文件
    %获得数据
    S=data(:,1:2);
    Y=data(:,3);
    %采用克里金插值法
        theta = [10 10]; lob = [1e-1 1e-1]; upb = [20 20];%参数
        %调用克里金插值算法工具箱
        %进行拟合操作
        [dmodel, perf] = dacefit(S, Y, @regpoly0, @corrspherical, theta, lob, upb) 
        m=100;
        %插值计算
        X = gridsamp([0 0;30000 20000], m);
        [YX MSE] = predictor(X, dmodel);
        %获得插值
        X1 = reshape(X(:,1),m,m); X2 = reshape(X(:,2),m,m);
        YX = reshape(YX, size(X1));
        %作图
        figure;
        surf(X1, X2, YX)
        hold on,  
        hold off
        xlabel('x/m')
        ylabel('y/m')
        zlabel('海拔')
        title('采样地形图')
        figure;
        contourf(X1,X2,YX)%做平面图
        [C,h] = contour(X1,X2,YX);  
        clabel(C,h)
        xlabel('x/m')
        ylabel('y/m')
        zlabel('海拔')
        title('采样地形图')
    
    %污染物浓度分布
    clc,clear
    b={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn'};
    nd=xlsread('zz.xls','附件2','B4:I322');%读取文件
    S=xlsread('zz.xls','附件1','B4:C322');%读取文件
    %循环读取数据
    for i=1:8
        Y=fix(nd(:,i));
        %采用克里金插值法
        theta = [10 10]; lob = [1e-1 1e-1]; upb = [20 20];%参数
        %调用克里金插值算法工具箱
        %进行拟合操作
        [dmodel, perf] = dacefit(S, Y, @regpoly0, @corrspherical, theta, lob, upb) 
        m=100;
         %插值计算
        X = gridsamp([0 0;30000 20000], m);
        [YX MSE] = predictor(X, dmodel);
        %获得插值
        X1 = reshape(X(:,1),m,m); X2 = reshape(X(:,2),m,m);
        YX = reshape(YX, size(X1));
        %作图
        figure;
        mesh(X1, X2, YX)
        hold on,  
        hold off
        xlabel('x/m')
        ylabel('y/m')
        zlabel('浓度')
        title([b(i)])
        figure;
        contourf(X1,X2,YX)
        %做平面图
        [C,h] = contour(X1,X2,YX);
        xlabel('x/m')
        ylabel('y/m')
        zlabel('浓度')
        title([b(i)])
    end

需要调用工具箱文件

在我的百度云盘里有，希望对你有帮助

链接：https://pan.baidu.com/s/1O-mqKowNBJ06llEldHu90A  
提取码：g4wl

转载于:https://www.cnblogs.com/zq-dmhy/p/11288080.html