本文已参加「新人创造礼」活动，一同敞开创造之路。

  本文介绍基于MATLAB求取空间数据的变异函数，并制作经历半方差图的办法。

  其间，因为本文所用的数据并不是我的，因而遗憾不能将数据同时展示给大家；可是根据本篇博客的思想与对代码的详细解说，大家用自己的数据，能够将空间数据变异函数核算与经历半方差图制作的悉数进程与剖析办法加以完好重现。

1 数据处理

1.1 数据读取

  本文中，我的初始数据为某区域658个土壤采样点的空间方位（X与Y，单位为米）、pH值、有机质含量与全氮含量。这些数据均存储于data.xls文件中；而后期操作多于MATLAB软件中进行。因而，首要需将源数据挑选性地导入MATLAB软件中。

  利用MATLAB软件中xlsread函数能够完成这一功能。详细代码附于本文的1.3 正态散布查验及转化处。

1.2 异常数据除掉

  得到的采样点数据因为采样记载、试验室测试等进程，可能具有必定差错，然后出现个别异常值。选用均匀值加规范差法对这些异常数据加以筛选、除掉。

  别离利用均匀值加规范差法中“2S”与“3S”办法加以处理，发现“2S”办法处理效果相对后者较好，故后续试验取“2S”办法处理成果继续进行。

  其间，“2S”办法是指将数值大于或小于其均匀值2倍规范差的部分视作异常值，“3S”办法则是指将数值大于或小于其均匀值3倍规范差的部分视作异常值。

  得到异常值后，将其从658个采样点中除掉；剩下的采样点数据继续后续操作。

  本部分详细代码附于1.3 正态散布查验及转化处。

1.3 正态散布查验及转化

  核算变异函数需建立在初始数据契合正态散布的假设之上；而采样点数据并不必定契合正态散布。因而，咱们需求对原始数据加以正态散布查验。

  一般地，正态散布查验能够经过数值查验与直方图、QQ图等图画加以直观判断。本文综合采取以上两种数值、图画查验办法，共同判断正态散布特性。

  针对数值查验办法，我在一开始准备挑选选用Kolmogorov-Smirnov查验办法；但因为了解到，这一办法只是适用于规范正态查验，因而随后改用Lilliefors查验。

  Kolmogorov-Smirnov查验经过样本的经历散布函数与给定散布函数的比较，推断该样本是否来自给定散布函数的总体；当其用于正态性查验时只能做规范正态查验。

  Lilliefors查验则将上述Kolmogorov-Smirnov查验改善，其可用于一般的正态散布查验。

  QQ图（Quantile Quantile Plot）是一种散点图，其横坐标表明某一样本数据的分位数，纵坐标则表明另一样本数据的分位数；横坐标与纵坐标组成的散点图代表同一个累计概率所对应的分位数。

  因而，QQ图具有这样的特点：针对y=x这一直线，若散点图中各点均在直线邻近散布，则说明两个样本为同等散布；因而，若将横坐标（纵坐标）表明为一个规范正态散布样本的分位数，则散点图中各点均在上述直线邻近散布能够说明，纵坐标（横坐标）表明的样本契合或根本近似契合正态散布。本文选用将横坐标表明为正态散布的方法。

  此外，PP图（Probability Probability Plot）相同能够用于正态散布的查验。PP图横坐标表明某一样本数据的累积概率，纵坐标则表明另一样本数据的累积概率；其根据变量的累积概率对应于所指定的理论散布累积概率并制作的散点图，用于直观地检测样本数据是否契合某一概率散布。和QQ图类似，假如被查验的数据契合所指定的散布，则其各点均在上述直线邻近散布。若将横坐标（纵坐标）表明为一个规范正态散布样本的分位数，则散点图中各点均在直线邻近散布能够说明，纵坐标（横坐标）表明的样本契合或根本近似契合正态散布。

  三种土壤特点，我挑选首要以pH数值为例进行操作。经过上述数值查验、图画查验办法，查验得到除掉异常值后的原始pH数值数据并不契合正态散布这一结论。因而，测验对原数据加以对数、开平方等转化处理；随后发现，原始pH值开平方数据的正态散布特征虽然仍旧无法经过较为严厉的Lilliefors查验，但其直方图、QQ图的图画查验成果较为挨近正态散布，并较之前二者更加显着。故后续取开平方处理成果继续进行。

  值得一提的是，本文后半部分得到pH值开平方数据的试验变异函数及其散点图后，在对其他两种空间特点数据（即有机质含量与全氮含量）进行相同的操作时，发现全氮含量数据在经过“2S”办法除掉异常值后，其原始形式的数据是能够经过Lilliefors查验的，且其直方图、QQ图散布特点十分挨近正态散布。

  我亦准备测验对空间特点数据进行反正弦转化。但随后发现，已有三种特点数值的原始数据并不严厉散布在-1至1的区间内，因而并未对其进行反正弦方法的转化。

  经过上述查验、转化处理往后的图画查验成果如下所示。

  以上部分代码如下：

