利用R語(yǔ)言進(jìn)行地統(tǒng)計(jì)學(xué)的變異函數(shù)教學(xué)

楊孟 刁一偉

1.無(wú)錫學(xué)院環(huán)境工程學(xué)院 江蘇無(wú)錫 214105；2.無(wú)錫學(xué)院大氣與遙感學(xué)院 江蘇無(wú)錫 214105

1 概述

地質(zhì)統(tǒng)計(jì)學(xué)，簡(jiǎn)稱地統(tǒng)計(jì)學(xué)，發(fā)源于礦業(yè)領(lǐng)域，但理論上可以對(duì)任何空間上即具有隨機(jī)性又具有結(jié)構(gòu)性的自然現(xiàn)象進(jìn)行研究[1]。近年來(lái)，地統(tǒng)計(jì)學(xué)的應(yīng)用更加廣泛，從傳統(tǒng)的地學(xué)和環(huán)境生態(tài)學(xué)領(lǐng)域進(jìn)一步擴(kuò)展到醫(yī)學(xué)領(lǐng)域。

地統(tǒng)計(jì)學(xué)研究的是空間變量，因此，無(wú)論是理論還是實(shí)踐教學(xué)都需要基于具有空間分析功能的軟件來(lái)開(kāi)展。在地統(tǒng)計(jì)學(xué)的教學(xué)中，通常采用的是ArcGIS軟件[2]。ArcGIS軟件作為一種圖形界面軟件，具有直觀和操作簡(jiǎn)易的特點(diǎn)。但是，ArcGIS是一種收費(fèi)軟件，需要找出一種免費(fèi)的替代產(chǎn)品。R語(yǔ)言[3]是統(tǒng)計(jì)學(xué)領(lǐng)域廣泛使用的軟件，它是基于通用公共許可協(xié)議(GNU，General Public License)進(jìn)行數(shù)據(jù)處理、可視化和統(tǒng)計(jì)學(xué)分析的開(kāi)源解釋型語(yǔ)言，正好是一種很好的替代軟件。

本文闡釋了“地統(tǒng)計(jì)學(xué)”課程知識(shí)內(nèi)容架構(gòu)，并結(jié)合案例介紹了R語(yǔ)言在變異函數(shù)的理論和實(shí)踐學(xué)中的應(yīng)用，旨在展現(xiàn)R語(yǔ)言在地統(tǒng)計(jì)學(xué)分析中的功能，激發(fā)學(xué)生借助開(kāi)源軟件進(jìn)行地統(tǒng)計(jì)分析的興趣。

2 “地統(tǒng)計(jì)學(xué)”課程內(nèi)容架構(gòu)

“地統(tǒng)計(jì)學(xué)”是我校地理信息科學(xué)專業(yè)的專業(yè)課，本課程包括五個(gè)部分的內(nèi)容。第一部分是回顧單變量和雙變量統(tǒng)計(jì)分析的基礎(chǔ)知識(shí)，介紹空間抽樣的類型及其特征。第二部分是空間分析，包括傳統(tǒng)空間分析以及空間自相關(guān)分析。第三部分是空間結(jié)構(gòu)建模，包括空間變異的通用模型，隨機(jī)場(chǎng)理論以及變異函數(shù)模型及其擬合。第四部分是克里金法，介紹常用克里金法的特點(diǎn)和基本假設(shè)。第五部分是通過(guò)交叉驗(yàn)證或者獨(dú)立驗(yàn)證數(shù)據(jù)集對(duì)模型進(jìn)行評(píng)價(jià)與優(yōu)化。

3 R語(yǔ)言在地統(tǒng)計(jì)分析中的優(yōu)勢(shì)特征

R語(yǔ)言為開(kāi)源免費(fèi)軟件，其基本內(nèi)核R及作為其集成開(kāi)發(fā)環(huán)境的R Studio，均可以免費(fèi)下載和使用。與商業(yè)化的GIS軟件相比，它的安裝程序小巧，對(duì)電腦硬件要求低，具有多種操作系統(tǒng)版本，在系統(tǒng)中裝載便捷且穩(wěn)定可靠。

