??在前期的博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/113943720）中希太，我們?cè)敿?xì)介紹了地學(xué)計(jì)算的幾個(gè)基本概念录淡，并對(duì)其數(shù)學(xué)推導(dǎo)公式加以了梳理。接下來(lái)阔涉，我將通過(guò)幾篇新的專(zhuān)題博客，對(duì)地學(xué)計(jì)算相關(guān)的代碼耳幢、操作加以實(shí)踐與詳細(xì)講解迂苛。本篇博客便是第一篇——基于MATLAB的空間數(shù)據(jù)變異函數(shù)計(jì)算與經(jīng)驗(yàn)半方差圖繪制。

??另一方面弃揽，由于上述博客所涉及的相關(guān)理論概念較為抽象脯爪，往往需要結(jié)合實(shí)踐才可以更好理解，因此大家可以將上述博客與本篇及后期的其它地學(xué)計(jì)算博客一同來(lái)看矿微，可以更好理解相關(guān)理論的含義痕慢。

??其中，由于本文所用的數(shù)據(jù)并不是我的涌矢，因此遺憾不能將數(shù)據(jù)一并展示給大家掖举；但是依據(jù)本篇博客的思想與對(duì)代碼的詳細(xì)解釋?zhuān)蠹矣米约旱臄?shù)據(jù)，可以將空間數(shù)據(jù)變異函數(shù)計(jì)算與經(jīng)驗(yàn)半方差圖繪制的全部過(guò)程與分析方法加以完整重現(xiàn)娜庇。

1 數(shù)據(jù)處理

1.1 數(shù)據(jù)讀取

??本文中塔次，我的初始數(shù)據(jù)為某區(qū)域658個(gè)土壤采樣點(diǎn)的空間位置（X與Y，單位為米）名秀、pH值励负、有機(jī)質(zhì)含量與全氮含量。這些數(shù)據(jù)均存儲(chǔ)于“data.xls”文件中泰偿；而后期操作多于MATLAB軟件中進(jìn)行熄守。因此，首先需將源數(shù)據(jù)選擇性地導(dǎo)入MATLAB軟件中。

??利用MATLAB軟件中xlsread函數(shù)可以實(shí)現(xiàn)這一功能裕照。具體代碼附于“1.3 正態(tài)分布檢驗(yàn)及轉(zhuǎn)換”處攒发。

1.2 異常數(shù)據(jù)剔除

??得到的采樣點(diǎn)數(shù)據(jù)由于采樣記錄、實(shí)驗(yàn)室測(cè)試等過(guò)程晋南，可能具有一定誤差惠猿，從而出現(xiàn)個(gè)別異常值。選用“平均值加標(biāo)準(zhǔn)差法”對(duì)這些異常數(shù)據(jù)加以篩選负间、剔除偶妖。

??分別利用“平均值加標(biāo)準(zhǔn)差法”中“2S”與“3S”方法加以處理，發(fā)現(xiàn)“2S”方法處理效果相對(duì)后者較好政溃，故后續(xù)實(shí)驗(yàn)取“2S”方法處理結(jié)果繼續(xù)進(jìn)行趾访。

??其中，“2S”方法是指將數(shù)值大于或小于其平均值±2倍標(biāo)準(zhǔn)差的部分視作異常值董虱，“3S”方法則是指將數(shù)值大于或小于其平均值±3倍標(biāo)準(zhǔn)差的部分視作異常值扼鞋。

??得到異常值后，將其從658個(gè)個(gè)采樣點(diǎn)中剔除；剩余的采樣點(diǎn)數(shù)據(jù)繼續(xù)后續(xù)操作。

??本部分具體代碼附于“1.3 正態(tài)分布檢驗(yàn)及轉(zhuǎn)換”處绍赛。

1.3 正態(tài)分布檢驗(yàn)及轉(zhuǎn)換

??計(jì)算變異函數(shù)需建立在初始數(shù)據(jù)符合正態(tài)分布的假設(shè)之上；而采樣點(diǎn)數(shù)據(jù)并不一定符合正態(tài)分布溃槐。因此，我們需要對(duì)原始數(shù)據(jù)加以正態(tài)分布檢驗(yàn)科吭。

??一般地昏滴，正態(tài)分布檢驗(yàn)可以通過(guò)數(shù)值檢驗(yàn)與直方圖、QQ圖等圖像加以直觀判斷砌溺。本文綜合采取以上兩種數(shù)值影涉、圖像檢驗(yàn)方法变隔，共同判斷正態(tài)分布特性规伐。

??針對(duì)數(shù)值檢驗(yàn)方法，我在一開(kāi)始準(zhǔn)備選擇采用Kolmogorov-Smirnov檢驗(yàn)方法匣缘；但由于了解到猖闪，這一方法僅僅適用于標(biāo)準(zhǔn)正態(tài)檢驗(yàn)，因此隨后改用Lilliefors檢驗(yàn)肌厨。