clc;clear;
info=xlsread('data.xls');
oPH=info(:,3);
oOM=info(:,4);
oTN=info(:,5);
mPH=mean(oPH);
sPH=std(oPH);
num2=find(oPH>(mPH+2*sPH)|oPH<(mPH-2*sPH));
num3=find(oPH>(mPH+3*sPH)|oPH<(mPH-3*sPH));
PH=oPH;
for i=1:length(num2)
    n=num2(i,1);
    PH(n,:)=[0];
end
PH(all(PH==0,2),:)=[];
%KSTest(PH,0.05)
H1=lillietest(PH);
for i=1:length(PH)
    lPH(i,:)=log(PH(i,:));
end
H2=lillietest(lPH);
for i=1:length(PH)
    sqPH(i,:)=(PH(i,:))^0.5;
end
H3=lillietest(sqPH);
% for i=1:length(PH)
%     arcPH(i,:)=asin(PH(i,:));
% end
% 
% H4=lillietest(arcPH);
subplot(2,3,1),histogram(PH),title("Distribution Histogram of pH");
subplot(2,3,2),histogram(lPH),title("Distribution Histogram of Natural Logarithm of pH");
subplot(2,3,3),histogram(sqPH),title("Distributio n Histogram of Square Root of pH");
subplot(2,3,4),qqplot(PH),title("Quantile Quantile Plot of pH");
subplot(2,3,5),qqplot(lPH),title("Quantile Quantile Plot of Natural Logarithm of pH");
subplot(2,3,6),qqplot(sqPH),title("Quantile Quantile Plot of Square Root of pH");

2 间隔量算

  接下来，需求对筛选出的采样点相互之间的间隔加以量算。这是一个杂乱的进程，需求凭借循环句子。

  本部分详细代码如下。

poX=info(:,1);
poY=info(:,2);
dis=zeros(length(PH),length(PH));
for i=1:length(PH)
    for j=i+1:length(PH)
        dis(i,j)=sqrt((poX(i,1)-poX(j,1))^2+(poY(i,1)-poY(j,1))^2);
    end
end

3 间隔分组

  核算得到悉数采样点相互之间的间隔后，咱们需求根据必定的规模划定原则，对间隔数值加以分组。

  间隔分组首要需求确认步长。经过试验发现，若将步长选取过大会导致得到的散点图精度较低，而若步长选取过小则可能会使得每组点对总数量较少。因而，这儿取步长为500米；其次确认最大滞后距，这儿以悉数采样点间最大间隔的一半为其值。随后核算各组对应的滞后级别、各组上下界规模等。

  本部分详细代码附于本文4 均匀间隔、半方差核算及其绘图处。

4 均匀间隔、半方差核算及其绘图

  别离核算各个组内对应的点对个数、点对间间隔总和以及点对间特点值差值总和等。随后，根据上述参数，终究求出点对间间隔均匀值以及点对间特点值差值均匀值。

  根据各组对应点对间间隔均匀值为横轴，各组对应点对间特点值差值均匀值为纵轴，制作出经历半方差图。

  本部分及上述部分详细代码如下。

madi=max(max(dis));
midi=min(min(dis(dis>0)));
radi=madi-midi;
ste=500;
clnu=floor((madi/2)/ste)+1;
ponu=zeros(clnu,1);
todi=ponu;
todiav=todi;
diff=ponu;
diffav=diff;
for k=1:clnu
    midite=ste*(k-1);
    madite=ste*k;
    for i=1:length(sqPH)
        for j=i+1:length(sqPH)
            if dis(i,j)>midite && dis(i,j)<=madite
                ponu(k,1)=ponu(k,1)+1;
                todi(k,1)=todi(k,1)+dis(i,j); diff(k,1)=diff(k,1)+(sqPH(i)-sqPH(j))^2;
            end
        end
    end
    todiav(k,1)=todi(k,1)/ponu(k,1);
    diffav(k,1)=diff(k,1)/ponu(k,1)/2;
end
plot(todiav(:,1),diffav(:,1)),title("Empirical Semivariogram of Square Root of pH");
xlabel("Separation Distance (Metre)"),ylabel("Standardized Semivariance");

5 绘图成果

  经过上述进程，得到pH值开平方后的试验变异函数折线图及散点图。

  能够看到，pH值开平方后的试验变异函数较契合于有基台值的球状模型或指数模型。函数数值在间隔为0至8000米区间内快速上升，在间隔为8000米后数值上升放缓，变程为25000米左右；即其“先快速上升，再增速减缓，后趋于平稳”的图画全体趋势较为显着。但其数值全体表现较低——块金常数为0.004左右，而基台值仅为0.013左右。为验证数值正确性，相同对有机质、全氮进行上述全程操作。

  得到二者对应变异函数折线图与散点图。

  由以上三组、合计六幅的pH值开平方、有机质与全氮对应的试验变异函数折线图与散点图可知，不同数值对应试验变异函数数值的数量级亦会有所不同；但其全体“先快速上升，再增速减缓，后趋于平稳”的图画全体趋势是十分共同的。

  此外，如上文所说到的，针对三种空间特点数据（pH值、有机质含量与全氮含量）中最契合正态散布，亦是三种特点数据各三种（原始值、取对数与开平方）、共九种数据状况中唯一一个经过Lilliefors正态散布查验的数值——全氮含量经过异常值除掉后的原始值，将其正态散布的图画查验成果特展示如下。

至此，咱们就完成了悉数的操作、剖析进程~