R語(yǔ)言具有多個(gè)地統(tǒng)計(jì)擴(kuò)展包。R語(yǔ)言作為開(kāi)源的解釋型語(yǔ)言，擁有數(shù)萬(wàn)個(gè)用戶制作和貢獻(xiàn)的擴(kuò)展程序包，其內(nèi)容涵蓋諸多領(lǐng)域，可以實(shí)現(xiàn)特定目的。這些擴(kuò)展程序包可以在官方網(wǎng)站或者國(guó)內(nèi)鏡像網(wǎng)站免費(fèi)下載。對(duì)于地統(tǒng)計(jì)分析，R語(yǔ)言最常用的擴(kuò)展程序包是gstat，能夠系統(tǒng)實(shí)現(xiàn)地統(tǒng)計(jì)分析的常用功能，還能與空間數(shù)據(jù)處理擴(kuò)展包sp和sf兼容。另一個(gè)實(shí)現(xiàn)地統(tǒng)計(jì)分析的擴(kuò)展包是georob，它具有以下特點(diǎn)：具有針對(duì)異常值的穩(wěn)健地統(tǒng)計(jì)方法；使用最大似然法對(duì)參數(shù)進(jìn)行估計(jì)。R語(yǔ)言中實(shí)現(xiàn)地統(tǒng)計(jì)功能的擴(kuò)展包還有g(shù)eoR和geoRglm。

地統(tǒng)計(jì)分析的采樣點(diǎn)數(shù)據(jù)通常為矢量數(shù)據(jù)，因此軟件對(duì)矢量數(shù)據(jù)的處理和分析能力尤為重要。R語(yǔ)言具有操作矢量數(shù)據(jù)的擴(kuò)展包sf，其引入了空間數(shù)量分析領(lǐng)域通用的sf(simple features)標(biāo)準(zhǔn)規(guī)范，降低了處理、轉(zhuǎn)化與繪制地理空間數(shù)據(jù)的復(fù)雜度。

4 基于R語(yǔ)言的變異函數(shù)案例教學(xué)

該教學(xué)內(nèi)容運(yùn)用R語(yǔ)言實(shí)現(xiàn)經(jīng)驗(yàn)變異函數(shù)的計(jì)算、可視化以及擬合。使用的R擴(kuò)展包有g(shù)stat、sf、sp以及ggplot2。首先，需要安裝所需的擴(kuò)展包及其依賴，并將擴(kuò)展包加載到工作區(qū)，代碼為：

4.1 加載和查看數(shù)據(jù)

在本案例中，將使用隨sp包分發(fā)的Meuse土壤污染數(shù)據(jù)集。該數(shù)據(jù)集由表層土樣本的重金屬含量以及樣本空間位置數(shù)據(jù)。利用data函數(shù)將Meuse數(shù)據(jù)集加載到工作區(qū)中：

4.2 將數(shù)據(jù)設(shè)置為空間顯式的R對(duì)象

Meuse數(shù)據(jù)集為類似屬性表的數(shù)據(jù)框格式，坐標(biāo)還只是數(shù)據(jù)框中的字段。需要用到sf包的st_as_sf函數(shù)指定代表坐標(biāo)的字段，將數(shù)據(jù)框轉(zhuǎn)換為空間顯式的sf對(duì)象；保留原始的數(shù)據(jù)框，并制作一個(gè)名為meuse.sf的sf對(duì)象：

從元數(shù)據(jù)中得知Meuse數(shù)據(jù)集的坐標(biāo)系的EPSG代碼為28992。使用st_crs函數(shù)為sf對(duì)象定義地理參考系統(tǒng)：

接下來(lái)，繪制鋅(zinc)含量的符號(hào)圖，即繪制樣點(diǎn)位置并用符號(hào)大小表示鋅含量大小，結(jié)果如圖1所示：

圖1 Meuse數(shù)據(jù)集中金屬鋅含量的符號(hào)圖

4.3 變異函數(shù)建模

4.3.1 計(jì)算點(diǎn)對(duì)的半方差

空間中的任意兩點(diǎn)構(gòu)成一個(gè)點(diǎn)對(duì)，計(jì)算Meuse數(shù)據(jù)集的點(diǎn)最多能形成多少點(diǎn)對(duì)，代碼如下所示，結(jié)果表明最多形成11935個(gè)點(diǎn)對(duì)。

半方差是指對(duì)構(gòu)成點(diǎn)對(duì)的兩個(gè)點(diǎn)之間屬性值差異的數(shù)學(xué)度量。每個(gè)點(diǎn)對(duì)都可以計(jì)算一個(gè)半方差(γ)：

