劉祖涵, 王莉莉, 詹新武, 張紅梅, 王文豐

(南昌工程學院 a.信息工程學院；b.理學院；c.水利與生態(tài)工程學院，江西 南昌 330099)

目前，地理統(tǒng)計技術已被廣泛地應用于空間數據分析中，涉及煤炭行業(yè)、地質、生物學、礦業(yè)、氣象、水文、土壤、環(huán)境、生態(tài)和流行病學等領域[1-8].地理統(tǒng)計技術發(fā)展之初的主要目的是，在某一區(qū)域內，以少數的樣本數據，依據數據的空間變異結構，估算出完整的空間分布情況.聯合克利金法(CoKriging插值法)能以不同的統(tǒng)計結構，利用樣本數較為充足的輔助隨機變數來推估取樣點較少的氣候水文資料，進而獲得較為準確的因子變化特征的空間分布情況.同一空間位置樣點的多個屬性之間的某個屬性的空間分布與其他屬性(例如氣象臺站的經度、緯度、海拔等)密切相關，且某些屬性不易獲得，而另一些屬性則易于獲取.如果兩種屬性空間相關，可以考慮選用聯合克利金法，把區(qū)域化變量的最佳估值方法從單一屬性發(fā)展到兩個以上的協(xié)同區(qū)域化屬性，但它在計算中要用到兩種屬性各自的半方差函數和交叉半方差函數，比較復雜，學生不易理解和掌握.為此，筆者在地理信息系統(tǒng)上機理論與實驗教學過程中，在全面分析地統(tǒng)計學基本原理的基礎上，設計了地統(tǒng)計學的CoKriging法，以塔里木河流域日平均降水的 Hurst指數H1與其他屬性(包括流域23個氣象臺站的經度、緯度、海拔)為實驗數據，在ArcGIS 10.2 軟件平臺，對指數H1進行了聯合克利金法的空間插值.

1 區(qū)域化變量理論

一事物若能以特定統(tǒng)計空間結構表示，則稱為區(qū)域化(Regionalized).若Z(x)定義為位置x的隨機量測值，則Z(x)稱為區(qū)域化變量(Regionalized variable)[9]，自然界的空間變量，如礦產分布、地形高程分布及降雨分布，均可視為區(qū)域化變量.區(qū)域化變量具有兩個特性，局部而言，點與點間呈不規(guī)則變化，可視為隨機變數；整體而言，可用某種統(tǒng)計結構來代表其平均結構[10].

隨機變數Z(x)與Z(x+h)隨著距離的增大而相關性降低，因此共變異數函數cov(h)也隨距離的增大而逐漸變小.通常相距很遠的Z(x)與Z(x+h)之間的空間相依性為0，若h大于某特定距離，Z(x)與Z(x+h)的相關性趨近于0，則此特定距離稱為影響范圍(Range)，可以代表相關性有無的分界距離.區(qū)域性變數的空間結構可用半變異圖(Semivariogram)表示.半變異數是一種表現區(qū)域化變量沿某一特定方向的不同位置變化率的值.一般常用實驗半變異圖(Experimental semivariogram)來進行空間結構性分析.

1.1 半變異圖理論

假設區(qū)域化變量符合定常性(Stationary)，即其平均值為一常數.沿某方向任取兩點樣本配對，兩點間的距離為h個單位，

則定義半變異數為：

(1)

其中，Z(xi)為第i點的區(qū)域變數值；Z(xi+h)為第i點相距h個間隔的區(qū)域變數值；Nh為取樣點的配對數.

在實踐中，樣本的長度是有限的.把有限實測樣本值構成的變異函數稱為實驗，根據不同的h值及所對應的γ(h)值可以繪出半變異圖.相距愈小的兩點，其半變異數愈小.隨著距離的增加，任意兩點間的空間相依性愈小，區(qū)域變數中的半變異數趨向一個穩(wěn)定的值，此穩(wěn)定值稱為總體半變異數(Sill)，而總體半變異數的最小h值稱為影響范圍(Range)，如圖1(a)所示.

