“地統計學方法”課程實驗設計

楊 勇, 楊 雪, 劉詩詩

(華中農業大學 資源與環境學院, 湖北 武漢 430070)

地統計學(Geostatistics)是以區域化變量理論為基礎,以變異函數為主要工具,研究那些在空間分布上既有隨機性，又有結構性,或空間相關性和依賴性的自然現象的科學[1]。它起源于20世紀50年代的南非采礦業中,目前已被廣泛應用于地質[2-3]、土壤[4-5]、生態環境[6-7]、氣象[8-9]等涉及地理屬性空間/時空變異及插值的多個領域,是集空間數據采集、統計分析、建模與計算、地圖制圖等多個技術環節理論與方法的集合。其核心技術包括基于空間樣點的地理變量試驗變異函數技術、理論變異函數擬合和地統計插值等,最常用的方法為普通克里格方法。此外,隨著各個領域應用的延伸,各種衍生方法不斷推出,如協同克里格[10]、回歸克里格[11]、泛克里格[12]、指示克里格[13]、析取克里格[14]以及各種隨機模擬方法[15-16]等,極大地豐富了地統計學的方法體系和內涵。

目前,有多個軟件包含地統計功能,如ArcGIS的Geostatistics模塊、GS+、GSLIB、Surfer等,各個軟件均能執行部分地統計學方法。這些軟件降低了使用門檻的同時,也促進了地統計學方法的普及,使得用戶在不了解地統計理論和技術的情況下,也能利用地統計學方法對所分析的地理屬性進行空間插值或不確定性評估。但從已發表的文獻來看,由于不掌握地統計學的基本理論知識,造成一些失敗的建模和計算結果,具體表現在以下幾個方面:

(1) 變異模型精度不高,造成最終克里格插值結果精度較低,這可能是在計算試驗變異函數時,選擇了不合適的滯后間距,或是忽略了變異函數不同方向的各向異性等原因造成的；

(2) 只依賴樣點數據,而未能有效地利用其他輔助數據來提供插值精度；

(3) 未能根據數據特點,選擇合適的地統計學方法。

鑒于地統計學方法的廣泛應用及其在使用過程中容易出現的失誤,為普及地統計學基本理論,介紹地統計學方法體系,我校為地理信息科學和生態學專業的本科生和研究生開設了“地統計學方法”課程。由于該課程內容多、算法復雜、理論性較強,單獨依靠理論課教學手段難以使學生理解抽象的概念和復雜的算法流程,因此,開設了相關實驗教學環節來加深學生對理論和算法的理解并提高學生學習興趣。

1 地統計學方法實驗設計

地統計學的特點之一就是方法多樣,但目前尚沒有一個軟件囊括了所有的方法,且有的方法實施過程(如回歸克里格)需要借助多個軟件完成,因此課程的實驗過程也是多個軟件(還包括非GIS或地統計軟件)綜合利用的過程。另外,本課程的主要授課對象是地理信息科學本科生,更注重學生對基本理論和方法的理解,若實驗過程只依賴已發布的地統計軟件或模塊,這些理論和算法對學生來講就是暗箱操作,因此,在部分實驗環節,設計了程序設計環節,由學生利用基本的程序設計語言完成關鍵算法。

為了便于方法比較,本實驗課程使用一套數據作為實驗對象,數據為湖北省某縣土壤樣點數據,數據項包括:[x,y, 有機質(z1),土壤全氮(z2),pH值(z3)],其中前2項為樣點的坐標數據,后3項為待分析的地理屬性項。本數據共計658個樣點,另外還包括該縣的行政邊界地形數據等。根據地統計學授課內容、實驗數據和上述分析,設計了6個實驗內容。

1.1 數據的基本統計分析和預處理