(1)

上式中，X表示觀測(cè)點(diǎn)的位置，z(X)表示觀測(cè)點(diǎn)的屬性值。利用下面的代碼計(jì)算數(shù)據(jù)集中前兩個(gè)點(diǎn)的距離和半方差，可知該點(diǎn)對(duì)的距離和半方差分別約為70.84m以及0.001(mg/kg)2。

4.3.2 繪制變異函數(shù)云圖

可以通過(guò)變異函數(shù)云展示所有點(diǎn)對(duì)的半方差。以下代碼以logZn為變量，利用variogram函數(shù)計(jì)算距離小于200m的點(diǎn)對(duì)的經(jīng)驗(yàn)變異函數(shù)云。可選參數(shù)cloud設(shè)置為TRUE時(shí)計(jì)算變差函數(shù)云；cutoff設(shè)置點(diǎn)對(duì)的最大距離。利用函數(shù)print繪制變異函數(shù)云圖，結(jié)果如圖2(a)所示。

(a)最大距離200m

當(dāng)最大距離增加時(shí)，見(jiàn)圖2(b)，在半方差低值區(qū)域的許多點(diǎn)幾乎相互重疊，此時(shí)很難通過(guò)變異函數(shù)云看出點(diǎn)對(duì)距離和半方差之間存在的關(guān)系。為了更好地可視化這種關(guān)系，常用方法是將點(diǎn)對(duì)進(jìn)行分組，計(jì)算每組的半方差的算術(shù)平均值，這就是經(jīng)典的馬瑟隆經(jīng)驗(yàn)變異函數(shù)：

(2)

上式中，從編號(hào)為1，2，…，i，…，j，…n的點(diǎn)構(gòu)建點(diǎn)對(duì)；符號(hào)(i，j)|hij=h表示編號(hào)為i和j的點(diǎn)對(duì)之間的分隔向量為h。

4.3.3 確定點(diǎn)對(duì)分組的步長(zhǎng)

式(2)中分隔向量h實(shí)際上是某個(gè)距離范圍，而非一個(gè)特定值。通常將距離在一定范圍內(nèi)的點(diǎn)對(duì)作為一組來(lái)考慮，用來(lái)分組的距離范圍稱為步長(zhǎng)。例如，將“距離為0～50m的點(diǎn)對(duì)”分為第1組，“距離為50～100m的點(diǎn)對(duì)”分為第2組，依次類推一直到最大距離1500m。此時(shí)，步長(zhǎng)為50m、步數(shù)(組數(shù))為30。

經(jīng)驗(yàn)變異函數(shù)受到點(diǎn)對(duì)分組的顯著影響。同樣的變異函數(shù)云，選擇的步長(zhǎng)不同，得到的經(jīng)驗(yàn)變異函數(shù)圖也不同。步長(zhǎng)越寬、步數(shù)越少的分組方式，具有越少的細(xì)節(jié)，越多的噪音被平滑。相反，步長(zhǎng)越窄、步數(shù)越多的分組，具有越多的細(xì)節(jié)和噪音。R語(yǔ)言的variogram函數(shù)設(shè)定步長(zhǎng)的經(jīng)驗(yàn)法則是用最大點(diǎn)對(duì)距離的1/3除以15得到步長(zhǎng)，根據(jù)該原則，Meuse數(shù)據(jù)集的步長(zhǎng)大小應(yīng)設(shè)置為100m。

以下代碼針對(duì)變量logZn，繪制了一個(gè)步長(zhǎng)為100m，最大距離1500m的經(jīng)驗(yàn)變異函數(shù)，結(jié)果如圖3(a)所示。其中，通過(guò)參數(shù)width設(shè)置步長(zhǎng)，cutoff設(shè)置最大距離(等于步長(zhǎng)乘以步數(shù))；可選參數(shù)plot.numbers設(shè)置為TRUE則在點(diǎn)旁標(biāo)準(zhǔn)每個(gè)組的點(diǎn)對(duì)數(shù)。