??Kolmogorov-Smirnov檢驗(yàn)通過(guò)樣本的經(jīng)驗(yàn)分布函數(shù)與給定分布函數(shù)的比較培慌，推斷該樣本是否來(lái)自給定分布函數(shù)的總體；當(dāng)其用于正態(tài)性檢驗(yàn)時(shí)只能做標(biāo)準(zhǔn)正態(tài)檢驗(yàn)柑爸。

??Lilliefors檢驗(yàn)則將上述Kolmogorov-Smirnov檢驗(yàn)改進(jìn)吵护，其可用于一般的正態(tài)分布檢驗(yàn)。

??QQ圖（Quantile Quantile Plot）是一種散點(diǎn)圖，其橫坐標(biāo)表示某一樣本數(shù)據(jù)的分位數(shù)馅而，縱坐標(biāo)則表示另一樣本數(shù)據(jù)的分位數(shù)祥诽；橫坐標(biāo)與縱坐標(biāo)組成的散點(diǎn)圖代表同一個(gè)累計(jì)概率所對(duì)應(yīng)的分位數(shù)。因此瓮恭，QQ圖具有這樣的特點(diǎn)：

??y=x

??針對(duì)這一直線雄坪，若散點(diǎn)圖中各點(diǎn)均在直線附近分布，則說(shuō)明兩個(gè)樣本為同等分布屯蹦；因此维哈，若將橫坐標(biāo)（縱坐標(biāo)）表示為一個(gè)標(biāo)準(zhǔn)正態(tài)分布樣本的分位數(shù)，則散點(diǎn)圖中各點(diǎn)均在上述直線附近分布可以說(shuō)明登澜，縱坐標(biāo)（橫坐標(biāo)）表示的樣本符合或基本近似符合正態(tài)分布阔挠。本文采用將橫坐標(biāo)表示為正態(tài)分布的方式。

??此外脑蠕，PP圖（Probability Probability Plot）同樣可以用于正態(tài)分布的檢驗(yàn)谒亦。PP圖橫坐標(biāo)表示某一樣本數(shù)據(jù)的累積概率，縱坐標(biāo)則表示另一樣本數(shù)據(jù)的累積概率空郊；其根據(jù)變量的累積概率對(duì)應(yīng)于所指定的理論分布累積概率并繪制的散點(diǎn)圖份招，用于直觀地檢測(cè)樣本數(shù)據(jù)是否符合某一概率分布。和QQ圖類(lèi)似狞甚，如果被檢驗(yàn)的數(shù)據(jù)符合所指定的分布锁摔，則其各點(diǎn)均在上述直線附近分布。若將橫坐標(biāo)（縱坐標(biāo)）表示為一個(gè)標(biāo)準(zhǔn)正態(tài)分布樣本的分位數(shù)哼审，則散點(diǎn)圖中各點(diǎn)均在直線附近分布可以說(shuō)明谐腰，縱坐標(biāo)（橫坐標(biāo)）表示的樣本符合或基本近似符合正態(tài)分布。

??三種土壤屬性涩盾，我選擇首先以pH數(shù)值為例進(jìn)行操作十气。通過(guò)上述數(shù)值檢驗(yàn)、圖像檢驗(yàn)方法春霍，檢驗(yàn)得到剔除異常值后的原始pH數(shù)值數(shù)據(jù)并不符合正態(tài)分布這一結(jié)論砸西。因此，嘗試對(duì)原數(shù)據(jù)加以對(duì)數(shù)址儒、開(kāi)平方等轉(zhuǎn)換處理芹枷；隨后發(fā)現(xiàn)，原始pH值開(kāi)平方數(shù)據(jù)的正態(tài)分布特征雖然依舊無(wú)法通過(guò)較為嚴(yán)格的Lilliefors檢驗(yàn)莲趣，但其直方圖鸳慈、QQ圖的圖像檢驗(yàn)結(jié)果較為接近正態(tài)分布，并較之前二者更加明顯喧伞。故后續(xù)取開(kāi)平方處理結(jié)果繼續(xù)進(jìn)行走芋。

??值得一提的是绩郎，本文后半部分得到pH值開(kāi)平方數(shù)據(jù)的實(shí)驗(yàn)變異函數(shù)及其散點(diǎn)圖后，在對(duì)其余兩種空間屬性數(shù)據(jù)（即有機(jī)質(zhì)含量與全氮含量）進(jìn)行同樣的操作時(shí)翁逞，發(fā)現(xiàn)全氮含量數(shù)據(jù)在經(jīng)過(guò)“2S”方法剔除異常值后嗽上，其原始形式的數(shù)據(jù)是可以通過(guò)Lilliefors檢驗(yàn)的，且其直方圖熄攘、QQ圖分布特點(diǎn)十分接近正態(tài)分布兽愤。

??我亦準(zhǔn)備嘗試對(duì)空間屬性數(shù)據(jù)進(jìn)行反正弦轉(zhuǎn)換。但隨后發(fā)現(xiàn)挪圾，已有三種屬性數(shù)值的原始數(shù)據(jù)并不嚴(yán)格分布在-1至1的區(qū)間內(nèi)浅萧，因此并未對(duì)其進(jìn)行反正弦方式的轉(zhuǎn)換。