理論上，當h=0 時，γ(0)應該為0.但實際應用中，常會出現當h趨近于0 時，γ(0)的值呈現不連續(xù)情況，以數學式表示為：

(2)

這種當距離趨近于0時半變異圖的不連續(xù)性，稱為塊金效應(Nugget effect)[11]，如圖1(b)所示，代表距離小于采樣間距內的變異性以及測量誤差.

(a) 理想情況

(b)有塊金效應的情況

1.2 常用的半變異圖模型

由于一般繪制的實驗半變異圖常呈現散亂的不規(guī)則點，無法發(fā)現一個較佳的趨勢以求出最佳的Range和Sill值，因此，常假設某一函數模型為實驗半變異圖的方程式，利用試誤法求出最佳擬合曲線，再由曲線求出相對應的半變異數及其影響范圍.

變異函數表征了在聯合克利金法插值中權值取決于變量的空間結構[12]，一般用變異曲線表示，常用的有球狀模型、高斯模型及指數模型.

(1)球形模型(Spherical model)

(3)

Sill=C0+C1， Range=a.

(2)指數模型(Exponential model)

(4)

Sill=C0+C1， Range=3a.

(3)高斯模型(Gaussian model)

(5)

(4)冪級數模型(Power model)

γ(h)=C0+C1hθθ<2,

(6)

no Sill,no Range.

2 聯合克利金法

聯合克利金法(CoKriging)的主要原理為，將兩個或兩個以上具有高度空間相關性的區(qū)域化隨機變數合并考慮，進行空間資料推估[13].分析過程大致與一般克利金法(Kriging)相似，著眼點在于利用取樣點較多的資料輔助推估取樣點較少的資料，降低其推估誤差，常用于主要變數取樣點較少時，使用與主要變數高度相關，且取樣點較多的次要變數的輔助推估.需符合三個條件.

(1)線性

估計值為觀測值的線性組合.

(7)

(2)不偏估

估計值的期望值等于觀測值的期望值.

(8)

將(7)式代入(8)式可得：

(9)

(3)優(yōu)化

估計值與觀測值差的變異數為最小，即

(10)

將(7)式代入(10)式，利用拉格朗日法引入拉格朗日參數μ1和μ2，則拉格朗日函數L為：

(11)

(12)

寫成矩陣形式為：

(13)

(14)

估計誤差(Estimation error)為：

(15)

文獻[15]利用高斯模型成功地分析了烏魯木齊河流域月平均降水數據的分布規(guī)律，并給該模型的參數賦予了明確的物理意義，認為該模型不僅能夠實現降水在時間和空間上的插值，而且能夠實現降水量和降水分布函數的相互轉換，因此本文用高斯模型對文獻[1]中的塔里木河流域51年的日降水量的Hurst指數H1進行聯合克利金插值.需要指出的是，Hurst指數H1一般介于0到1之間，它反映的是非線性時間序列統(tǒng)計特征量的標度不變性和表征序列的長記憶性.

3 實驗過程及結果

在地理信息系統(tǒng)實驗課上，以ArcGIS 10.2軟件為平臺，對塔里木河流域51年日降水量的Hurst指數H1進行CoKriging空間插值舉例說明.

第一步，點擊地統(tǒng)計向導.

ArcGIS的地統(tǒng)計分析有一個地統(tǒng)計向導，按照這個向導一步一步就可以實現CoKriging空間插值.點擊向導，截圖如圖2所示.

第二步，選擇輸入的數據和屬性，采用CoKriging方法選擇插值模型，如圖3所示，然后點Next按鈕即可.

圖2地統(tǒng)計向導截圖

圖3 選擇插值模型截圖

第三步，選擇插值的主要變量(降水量)及三個次要變量(23個氣象臺站的經度、緯度及海拔高度)，然后點擊Next按鈕，如圖4所示.

圖4選擇插值的主要變量截圖

第四步，確定高斯模型.

這一步對插值的精度影響很大.根據左上角的高斯函數云圖，在右邊模型欄內，選擇對其擬合較好的模型.通過調整搜尋角度和搜尋半徑，確定一個合適的模型，該模型的參數在左下角框內給出，如圖5所示，然后點擊Next按鈕.同樣，可以選擇默認設置.