通過數據的基本統計分析能夠了解待插值地理屬性的基本統計特征,了解其概率分布特點,為正式計算前是否要對數據做某些預處理提供參考依據。基本統計分析內容包括最大值、最小值、均值、中值、標準差、變異系數、峰度系數、偏度系數、K-S檢驗等。其中,聯合峰度系數、偏度系數、K-S檢驗的結果判斷數據是否接近正態分布,從而是否要對屬性數據進行轉換,以使轉換后的數據接近正態分布。數據預處理主要包括2個部分:一是待分析屬性數據的異常值檢驗及處理,鑒于實驗數據量大,推薦使用2倍標準差或3倍標準差法；二是數據轉換,即若上一步判斷原始數據不接近正態分布,則要對數據進行轉換,轉換方法包括對數、平方、反正弦轉換等。本實驗利用Excel和SPSS軟件聯合完成,要求學生以圖表的形式提供結果,并為下一個實驗提供接近正態分布的屬性數據。

1.2 基于樣點的試驗變異函數計算

定義z(s)為在位置s=(x,y)處的地理屬性值,即地理屬性Z在位置s處的觀測值,則試驗變異函數的計算公式為:

(1)

其中h為給定的空間距離,N(h)為樣點中,符合空間距離h的點對數量,z(si)和z(si+h)為空間距離為h的點對。取若干個等間隔的h,則能計算得到若干個γ(h)*,形成[h,γ(h)*]數據對,這就是基于樣點的試驗變異函數計算結果。當樣點為規則取樣時,只要以取樣間距為距離間距,則能很容易地獲取各距離h下的若干點對。但大部分的樣點并非規則布設,實際計算時若嚴格按給定的距離尋找樣點對,則出現樣點對極少甚至可能找不到的可能。因此指導學生按照以下步驟進行計算:

(1) 假設步長為lag,當前滯后級別為n(n為正整數),則h=n*lag；

(2) 對所有樣點,找到點對(si,sj),其符合條件:(n-0.5)*lag

(3) 計算[z(si)-z(sj)]2,記為Si；

(6) 將各個滯后距級別上的(havg,γ(havg)*)繪制成圖,形成試驗變異(半方差)圖。

