克里金插值 調用matlab工具箱

克里金插值

克里金插值是依據協方差函數對隨機過程或隨機場進行空間建模和插值的回歸算法。

克里金插值法的公式為：

 

式中為待插入的各點的重金屬污染值，為已知點的重金屬污染值，為每個點的權重值。

用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