第五步，搜尋鄰居的調整.

選擇默認設置，截圖如圖6所示.

圖5 確定高斯模型的截圖

圖6 搜尋鄰居的調整截圖

第六步，交叉驗證結果.

這一步給出了上述設置計算結果的交叉驗證值，通過它可以分析各種驗證的精度.如果誤差很大，不能通過檢驗，則需要重新設置，如圖7所示.

圖7獲得交叉驗證結果截圖

點擊Finish按鈕，會得到這個過程的一個總結，如圖8所示.例如數據預處理過程、采用的模型等.

得到插值的結果如圖9所示.

圖8 過程總結截圖

圖9 插值結果截圖

結果表明，聯合克利金法可以精確顯示塔里木河流域的降水變化的長記憶性空間分布規(guī)律.這種方法有利于學生在上機實驗過程中理解和掌握相關知識，明確在實際的工程項目中，研究要素不僅受大尺度因素的制約，而且還受小尺度因素的影響，具有空間分布的不確定性.

參考文獻：

[1] 劉祖涵.塔里木河流域氣候——水文過程的復雜性與非線性研究[D].上海：華東師范大學,2014.

[2] HENGL T,MINASNY B,GOULD M.A Geostatistical Analysis of Geostatistics[J].Scientometrics,2009,80(2):491-514.

[3] SRIVASTAVA RM.Geostatistics:A Toolkit for Data Analysis,Apatial Prediction and Risk Management in the Coal Industry[J].International Journal of Coal Geology,2013,112(2):2-13.

[4] ASSARI A,MOHAMMADI Z.Combined Use of Geostatistics and Multi-criteria Decision Analysis to Determine New Pumping Well Locations in the Gol-gohar Open Pit Mine,Iran[J].Mine Water & the Environment,2017,36(2):1-16.

[5] LIU M,LEI LP,LIU D,et al.Geostatistical Analysis of CH4Columns over Monsoon Asia Using Five Years of Goast Observations[J].Remote Sensing,2016,8(5):361.

[6] BACHIR H,SEMAR A,MAZARI A.Statistical and Geostatistical Analysis Related to Geographical Parameters for Spatial and Temporal Representation of Rainfall in Semi-arid Environments:the Case of Algeria[J].Arabian Journal of Geosciences,2016,9(7):1-12.

[7] LIU RM,XU F,YU WW,et al.Analysis of Field-scale Spatial Correlations and Variations of Soil Nutrients Using Geostatistics[J].Environmental Monitoring & Assessment,2016,188(2):126.

[8] ANDERSON F.Application of Multivariate Geostatistics in Environmental Epidemiology:Case Study from Houston Texas[J].Occupational Diseases & Environmental Medicine,2016,(4):110-115.

[9] HUANG D,WANG G.Stochastic Simulation of Regionalized Ground Motions Using Wavelet Packets and Cokriginganalysis[J].Earthquake Engineering & Structural Dynamics,2015,44(5):775-794.

[10] JOURNEL A G,HUI JBREGTS C J.Mining Geostatistics[D].New York:Academic Press,1978.

[11] YIN J,NG SH,NG KM.Kriging Metamodel with Modified Nugget-effect:The Heteroscedastic Variancecase[J].Computers & Industrial Engineering,2011,61(3):760-777.

[12] 季青,余明.基于協(xié)同克里格插值法的年均溫空間插值的參數選擇研究[J].首都師范大學學報:自然科學版,2010,31(4):81-87.

[13] NERINI D,MONESTIEZ P,MANTé C.Cokriging for Spatial Functional Data[J].Journal of Multivariate Analysis,2010,101(2):409-418.

[14] ZHANG C,LI W,TRAVIS D.Gaps-fill of SLC-off Landsat ETM+ Satellite Image Using a Geostatistical Approach[J].International Journal of Remote Sensing,2007,28(22):5103-5122.

[15] 張小詠,劉耕年,李永化,等.高斯函數參量法及其在山區(qū)降水計算中的應用[J].地理研究,2008,27(3):594-602.