此實驗要求學生根據上述算法,基于程序開發語言(如C#或Matlab)完成。

1.3 理論變異函數擬合

基于試驗變異函數散點圖,由學生根據圖形特點選擇一個合適的理論變異函數模型(如球狀、高斯、指數、線性模型等)后,利用最小二乘法估計模型中的參數。但由于大部分理論模型均不是線性模型,因此,在參數估計前,需進行模型變換,即將非線性模型變換成線性模型,待擬合好線性模型中的參數后,再反算得出理論模型中的參數解[1]。另外,若選擇套合模型作為理論模型,由于待估計的參數數量較多,可引導學生使用智能算法(如遺傳算法)進行模型參數的求解[17]。模型求解后,要求學生利用試驗變異函數散點的實際值及模型擬合值之間的相關系數和均方根誤差進行精度驗證,若精度較低,則重新進行參數估計,直到滿足給定的精度要求為止。此實驗要求學生根據算法,利用程序開發語言(如C#或Matlab)完成。

1.4 普通克里格算法實現及地圖制圖

普通克里格方法是目前運用最廣泛的克里格方法,也是目前各領域最常用的地理屬性空間插值算法,其核心思想是利用待估位置周圍的若干個已知樣點對地理屬性值進行線性無偏最優估計,并給出估計方差。利用這個思想,推導出克里格方程組:

(2)

(1) 設定待插值柵格矩陣,主要是確定柵格單元大小,為制圖美觀,柵格單元大小不應大于50 m。設定好柵格矩陣后,每個柵格單元的中心點坐標作為該柵格的待插值位置,其插值結果代表了整個柵格單元的克里格估值結果；

(2) 按一定順序遍歷每個柵格,尋找離柵格單元的中心點最近的n個樣點(一般為8個),構建克里格方程組,解該方程組,得到各樣點的權重系數,進而得到待估值位置的克里格插值結果；

(3) 對柵格矩陣中每個柵格估值完后,將結果存儲成ASC格式文件,導入到ArcGIS軟件中,成為柵格圖層。利用區域邊界對柵格圖層進行裁剪,只保留區域邊界之內的部分,再添加指北針、比例尺、圖例等地圖整飾要素,最終形成區域土壤屬性空間分布專題地圖。

1.5 指示克里格實驗

與其他克里格方法所得結果為地理屬性值不同,指示克里格方法所得結果為地理屬性超過或不超過某閾值概率的空間分布,常用于環境風險評價等領域。指示克里格方法步驟與普通克里格法大體相同,區別在于計算前,要將樣點數據按照某個閾值做指示變換,具體實驗步驟如下:

(1) 設定閾值c,將原始樣點的屬性數據做如下變換:

(3)

即將原始地理屬性含量值按照閾值c變換成指示值；

(2) 基于上一步獲得的指示值,計算試驗變異函數,并擬合理論變異函數模型；

(3) 構建區域柵格矩陣,基于指示值和理論變異函數模型,對各柵格中心位置進行普通克里格估值,得到每個柵格超過閾值c的概率；

(4) 利用區域邊界數據進行邊界裁剪,并添加指北針、比例尺、圖例等地圖整飾。

該實驗可在ArcGIS中的Geostatistics模塊上完成。

1.6 回歸克里格實驗

(1) 在ArcGIS中,提取樣點上各環境變量的值；

(2) 在SPSS軟件中,計算土壤有機質與各環境變量的相關系數；

(3) 在SPSS中,用逐步回歸法建立土壤有機質與各環境變量的線性回歸方程,即得到趨勢模型m(s)；

(4) 使用ArcGIS的地圖代數功能,利用上一步得到的趨勢模型,計算各位置上的趨勢值；

(5) 提取趨勢值到樣點圖層中,可在Excel中得到各樣點的殘差值；

(6) 利用ArcGIS中的Geostatistics模塊完成樣點殘差的克里格插值；

(7) 使用ArcGIS的地圖代數功能,將各位置上的趨勢值與殘差的克里格插值結果相加,得到回歸克里格預測結果；

(8) 進行邊界裁剪,添加地圖整飾,完成地圖制圖并提交。

2 總結與討論

綜上所述,地統計學實驗課程分為6個內容,其中第2、3、4個內容需要學生利用程序設計語言自行開發,通過這些實驗過程,能加深學生對地統計學最基本理論和最核心算法的理解和掌握,但也對學生提出了更高的要求,即需要有良好的程序設計基礎。第1、5、6個實驗需要學生綜合利用多個軟件完成,其中第5、6個實驗可以看作是對基本的統計學方法的擴展,拓寬學生視野的同時,鍛煉了學生多種工具軟件的綜合使用能力。而在第6個實驗的教學過程中,只介紹回歸克里格的思想,但不明確給出具體的實驗步驟和所使用的軟件,由學生自行探索,實驗完成后,提升了學生解決問題的能力。另外,第4、5、6個實驗中地圖制圖的內容能夠鞏固地理信息科學專業的學生地圖學的知識和技能。

此實驗方案已實施3年,并不斷地修改和完善,根據實驗后學生的反饋,通過實驗獲得的收獲包括:

(1) 通過對變異函數和普通克里格插值的算法實現,深刻地理解了地統計學的基本理論和算法,同時在空間數據結構、程序設計語言、數據和文件間的I/O等方面均有深刻體會和提升；

(2) 同歸指示克里格和回歸克里格實驗,印證了課堂的理論教學內容,拓展了視野；

(3) 通過地圖制圖,鞏固了地圖學中關于地圖整飾和地圖符號設計知識；

(4) 通過觀察所得土壤屬性空間分布結果,理解了地理屬性空間分異規律和特征。

但該實驗方案還存在一些問題,主要是第2、3、4內容難度太大,對一些程序設計基礎不好的學生來講,難以開展。這需要在實驗前對學生的技術能力進行摸底調查,如果大部分學生技術力量薄弱,那么可以通過實驗分組來解決。即將多個學生分成一個組,每個學生只負責其中的一部分,共同協調合作完成這部分實驗內容。