技術文章 » MATLAB計算變異函數並繪製經驗半方差圖

MATLAB計算變異函數並繪製經驗半方差圖

2023-04-03 15:00:38

  本文介紹基於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正態分佈檢驗的數值——全氮含量經過異常值剔除後的原始值，將其正態分佈的影象檢驗結果特展示如下。

至此，我們就完成了全部的操作、分析過程~