首页天道酬勤遥感茶的半方差函数(半方差定义)

遥感茶的半方差函数(半方差定义)

admin 12-04 16:52 372次浏览

在之前的头条文章《插值、变异函数、克里格、线性无偏最优…地学计算概念及公式推导》中,我们详细介绍了地学计算的几个基本概念,并梳理了其数学推导公式。接下来,我将通过几篇新的标题文章详细练习和解释与地理计算相关的代码和操作。本文为第一章——:基于MATLAB的空间数据变异函数计算及经验半方差图绘制。

另一方面,由于上述文章涉及的相关理论概念相对抽象,往往需要结合实际才能更好地理解。因此,我们可以把上面提到的文章和这篇以及其他后来的地球科学计算文章一起看,以更好地理解相关理论的含义。

其中,因为本文使用的数据不是我的,所以很遗憾不能将数据展示给大家。但按照本文的思路和代码的详细讲解,完全可以用自己的数据重现空间数据变异函数计算和经验半变异函数绘制的全过程和分析方法。

1 数据处理

10-1010在本文中,我的初始数据是某地区658个土壤采样点的空间位置(x和y,以米为单位)、pH值、有机质含量和全氮含量。这些数据存储在“data.xls”文件中;而且后期的操作也比MATLAB软件中的多。因此,首先需要有选择地将源数据导入到MATLAB软件中。

这个功能可以通过在MATLAB软件中使用xlsread函数来实现。具体代码附在“1.3正态分布检验与变换”中。

1.1 数据读取

采样点数据可能因采样记录、实验室检测等过程存在一定误差,导致个别值出现异常。选择“平均加标准差法”来过滤和消除这些异常数据。

分别使用“平均加标准差法”中的“2S”和“3S”方法,发现“2S”方法优于后一种方法,因此后续实验继续使用“2S”方法。

其中,“2S”法是指数值大于或小于其平均值标准偏差2倍的部分视为异常值,“3S”法是指数值大于或小于其平均值标准偏差3倍的部分视为异常值。

得到异常值后,从658个采样点中剔除;剩余的采样点数据继续后续操作。

该特定代码的一部分附在“1.3正态分布检查和转换”中。

1.2 异常数据剔除

变差函数的计算应基于初始数据符合正态分布的假设;但采样点的数据不一定符合正态分布。因此,我们需要检验原始数据的正态分布。

正态分布检验一般可以通过数值检验、直方图、QQ图等图像直观判断。本文综合采用上述两种数值和图像检验方法,共同判断正态分布的特征。

对于数值试验方法,我一开始准备选择Kolmogorov-Smirnov试验方法;但知道这种方法只适用于标准正态试验,后来采用了Lilliefors试验。

Kolmogorov-Smirnov检验通过比较样本的经验分布函数和给定分布函数,推导出样本是否来自给定分布函数的总体。用于正态性检验时,只能做标准正态性检验。

Lilliefors检验改进了上述Kolmogorov-Smirnov检验,可用于一般正态分布检验。

QQ图(分位数分位数图)是散点图,横坐标代表一个样本数据的分位数,纵坐标代表另一个样本数据的分位数;由横坐标和纵坐标组成的散点图表示对应于相同累积概率的分位数。因此,QQ地图具有以下特点:

y=x

对于这条直线,如果散点图中的所有点都分布在直线附近,则两个样本的分布相等。因此,如果横坐标(纵坐标)表示为标准正态分布样本的分位数,那么散点图中的所有点都分布在上述直线附近,这表明纵坐标(纵坐标)表示的样本符合或近似符合正态分布。本文将横坐标表示为正态分布。

另外,PP图(概率概率图)也可以用来检验正态分布。PP的横坐标表示一个样本数据的累积概率,纵坐标表示另一个样本数据的累积概率。根据变量的累积概率,对应指定理论分布的累积概率,绘制散点图,用于直观检测样本数据是否满足一定的概率分布。与QQ图类似,如果测试数据符合规定的分布,其点都分布在上述直线附近。如果横坐标(纵坐标)表示为标准正态分布样本的分位数,那么散点图中的所有点都分布在直线附近,这表明纵坐标(纵坐标)表示的样本符合或基本接近正态分布。

三种土壤性质中,我选择以pH值为例先操作。通过上述数值测试和图像测试方法,发现没有异常值的原始pH值数据不符合正态分布。因此,尽量将原始数据转换成对数和平方根。随后发现,尽管原始pH平方数据的正态分布特征仍未能通过严格的Lilliefors检验,但其直方图和QQ图的图像检验结果更接近正态分布,也比前两者更为明显。因此,后续平方处理结果继续。

值得一提的是,本文后半部分在获得pH值平方根数据的实验变异函数和散点图后,对另外两类空间属性数据(即有机质含量和全氮含量)做同样的运算时,发现用“2S”法排除了全氮含量数据。

值后,其原始形式的数据是可以通过Lilliefors检验的,且其直方图、QQ图分布特点十分接近正态分布。

  我亦准备尝试对空间属性数据进行反正弦转换。但随后发现,已有三种属性数值的原始数据并不严格分布在-1至1的区间内,因此并未对其进行反正弦方式的转换。

  经过上述检验、转换处理过后的图像检验结果如下所示。

  以上部分代码如下:

1clc;clear; 2info=xlsread('data.xls'); 3oPH=info(:,3); 4oOM=info(:,4); 5oTN=info(:,5); 6 7mPH=mean(oPH); 8sPH=std(oPH); 9num2=find(oPH>(mPH+2*sPH)|oPH<(mPH-2*sPH)); 10num3=find(oPH>(mPH+3*sPH)|oPH<(mPH-3*sPH)); 11PH=oPH; 12for i=1:length(num2) 13 n=num2(i,1); 14 PH(n,:)=[0]; 15end 16PH(all(PH==0,2),:)=[]; 17 18%KSTest(PH,0.05) 19H1=lillietest(PH); 20 21for i=1:length(PH) 22 lPH(i,:)=log(PH(i,:)); 23end 24 25H2=lillietest(lPH); 26 27for i=1:length(PH) 28 sqPH(i,:)=(PH(i,:))^0.5; 29end 30 31H3=lillietest(sqPH); 32 33% for i=1:length(PH) 34% arcPH(i,:)=asin(PH(i,:)); 35% end 36% 37% H4=lillietest(arcPH); 38 39subplot(2,3,1),histogram(PH),title("Distribution Histogram of pH"); 40subplot(2,3,2),histogram(lPH),title("Distribution Histogram of Natural Logarithm of pH"); 41subplot(2,3,3),histogram(sqPH),title("Distributio n Histogram of Square Root of pH"); 42subplot(2,3,4),qqplot(PH),title("Quantile Quantile Plot of pH"); 43subplot(2,3,5),qqplot(lPH),title("Quantile Quantile Plot of Natural Logarithm of pH"); 44subplot(2,3,6),qqplot(sqPH),title("Quantile Quantile Plot of Square Root of pH");

2 距离量算

  接下来,需要对筛选出的采样点相互之间的距离加以量算。这是一个复杂的过程,需要借助循环语句。

  本部分具体代码如下。

1poX=info(:,1); 2poY=info(:,2); 3dis=zeros(length(PH),length(PH)); 4for i=1:length(PH) 5 for j=i+1:length(PH) 6 dis(i,j)=sqrt((poX(i,1)-poX(j,1))^2+(poY(i,1)-poY(j,1))^2); 7 end 8end

3 距离分组

  计算得到全部采样点相互之间的距离后,我们需要依据一定的范围划定原则,对距离数值加以分组。

  距离分组首先需要确定步长。经过实验发现,若将步长选取过大会导致得到的散点图精度较低,而若步长选取过小则可能会使得每组点对总数量较少。因此,这里取步长为500米;其次确定最大滞后距,这里以全部采样点间最大距离的一半为其值。随后计算各组对应的滞后级别、各组上下界范围等。

  本部分具体代码附于本文“4 平均距离、半方差计算及其绘图”处。

4 平均距离、半方差计算及其绘图

  分别计算各个组内对应的点对个数、点对间距离总和以及点对间属性值差值总和等。随后,依据上述参数,最终求出点对间距离平均值以及点对间属性值差值平均值。

  依据各组对应点对间距离平均值为横轴,各组对应点对间属性值差值平均值为纵轴,绘制出经验半方差图。

  本部分及上述部分具体代码如下。

1madi=max(max(dis)); 2midi=min(min(dis(dis>0))); 3radi=madi-midi; 4ste=500; 5clnu=floor((madi/2)/ste)+1; 6ponu=zeros(clnu,1); 7todi=ponu; 8todiav=todi; 9diff=ponu; 10diffav=diff; 11for k=1:clnu 12 midite=ste*(k-1); 13 madite=ste*k; 14 for i=1:length(sqPH) 15 for j=i+1:length(sqPH) 16 if dis(i,j)>midite && dis(i,j)<=madite 17 ponu(k,1)=ponu(k,1)+1; 18 todi(k,1)=todi(k,1)+dis(i,j); diff(k,1)=diff(k,1)+(sqPH(i)-sqPH(j))^2; 19 end 20 end 21 end 22 todiav(k,1)=todi(k,1)/ponu(k,1); 23 diffav(k,1)=diff(k,1)/ponu(k,1)/2; 24end 25plot(todiav(:,1),diffav(:,1)),title("Empirical Semivariogram of Square Root of pH"); 26xlabel("Separation Distance (Metre)"),ylabel("Standardized Semivariance");

5 绘图结果

  通过上述过程,得到pH值开平方后的实验变异函数折线图及散点图。

  可以看到,pH值开平方后的实验变异函数较符合于有基台值的球状模型或指数模型。函数数值在距离为0至8000米区间内快速上升,在距离为8000米后数值上升放缓,变程为25000米左右;即其“先快速上升,再增速减缓,后趋于平稳”的图像整体趋势较为明显。但其数值整体表现较低——块金常数为0.004左右,而基台值仅为0.013左右。为验证数值正确性,同样对有机质、全氮进行上述全程操作。  得到二者对应变异函数折线图与散点图。

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

  此外,如上文所提到的,针对三种空间属性数据(pH值、有机质含量与全氮含量)中最符合正态分布,亦是三种属性数据各三种(原始值、取对数与开平方)、共九种数据状态中唯一一个通过Lilliefors正态分布检验的数值——全氮含量经过异常值剔除后的原始值,将其正态分布的图像检验结果特展示如下。

至此,我们就完成了全部的操作、分析过程~

怎么使用jib插件为Java应用构建镜像让混合云变成“一朵云”React中的生命周期和子组件是什么怎么使用uniapp自定义弹框然后返回undefined
啥啥真不少(低学历适合的技术工种) ()
相关内容