??經(jīng)過(guò)上述檢驗(yàn)哲思、轉(zhuǎn)換處理過(guò)后的圖像檢驗(yàn)結(jié)果如下所示洼畅。

image

??以上部分代碼如下：

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 距離量算

??接下來(lái)，需要對(duì)篩選出的采樣點(diǎn)相互之間的距離加以量算棚赔。這是一個(gè)復(fù)雜的過(guò)程帝簇，需要借助循環(huán)語(yǔ)句。

??本部分具體代碼如下靠益。

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 距離分組

??計(jì)算得到全部采樣點(diǎn)相互之間的距離后丧肴，我們需要依據(jù)一定的范圍劃定原則，對(duì)距離數(shù)值加以分組胧后。

??距離分組首先需要確定步長(zhǎng)芋浮。經(jīng)過(guò)實(shí)驗(yàn)發(fā)現(xiàn)，若將步長(zhǎng)選取過(guò)大會(huì)導(dǎo)致得到的散點(diǎn)圖精度較低壳快，而若步長(zhǎng)選取過(guò)小則可能會(huì)使得每組點(diǎn)對(duì)總數(shù)量較少纸巷。因此，這里取步長(zhǎng)為500米眶痰；其次確定最大滯后距瘤旨，這里以全部采樣點(diǎn)間最大距離的一半為其值。隨后計(jì)算各組對(duì)應(yīng)的滯后級(jí)別竖伯、各組上下界范圍等存哲。

??本部分具體代碼附于本文“4 平均距離、半方差計(jì)算及其繪圖”處黔夭。

4 平均距離宏胯、半方差計(jì)算及其繪圖

??分別計(jì)算各個(gè)組內(nèi)對(duì)應(yīng)的點(diǎn)對(duì)個(gè)數(shù)、點(diǎn)對(duì)間距離總和以及點(diǎn)對(duì)間屬性值差值總和等本姥。隨后，依據(jù)上述參數(shù)杭棵，最終求出點(diǎn)對(duì)間距離平均值以及點(diǎn)對(duì)間屬性值差值平均值婚惫。

??依據(jù)各組對(duì)應(yīng)點(diǎn)對(duì)間距離平均值為橫軸氛赐，各組對(duì)應(yīng)點(diǎn)對(duì)間屬性值差值平均值為縱軸，繪制出經(jīng)驗(yàn)半方差圖先舷。

??本部分及上述部分具體代碼如下艰管。

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 繪圖結(jié)果

??通過(guò)上述過(guò)程，得到pH值開(kāi)平方后的實(shí)驗(yàn)變異函數(shù)折線圖及散點(diǎn)圖蒋川。

image

image

??可以看到牲芋，pH值開(kāi)平方后的實(shí)驗(yàn)變異函數(shù)較符合于有基臺(tái)值的球狀模型或指數(shù)模型。函數(shù)數(shù)值在距離為0至8000米區(qū)間內(nèi)快速上升捺球，在距離為8000米后數(shù)值上升放緩缸浦，變程為25000米左右；即其“先快速上升氮兵，再增速減緩裂逐，后趨于平穩(wěn)”的圖像整體趨勢(shì)較為明顯。但其數(shù)值整體表現(xiàn)較低——塊金常數(shù)為0.004左右泣栈，而基臺(tái)值僅為0.013左右卜高。為驗(yàn)證數(shù)值正確性，同樣對(duì)有機(jī)質(zhì)南片、全氮進(jìn)行上述全程操作掺涛。
??得到二者對(duì)應(yīng)變異函數(shù)折線圖與散點(diǎn)圖。

image

image

image

image

??由以上三組疼进、共計(jì)六幅的pH值開(kāi)平方鸽照、有機(jī)質(zhì)與全氮對(duì)應(yīng)的實(shí)驗(yàn)變異函數(shù)折線圖與散點(diǎn)圖可知，不同數(shù)值對(duì)應(yīng)實(shí)驗(yàn)變異函數(shù)數(shù)值的數(shù)量級(jí)亦會(huì)有所不同颠悬；但其整體“先快速上升矮燎，再增速減緩，后趨于平穩(wěn)”的圖像整體趨勢(shì)是十分一致的赔癌。

??此外诞外，如上文所提到的，針對(duì)三種空間屬性數(shù)據(jù)（pH值灾票、有機(jī)質(zhì)含量與全氮含量）中最符合正態(tài)分布峡谊，亦是三種屬性數(shù)據(jù)各三種（原始值、取對(duì)數(shù)與開(kāi)平方）刊苍、共九種數(shù)據(jù)狀態(tài)中唯一一個(gè)通過(guò)Lilliefors正態(tài)分布檢驗(yàn)的數(shù)值——全氮含量經(jīng)過(guò)異常值剔除后的原始值既们，將其正態(tài)分布的圖像檢驗(yàn)結(jié)果特展示如下。

image

至此正什，我們就完成了全部的操作啥纸、分析過(guò)程~