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

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

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

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

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

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

1.1 半變異圖理論

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

則定義半變異數(shù)為：

(1)

其中，Z(xi)為第i點(diǎn)的區(qū)域變數(shù)值；Z(xi+h)為第i點(diǎn)相距h個(gè)間隔的區(qū)域變數(shù)值；Nh為取樣點(diǎn)的配對(duì)數(shù).

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

理論上，當(dāng)h=0 時(shí)，γ(0)應(yīng)該為0.但實(shí)際應(yīng)用中，常會(huì)出現(xiàn)當(dāng)h趨近于0 時(shí)，γ(0)的值呈現(xiàn)不連續(xù)情況，以數(shù)學(xué)式表示為：

(2)

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

(a) 理想情況

(b)有塊金效應(yīng)的情況

1.2 常用的半變異圖模型

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

變異函數(shù)表征了在聯(lián)合克利金法插值中權(quán)值取決于變量的空間結(jié)構(gòu)[12]，一般用變異曲線表示，常用的有球狀模型、高斯模型及指數(shù)模型.

(1)球形模型(Spherical model)

(3)

Sill=C0+C1， Range=a.

(2)指數(shù)模型(Exponential model)

(4)

Sill=C0+C1， Range=3a.

(3)高斯模型(Gaussian model)

(5)

(4)冪級(jí)數(shù)模型(Power model)

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

(6)

no Sill,no Range.

2 聯(lián)合克利金法

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

(1)線性

估計(jì)值為觀測值的線性組合.

(7)

(2)不偏估

估計(jì)值的期望值等于觀測值的期望值.

(8)

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

(9)

(3)優(yōu)化

估計(jì)值與觀測值差的變異數(shù)為最小，即

(10)

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

(11)

(12)

寫成矩陣形式為：

(13)

(14)

估計(jì)誤差(Estimation error)為：

(15)

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

3 實(shí)驗(yàn)過程及結(jié)果

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

第一步，點(diǎn)擊地統(tǒng)計(jì)向?qū)?

ArcGIS的地統(tǒng)計(jì)分析有一個(gè)地統(tǒng)計(jì)向?qū)?，按照這個(gè)向?qū)б徊揭徊骄涂梢詫?shí)現(xiàn)CoKriging空間插值.點(diǎn)擊向?qū)В貓D如圖2所示.

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

圖2地統(tǒng)計(jì)向?qū)Ы貓D

圖3 選擇插值模型截圖

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

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

第四步，確定高斯模型.

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

第五步，搜尋鄰居的調(diào)整.

選擇默認(rèn)設(shè)置，截圖如圖6所示.

圖5 確定高斯模型的截圖

圖6 搜尋鄰居的調(diào)整截圖

第六步，交叉驗(yàn)證結(jié)果.

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

圖7獲得交叉驗(yàn)證結(jié)果截圖

點(diǎn)擊Finish按鈕，會(huì)得到這個(gè)過程的一個(gè)總結(jié)，如圖8所示.例如數(shù)據(jù)預(yù)處理過程、采用的模型等.

得到插值的結(jié)果如圖9所示.

圖8 過程總結(jié)截圖

圖9 插值結(jié)果截圖

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

參考文獻(xiàn)：

[1] 劉祖涵.塔里木河流域氣候——水文過程的復(fù)雜性與非線性研究[D].上海：華東師范大學(xué),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é)同克里格插值法的年均溫空間插值的參數(shù)選擇研究[J].首都師范大學(xué)學(xué)報(bào):自然科學(xué)版,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] 張小詠,劉耕年,李永化,等.高斯函數(shù)參量法及其在山區(qū)降水計(jì)算中的應(yīng)用[J].地理研究,2008,27(3):594-602.