圖3中還給出了相同數(shù)據(jù)、設(shè)置不同步長(zhǎng)時(shí)的經(jīng)驗(yàn)變異函數(shù)圖。由圖可知，在相同的最大距離下，不同的步長(zhǎng)大小對(duì)經(jīng)驗(yàn)變異函數(shù)具有顯著影響。步長(zhǎng)100m最為合適，可以在沒(méi)有噪音的情況下清楚地看出基臺(tái)值、變程和塊金值。太大的步長(zhǎng)(300m)導(dǎo)致細(xì)節(jié)信息被過(guò)度平滑。步長(zhǎng)太小(50m)則噪音太多，變異函數(shù)的基本特征幾乎難以識(shí)別，且有些組內(nèi)的點(diǎn)對(duì)數(shù)太少。據(jù)此，選擇100m作為最佳的步長(zhǎng)。

(a)100m

4.3.4 通過(guò)經(jīng)驗(yàn)變異函數(shù)分析空間自相關(guān)性

經(jīng)驗(yàn)變異函數(shù)圖圖3(a)提供了空間自相關(guān)的以下證據(jù)：點(diǎn)對(duì)間距越小，半方差越小；在某個(gè)距離范圍內(nèi)，半方差隨著距離增加而增大，直到一個(gè)平穩(wěn)狀態(tài)，此時(shí)的距離稱為變程；超過(guò)變程不再有空間自相關(guān)；塊金值與基臺(tái)值的比值(塊基比)表示樣本數(shù)據(jù)無(wú)法解釋的空間變異的占比；塊基比約為0.17，表明結(jié)構(gòu)基臺(tái)值的占比較大，即空間自相關(guān)的強(qiáng)度較大。

4.3.5 經(jīng)驗(yàn)變異函數(shù)的擬合

接下來(lái)，需要用變異函數(shù)模型總結(jié)空間自相關(guān)性。變異函數(shù)模型是距離和半方差的連續(xù)函數(shù)，該案例選擇球狀模型。一旦選定了模型形式，則需要調(diào)整模型參數(shù)以獲得經(jīng)驗(yàn)變異函數(shù)的最佳擬合。首先通過(guò)目視初步估計(jì)模型參數(shù)，然后通過(guò)自動(dòng)擬合對(duì)參數(shù)進(jìn)行調(diào)整：

上述代碼中，vgm函數(shù)用于指定一個(gè)變異函數(shù)模型。通過(guò)查看經(jīng)驗(yàn)變異函數(shù)圖圖3(a)，目視估計(jì)模型參數(shù)：偏基臺(tái)值psill為0.12；變程range為850；塊金值nugget為0.01；指定模型形式為球狀“Sph”。函數(shù)plot用于繪圖，結(jié)果如圖4(a)所示。

(a)目視擬合

由目視擬合結(jié)果可以看到，擬合曲線與經(jīng)驗(yàn)變異函數(shù)不太吻合，繼而使用fit.variogram函數(shù)進(jìn)行自動(dòng)調(diào)整：

擬合得到最終的變異函數(shù)模型圖4(b)表明：變程≈943m，偏基臺(tái)值≈0.012，基臺(tái)值≈0.111，塊基比≈0.11。

結(jié)語(yǔ)

“地統(tǒng)計(jì)學(xué)”可以廣泛應(yīng)用于解決多個(gè)領(lǐng)域的空間數(shù)據(jù)預(yù)測(cè)和插值問(wèn)題。變異函數(shù)是地統(tǒng)計(jì)學(xué)教學(xué)中的重點(diǎn)和難點(diǎn)。本文主要利用R語(yǔ)言的地統(tǒng)計(jì)學(xué)擴(kuò)展包gstat，以及隨sp包分發(fā)的Meuse數(shù)據(jù)集，展示了如何通過(guò)R語(yǔ)言進(jìn)行變異函數(shù)的理論和實(shí)踐教學(xué)。涉及以下步驟：計(jì)算點(diǎn)對(duì)半方差、繪制變異函數(shù)、確定點(diǎn)對(duì)分組的步長(zhǎng)，繪制經(jīng)驗(yàn)編譯函數(shù)圖并分析空間自相關(guān)性，最后擬合得到理論變異函數(shù)。該案例表明，R語(yǔ)言完全能夠?qū)崿F(xiàn)變異函數(shù)分析所需的功能，而且與ArcGIS這樣的商業(yè)化圖形界面軟件相比，R語(yǔ)言具有更強(qiáng)的擴(kuò)展性和靈活性，更有利于學(xué)生系統(tǒng)掌握地統(tǒng)計(jì)學(xué)的理論知識(shí)和實(shí)踐技能。