Matlab計(jì)算變異函數(shù)并繪制經(jīng)驗(yàn)半方差圖詳解
本文介紹基于MATLAB求取空間數(shù)據(jù)的變異函數(shù),并繪制經(jīng)驗(yàn)半方差圖的方法。
由于本文所用的數(shù)據(jù)并不是我的,因此遺憾不能將數(shù)據(jù)一并展示給大家;但是依據(jù)本篇博客的思想與對(duì)代碼的詳細(xì)解釋,大家用自己的數(shù)據(jù),可以將空間數(shù)據(jù)變異函數(shù)計(jì)算與經(jīng)驗(yàn)半方差圖繪制的全部過程與分析方法加以完整重現(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è)試等過程,可能具有一定誤差,從而出現(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è)采樣點(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)可以通過數(shù)值檢驗(yàn)與直方圖、QQ圖等圖像加以直觀判斷。本文綜合采取以上兩種數(shù)值、圖像檢驗(yàn)方法,共同判斷正態(tài)分布特性。
針對(duì)數(shù)值檢驗(yàn)方法,我在一開始準(zhǔn)備選擇采用Kolmogorov-Smirnov檢驗(yàn)方法;但由于了解到,這一方法僅僅適用于標(biāo)準(zhǔn)正態(tài)檢驗(yàn),因此隨后改用Lilliefors檢驗(yàn)。
Kolmogorov-Smirnov檢驗(yàn)通過樣本的經(jīng)驗(yàn)分布函數(shù)與給定分布函數(shù)的比較,推斷該樣本是否來自給定分布函數(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):針對(duì)y=x
這一直線,若散點(diǎn)圖中各點(diǎn)均在直線附近分布,則說明兩個(gè)樣本為同等分布;因此,若將橫坐標(biāo)(縱坐標(biāo))表示為一個(gè)標(biāo)準(zhǔn)正態(tài)分布樣本的分位數(shù),則散點(diǎn)圖中各點(diǎn)均在上述直線附近分布可以說明,縱坐標(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圖類似,如果被檢驗(yàn)的數(shù)據(jù)符合所指定的分布,則其各點(diǎn)均在上述直線附近分布。若將橫坐標(biāo)(縱坐標(biāo))表示為一個(gè)標(biāo)準(zhǔn)正態(tài)分布樣本的分位數(shù),則散點(diǎn)圖中各點(diǎn)均在直線附近分布可以說明,縱坐標(biāo)(橫坐標(biāo))表示的樣本符合或基本近似符合正態(tài)分布。
三種土壤屬性,我選擇首先以pH數(shù)值為例進(jìn)行操作。通過上述數(shù)值檢驗(yàn)、圖像檢驗(yàn)方法,檢驗(yàn)得到剔除異常值后的原始pH數(shù)值數(shù)據(jù)并不符合正態(tài)分布這一結(jié)論。因此,嘗試對(duì)原數(shù)據(jù)加以對(duì)數(shù)、開平方等轉(zhuǎn)換處理;隨后發(fā)現(xiàn),原始pH值開平方數(shù)據(jù)的正態(tài)分布特征雖然依舊無法通過較為嚴(yán)格的Lilliefors檢驗(yàn),但其直方圖、QQ圖的圖像檢驗(yàn)結(jié)果較為接近正態(tài)分布,并較之前二者更加明顯。故后續(xù)取開平方處理結(jié)果繼續(xù)進(jìn)行。
值得一提的是,本文后半部分得到pH值開平方數(shù)據(jù)的實(shí)驗(yàn)變異函數(shù)及其散點(diǎn)圖后,在對(duì)其余兩種空間屬性數(shù)據(jù)(即有機(jī)質(zhì)含量與全氮含量)進(jìn)行同樣的操作時(shí),發(fā)現(xiàn)全氮含量數(shù)據(jù)在經(jīng)過“2S”方法剔除異常值后,其原始形式的數(shù)據(jù)是可以通過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)過上述檢驗(yàn)、轉(zhuǎn)換處理過后的圖像檢驗(yàn)結(jié)果如下所示。
以上部分代碼如下:
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 距離量算
接下來,需要對(duì)篩選出的采樣點(diǎn)相互之間的距離加以量算。這是一個(gè)復(fù)雜的過程,需要借助循環(huán)語句。
本部分具體代碼如下。
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)過實(shí)驗(yàn)發(fā)現(xiàn),若將步長(zhǎng)選取過大會(huì)導(dǎo)致得到的散點(diǎn)圖精度較低,而若步長(zhǎng)選取過小則可能會(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é)果
通過上述過程,得到pH值開平方后的實(shí)驗(yàn)變異函數(shù)折線圖及散點(diǎn)圖。
可以看到,pH值開平方后的實(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)圖。
由以上三組、共計(jì)六幅的pH值開平方、有機(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ù)與開平方)、共九種數(shù)據(jù)狀態(tài)中唯一一個(gè)通過Lilliefors正態(tài)分布檢驗(yàn)的數(shù)值——全氮含量經(jīng)過異常值剔除后的原始值,將其正態(tài)分布的圖像檢驗(yàn)結(jié)果特展示如下。
至此,我們就完成了全部的操作、分析過程~
到此這篇關(guān)于Matlab計(jì)算變異函數(shù)并繪制經(jīng)驗(yàn)半方差圖詳解的文章就介紹到這了,更多相關(guān)Matlab計(jì)算變異函數(shù) 繪制經(jīng)驗(yàn)半方差圖內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家!
相關(guān)文章
C語言中你容易忽略的知識(shí)點(diǎn)與技巧總結(jié)
這篇文章主要給大家介紹了關(guān)于C語言中你容易忽略的知識(shí)點(diǎn)與技巧,文中通過實(shí)例代碼以及圖文介紹的非常詳細(xì),對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友可以參考下2022-03-03詳解C語言中accept()函數(shù)和shutdown()函數(shù)的使用
這篇文章主要介紹了詳解C語言中accept()函數(shù)和shutdown()函數(shù)的使用,用來操作socket相關(guān)的網(wǎng)絡(luò)通信,需要的朋友可以參考下2015-09-09C語言重難點(diǎn)之內(nèi)存對(duì)齊和位段
這篇文章主要介紹了C語言重難點(diǎn)之內(nèi)存對(duì)齊和位段,文中通過示例代碼介紹的非常詳細(xì),對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2021-05-05C語言實(shí)現(xiàn)“幸運(yùn)數(shù)”的實(shí)例詳解
這篇文章主要介紹了C語言實(shí)現(xiàn)“幸運(yùn)數(shù)”的實(shí)例詳解的相關(guān)資料,需要的朋友可以參考下2017-07-07C語言實(shí)現(xiàn)經(jīng)典排序算法的示例代碼
這篇文章主要為大家詳細(xì)介紹了如何利用C語言實(shí)現(xiàn)經(jīng)典排序算法中的冒泡排序、選擇排序、插入排序、希爾排序,文中的示例代碼講解詳細(xì),需要的可以參考一下2022-08-08C++類型轉(zhuǎn)換運(yùn)算符的實(shí)例詳解
這篇文章主要介紹了C++類型轉(zhuǎn)換運(yùn)算符的實(shí)例詳解的相關(guān)資料,希望通過本文大家能夠掌握這部分內(nèi)容,需要的朋友可以參考下2017-09-09