胡丹桂, 舒 紅, 胡泓達

(1.武漢大學(xué) 測繪遙感信息工程國家重點實驗室， 武漢 430079；2.武漢職業(yè)技術(shù)學(xué)院 計算機學(xué)院， 武漢 430074)

?

時空CoKriging的變異函數(shù)建模

胡丹桂1，2, 舒 紅1*, 胡泓達1

(1.武漢大學(xué) 測繪遙感信息工程國家重點實驗室， 武漢 430079；2.武漢職業(yè)技術(shù)學(xué)院 計算機學(xué)院， 武漢 430074)

在對地觀測中，所研究的地學(xué)變量不僅具有時間、空間特征，還受其他變量的影響，采用多元時空相關(guān)數(shù)據(jù)，可以提高時空估值的精度.時空CoKriging是多元時空插值中一種常用的方法，建立時空變異函數(shù)和協(xié)變異函數(shù)是時空CoKriging插值的關(guān)鍵一步.以東北三省為試驗區(qū)，利用該地區(qū)氣象站觀測數(shù)據(jù)的月平均空氣相對濕度為主變量，引入同時間同位置的月平均空氣溫度作為協(xié)變量，對空氣相對濕度和空氣溫度進行時空變異函數(shù)和時空協(xié)變異函數(shù)建模.實驗結(jié)果表明，采用和度量時空模型的時空變異函數(shù)的隨機性空間結(jié)構(gòu)建模的實際擬合效果較好.

多元時空變量； 和度量模型； 變異函數(shù)； 空氣相對濕度

地統(tǒng)計學(xué)是以區(qū)域化變量理論為基礎(chǔ)，以變異函數(shù)為工具，研究那些在空間分布上既有隨機性又有結(jié)構(gòu)性的自然現(xiàn)象的科學(xué).在自然科學(xué)及社會科學(xué)領(lǐng)域，有些變量不僅具有空間特征，而且具有時間特征，這時要把所研究的變量看成是時空隨機函數(shù).時空插值方法已成為解決時空離散點到連續(xù)體上的一種必不可少的映射方法，時空克里金是時空插值法中常用的一種方法[1].Marc G.Genton[2]，Matthew W. Mitchell[3]等對可分離時空協(xié)方差函數(shù)進行了理論研究和測試，對其優(yōu)缺點進行了深入的探討；Cressie and Huang[4]，Chunsheng Ma[5]，E. Porcu[6]，Gneiting[7]等提出了不可分離的，平穩(wěn)的協(xié)方差函數(shù)，允許時空交互；徐愛萍[8]，L.De Cesare[9]，D.E.Myers[10]，S. De Iaco[11]等將積和模型應(yīng)用于時空變異函數(shù)模型的建立中.他們雖然考慮到了時間空間特征，但僅局限于單變量.相對照，傳統(tǒng)的多元統(tǒng)計分析方法雖然考慮了多變量，卻很少考慮到多元信息的空間特征.但是，在諸如水文、石油、土壤、農(nóng)林、大氣和環(huán)保等科學(xué)研究中，所研究的變量不僅具有時間、空間特征，而且還在時間空間域中受其他相關(guān)變量的影響.將克里金法應(yīng)用于時空多元變量的插值研究中，一方面需要將克里金法進行時空擴展，另一方面將單變量插值的克里金擴展到多變量協(xié)同克里金.建立有效的多元時空協(xié)方差變異函數(shù)模型是研究時空協(xié)同克里金插值的關(guān)鍵一步.因此，本文將傳統(tǒng)地統(tǒng)計學(xué)擴展到多元時空地統(tǒng)計學(xué)，構(gòu)建多元時空變異函數(shù)及協(xié)變異函數(shù).

以東北三省72個氣象站點1996年～2005年的月均空氣相對濕度數(shù)據(jù)為例，以同時間同位置的月均空氣溫度為協(xié)變量，為時空協(xié)同克里金插值建立時空變異函數(shù)和時空協(xié)變異函數(shù).采用了和度量變異函數(shù)模型來擬合變量的時空變異結(jié)構(gòu)，將氣象要素模擬從空間維擴展至?xí)r空維.同時，考慮了氣溫數(shù)據(jù)對相對濕度的影響，在擬合相對濕度變異函數(shù)的過程中，加入了氣溫作為協(xié)變量，將時空克里金插值擴展到了時空協(xié)同克里金插值.

1 研究區(qū)域介紹及數(shù)據(jù)處理

1.1 區(qū)域和數(shù)據(jù)介紹

實驗研究的是我國東北三省(黑龍江、吉林、遼寧)，站點地處38.90°～52.97°N，119.70°～132.97°E之間.地面觀測臺站所觀測的資料來自1996年1月～2005年12月東北三省月空氣相對濕度.本實驗數(shù)據(jù)通過中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)獲取，共有觀測站點87個，由于部分站點數(shù)據(jù)嚴(yán)重缺失，實驗只采用72個觀測站數(shù)據(jù)，站點分布如圖1.實驗數(shù)據(jù)為1996年1月～2005年12月的月平均相對濕度和月平均氣溫值.

由于空氣相對濕度與空氣溫度相關(guān)性較好，我國東北地區(qū)，冬季干燥、夏季濕潤，空氣溫度對空氣濕度有著一定程度的相關(guān)影響，故選取空氣溫度為協(xié)變量.在實驗中，對空氣濕度和空氣溫度分別進行時空變異函數(shù)及兩者的協(xié)變異函數(shù)擬合.

圖1 東北三省觀測站點分布圖Fig.1 Monitoring station distribution in the three provinces of Northeast China

1.2 數(shù)據(jù)預(yù)處理

時空變異函數(shù)構(gòu)造的重要前提是假設(shè)時空變量滿足二階平穩(wěn).時空空氣相對濕度、空氣溫度數(shù)據(jù)可以看作為空間站點上的時間序列，時間序列一般包括：周期項、趨勢項、隨機項.因此空氣相對濕度變量的時間序列分解為：

(1)

為保持?jǐn)?shù)據(jù)連續(xù)性和整體趨勢變化，時間序列分解后的趨勢項和隨機項保留，因為趨勢項不明顯，對數(shù)據(jù)平穩(wěn)性影響不大.采用時序分解中的加法模型，將變量的周期項提取出來，余下的部分用于時空插值實驗.用自相關(guān)圖檢法[12]判斷去周期數(shù)據(jù)的平穩(wěn)性，如圖 4圖 5所示，空氣相對濕度和空氣溫度他們的時間序列的自相關(guān)系數(shù)(ACF)均隨延遲時期數(shù)很快衰減到±0.1以內(nèi)，圖示表明去周期項后數(shù)據(jù)近似平穩(wěn).

圖2 塔河空氣濕度時間序列綜合分析圖Fig.2 Air humidity time series decomposition in Tahe

圖3 塔河空氣溫度時間序列綜合分析圖Fig.3 Air temperature time series decomposition in Tahe

圖4 塔河空氣相對濕度時間序列隨機項自相關(guān)法檢測平穩(wěn)性Fig.4 The test of stationarity using autocorrelation method for Tahe Random

圖5 塔河空氣溫度時間序列隨機項自相關(guān)法檢測平穩(wěn)性Fig.5 The test of stationarity using autocorrelation method for Tahe random

空氣濕度和溫度變量在空間上也要滿足二階平穩(wěn)，將每一測站上的時間序列經(jīng)過去周期項后計算移動平均值，根據(jù)移動平均值探索空間空氣濕度和溫度的變化趨勢.用空間散點圖展示空間分布趨勢，圖 6是對濕度空間站點數(shù)據(jù)趨勢面的擬合.通過實驗對比，發(fā)現(xiàn)不管是濕度還是溫度均采用3次多項式(2)進行擬合，具有較好的擬合效果.

(2)

式中，value為月均空氣濕度，x，y為空間坐標(biāo)信息，a，b，c，d，e，f，g為3次多項式通過擬合得到的系數(shù).去除趨勢后空間數(shù)據(jù)分布呈現(xiàn)平穩(wěn)性，圖 7是月平均空氣濕度站點數(shù)據(jù)去空間趨勢后擬合面.分別對空氣濕度和空氣溫度各站的時間序列進行平穩(wěn)性處理再對整個空間數(shù)據(jù)進行平穩(wěn)處理，最終滿足構(gòu)造時空變異函數(shù)的前提條件，即隨機變量近似時空上的二階平穩(wěn).

圖6 月平均濕度空間站點數(shù)據(jù)趨勢面擬合Fig.6 Fitting trend surface for monthly average humidity data

圖7 月平均濕度站點數(shù)據(jù)去空間趨勢后擬合面Fig.7 Fitting surface for monthly average humidity data removed spatial trending

2 時空CoKriging變異函數(shù)建模方法

(3)

時空協(xié)同克里金(CoKriging)法的插值公式：

(4)

2.1 時空直接變異函數(shù)建模

假設(shè)Z={Z(s，t)，s、t∈Rd+1}(d+1表示歐式空間維加時間維)表示為時空隨機場，設(shè)定時空隨機場中兩位置的時空距離h=(hs，ht)，hs為矢量，代表樣本空間距離同時也包含方向信息，ht為樣本時間距離.當(dāng)Z(s，t)滿足二階平穩(wěn)時可定義其協(xié)方差函數(shù)為[9]：

C(hs，ht)=Cov(Z(s+hs，t+ht)，Z(s，t)).

(5)

顯然，協(xié)方差函數(shù)只與距離有關(guān)，與空間和時間位置無關(guān).在地統(tǒng)計學(xué)中，隨機變量Z(s，t)的期望不變(時空二階平穩(wěn)性)，協(xié)方差矩陣具有對稱正定性.理論上，滿足下列要求的連續(xù)函數(shù)可以定義為有效變異函數(shù)[9]：

σ2-C(hs，ht),

(6)式中,σ2為Z(s，t)的方差，正定條件即任意ai∈R，i=1，…，n，和任意正整數(shù)n，C必須符合：

(7)

1)可分離時空協(xié)方差模型

地統(tǒng)計學(xué)中協(xié)方差函數(shù)模型(如球狀模型、指數(shù)模型和高斯模型等)不能直接用于時空隨機變量的統(tǒng)計分析，必需進行時空擴展.如果一個隨機場Z的時空協(xié)方差函數(shù)可完全表示成純空間和純時間協(xié)方差函數(shù)相乘，則該時空協(xié)方差函數(shù)被認(rèn)為是可分離的[14].Mitchell Genton， Gumpertz(2005)[3]，Montserrat Fuentes(2004)[15]等文獻中提出了可分離型模型.這類模型相對容易構(gòu)建，且計算機實現(xiàn)插值效率較高，但往往要求一些比較理想的假設(shè)[16-17]，損失了精細(xì)的時空結(jié)構(gòu)信息或丟失了時空交互信息.

可分離型時空協(xié)方差函數(shù)可表示為：

C(hs，ht)=Cs(hs)·Ct(ht),

(8)

式中,C(hs，ht)是變量的時空協(xié)方差函數(shù)，Cs(hs)、Ct(ht)分別是純空間協(xié)方差和純時間協(xié)方差函數(shù).Cs(hs)、Ct(ht)的具體形式有多種選擇，或者同為一種模型，比如高斯模型，或者為不同模型，比如Cs(hs)為球形模型，Ct(ht)為指數(shù)模型，或是其他滿足正定條件的函數(shù).

2)時空積和模型

另一類稱為不可分離型時空協(xié)方差函數(shù)，它是將已知的純空間協(xié)方差與純時間協(xié)方差函數(shù)通過加乘、混合、積分等變換得到.積和式變異函數(shù)來擬合時空地理數(shù)據(jù)的時空變異結(jié)構(gòu)，協(xié)方差函數(shù)和變異函數(shù)如下[9]：

C(hs，ht)=k1Cs(hs)Ct(ht)+

k2Cs(hs)+k3Ct(ht),

(9)

γ(hs，ht)=(k1Ct(0)+k2)γs(hs)+

(k1Cs(0)+k3)γt(ht)-k1γs(hs)γt(ht),

(10)

式中，C(hs，ht)為時空協(xié)方差，Cs(hs)為空間協(xié)方差，Ct(ht)為時間協(xié)方差，γ(hs，ht)、γs(hs)、γt(ht)分別是對應(yīng)的時空、空間、時間變異函數(shù)，而C(0，0)、Cs(0)、Ct(0)分別是對應(yīng)的基臺值.這里，假定γ(0，0)=γs(0)=γt(0)=0，滿足k1>0、k2≥0、k3≥0，并通過正定條件由下式?jīng)Q定，推導(dǎo)過程見文獻[18]：

(11)

3)時空和度量模型

時空和度量模型也是不可分離時空變異函數(shù)模型.采用和度量模型變異函數(shù)來擬合時空地理數(shù)據(jù)的時空變異結(jié)構(gòu)，協(xié)方差函數(shù)可以表示為[19]：

(12)

式中，C(hs，ht)是空間間隔距離為hs，時間間隔距離為ht的協(xié)方差函數(shù)；

時空和度量變異函數(shù)為[19]：

(13)式中，γ(hs，ht)、γs(hs)、γt(ht) 分別為時空變異函數(shù)、純空間變異函數(shù)、純時間變異函數(shù).這個模型的缺點是求出這10個參數(shù)的值存在一定的復(fù)雜性.

2.2 時空協(xié)變異函數(shù)建模

在使用時空CoKriging來研究變量的時空變異性時，關(guān)鍵的一步就是決定變量之間的時空協(xié)變異函數(shù).它用來描述兩變量之間交叉的時空連續(xù)性.協(xié)變異函數(shù)定義為：

(14)

(15)

即在同一個位置上，將兩個變量的樣本數(shù)值相加，所得之和即是新變量在該點的樣本值，然后計算新變量的變異函數(shù)：

(16)

而這一新變量的變異函數(shù)與原變量的變異函數(shù)和協(xié)變異函數(shù)有如下關(guān)系：

(17)

因此可得：

(18)

3 時空CoKriging變異函數(shù)的擬合結(jié)果

分別采用時空可分離模型、時空積和模型、時空和度量模型3種不同的時空結(jié)合方法，對兩個時空變量空氣濕度和空氣溫度，采用時間去周期空間去趨勢后的殘差，來擬合時空變異函數(shù).如圖 8是空氣相對濕度的殘差的時空經(jīng)驗變異函數(shù)和理論變異函數(shù).圖 9是空氣溫度的殘差的時空經(jīng)驗變異函數(shù)和理論變異函數(shù).

圖8 空氣相對濕度殘差的時空變異函數(shù)Fig.8 Empirical and fitted space-time variograms of the humidity residuals

圖9 空氣溫度殘差的時空變異函數(shù)Fig.9 Empirical and fitted space-time variograms of the temperature residuals

并采用均方根差RMSD(root-mean-squared-difference)指標(biāo)來判斷一個經(jīng)驗變異函數(shù)擬合的好壞.表1分別是3種不同方法得到的模型的RMSD值，從表 1可以看出，針對這兩個時空變量，和度量時空模型的擬合效果最好.本實驗中，采用和度量時空模型來建立空氣濕度的殘差和空氣溫度的殘差的時空變異函數(shù).

圖10是將空氣濕度的殘差和溫度的殘差之和作為一個新變量，這個新變量的經(jīng)驗變異函數(shù)和經(jīng)和度量時空模型擬合的理論變異函數(shù).

表1 經(jīng)驗變異函數(shù)與理論變異函數(shù)之間的RMSD值Tab.1 RMSD between the empirical and fitted variograms

圖10 空氣濕度殘差與溫度殘差之和的時空經(jīng)驗變異函數(shù)和理論變異函數(shù)Fig.10 Empirical and fitted space-time variograms of the sum of humidity residuals and temperature residuals

使用和度量模型中的時空變異函數(shù)模型擬合數(shù)據(jù)，其中球狀模型被選為空間和時間變異函數(shù)，等式右邊第3項度量時空變異模型為指數(shù)模型[19].這3個組件都具有塊金(C0)，偏基(C1)和變程(r)的參數(shù)，度量時空變異模型還有各向異性比這個參數(shù)，一共10個參數(shù).表 2是空氣相對濕度的殘差，使用和度量模型擬合時空變異函數(shù)的10個不同參數(shù)的值.表 3是空氣溫度的殘差，使用和度量模型擬合時空變異函數(shù)的10個不同參數(shù)的值.表 4是空氣濕度殘差與溫度殘差之和，使用和度量模型擬合時空變異函數(shù)的10個不同參數(shù)的值.

表2 相對濕度殘差和度量變異函數(shù)模型擬合參數(shù)Tab.2 Fitting parameters of the relative humidity residuals variogram using sum metric model

經(jīng)過和度量時空模型擬合后，空氣濕度殘差的理論變異函數(shù)γ11為：

γ11(hs，ht)=1.01+7.41×

(19) 表3 空氣溫度殘差和度量變異函數(shù)模型擬合參數(shù)Tab.3 Fitting parameters of the air temperature residuals variogram using sum metric model

經(jīng)過和度量時空模型擬合后，空氣溫度殘差的理論變異函數(shù)γ22為：

(20) 表4 空氣濕度殘差與溫度殘差之和的 和度量變異函數(shù)模型擬合參數(shù)Tab.4 Fitting parameters of the sum of the relative humidity residuals and air temperature residuals variogram using sum metric model

(21)

所以，根據(jù)式(18)，空氣濕度的殘差和空氣溫度的殘差時空協(xié)變異函數(shù)γ12為：

(22)

4 結(jié)論

有效的變異函數(shù)模型是克里金插值的基礎(chǔ)，本文的重點對時空協(xié)同克里金的變異函數(shù)進行建模分析.在擬合變異函數(shù)之前，首先對變量進行時間上去周期，空間上去趨勢處理，以保證時空變量的平穩(wěn)性.得到插值結(jié)果后，必須先將其加上之前去除的周期項和趨勢項，才是最終估計結(jié)果.結(jié)果表明，空氣相對濕度和空氣溫度在純空間和純時間上均符合球狀模型，并且分別用3種不同的方法對這兩個變量進行時空變異函數(shù)建模，發(fā)現(xiàn)在可分離模型、積和模型、和度量模型這3種模型中，和度量模型的擬合效果最好.完成時空直接變異函數(shù)擬合后，最后根據(jù)D.E.Myers[21]提出的方法，進行時空協(xié)變異函數(shù)建模.

本文研究了多元時空數(shù)據(jù)進行時空直接變異函數(shù)和時空協(xié)變異函數(shù)的建模，不僅考慮了時間空間信息，而且還考慮了其他協(xié)變量.時空變異函數(shù)是時空克里金插值模型的權(quán)重參數(shù)構(gòu)建的基礎(chǔ)，對后續(xù)的時空協(xié)同克里金插值精度影響非常大.后續(xù)研究中，將研究大樣本下的時空協(xié)同克里金插值結(jié)果(估計量)的統(tǒng)計性質(zhì)(無偏、最優(yōu)、相合性、線性插值結(jié)果的漸近正態(tài)性)，發(fā)展適應(yīng)性多元時空統(tǒng)計模型應(yīng)用于多個地學(xué)領(lǐng)域時空數(shù)據(jù)分析并評價插值精度.

[1] 李 莎， 舒 紅， 徐正全. 利用時空Kriging進行氣溫插值研究[J].武漢大學(xué)學(xué)報:信息科學(xué)版， 2012， 37(2)：237-241.

[2]GentonMG.Separableapproximationsofspace-timecovariancematrices[J].Environmetrics， 2007， 18: 681-695.

[3]MitchellMW，GentonMG，GumpertzML.TestingforSeparabilityofspace-timecovariances[J].Environmetrics， 2005， 16: 819-831.

[4]CressieN，HuangHC.Classesofnonseparablespatio-temporalstationarycovariancefunctions[J].JournaloftheAmericanStatisticalAssociation， 1999， 94:1330-1340．

[5]MaCS.Spatio-temporalstationarycovariancemodels[J].JournalofMultivariateAnalysis， 2003， 86: 97-107.

[6]PorcuE，GregoriP，MateuJ.Nonseparablestationaryanisotropicspace-timecovariancefunctions[J].EnvironmentalResearchandRiskAssessment， 2006， 21: 113-122.

[7]GneitingT.Nonseparable，stationarycovariancefunctionsforspace-timedata[J].JournaloftheAmericanStatisticalAssociation， 2002， 97: 590-600．

[8] 徐愛萍， 圣文順， 舒 紅. 時空積和模型的數(shù)據(jù)插值與交叉驗證[J]. 武漢大學(xué)學(xué)報：信息科學(xué)版， 2012， 37(7)：766-770.

[9]DeCesareL，MyersDE，PosaD.Product-Sumcovarianceforspace-timemodeling:anenvironmentalapplication[J].Environmetrics， 2001， 12: 11-23．

[10]Myers.DE.Space-timecorrelationmodelsandcontaminantplumes[J].Environmetrics， 2002， 13:535-553.

[11]DeIacoS，MyersDE，PosaD.Space-timeanalysisusingageneralproduct-summodel[J].Statistics&ProbabilityLetters， 2001， 52 (1):21-28.

[12] 王 燕.應(yīng)用時間序列分析 [M].第2版.北京：中國人民大學(xué)出版社， 2008: 111-130.

[13]DeIacoS，PalmaM，PosaD.Modelingandpredictionofmultivariatespace-timerandomfields[J].ComputationalStatistics&DataAnalysis， 2005， 48: 525-547.

[14]MateuJ，PorcuE，GregoriP.Recentadvancestomodelanisotropicspace-timedata[J].StatisticalMethodsandApplications， 2008， 17: 209-23.

[15]FuentesM.Testingforseparabilityofspatial-temporalcovariancefunctions[J].JournalofStatisticalPlanningandInference， 2006， 136: 447-466.

[16]SteinM.Space-timecovariancefunctions[J].JournaloftheAmericanStatisticalAssociation， 2005， 100:310-321．

[17]FuentesM，ChenL，DavisJM.Aclassofnonseparableandnonstationaryspatialtemporalcovariancefunctions[J].Environmetrics， 2008， 19: 487-507.

[18]DeCesareL，MyersDE，PosaD.Estimatingandmodelingspace-timecorrelationstructures[J].Statistics&ProbabilityLetters， 2001， 51(1):9-14．

[19]DanielA.GriffithGB,HeuvelinkM.Derivingspace-timevariogramsfromSpace-TimeAutoregressive(STAR)modelspecifications[J].TheInternationalArchivesofthePhotogrammetry，RemoteSensingandSpatialInformationSciences， 2009， 38(2):15-19.

[20]DimitrakopoulosR，LuoX.SpatiotemporalModeling:CovariancesandOrdinaryKrigingSystems[C]//DimitrakopoulosR.GeostatisticsfortheNextCentury.SpringerNetherlands:KluwerAcademicPublishers， 1994(6):88-93．

[21]MyersDE.MatrixformulationofCokriging[J].MathematicalGeology，1982， 14 (3):249-257.

Variogram modeling in space-time CoKriging

HU Dangui1，2, SHU Hong1, HU Hongda1

( 1.State Key Laboratory of Information Engineering in Surveying，Mapping and Remote Sensing， Wuhan University， Wuhan 430079;2. Institute of Computer， Wuhan Polytechnic， Wuhan 430074 )

The variables used in various environmental studies not only have the spatial temporal dependence， but also are affected by other variables. However， the multivariate space-time correlated data can be used to improve the accuracy of space-time prediction and interpolation. The space-time CoKriging is a commonly used method of multivariate space-time interpolation， with modeling a space-time direct variogram and a cross variogram as a key step. Therefore， taking the three provinces in Northeastern China as study area. this paper modeled the space-time direct variograms and cross variograms for air relative humidity and air temperature by using the monthly mean relative air humidity observed in meteorological stations as a main variable and the air temperature in the same position and at the same time as a covariate variable. Inverse distance weighted interpolation was used to get the current meteorological observation data. The experimental results show that the space-time variogram using sum-metric model is well fitted， compared to the deterministic structure modeling spatial distance inversely．

multivariate space-time variables; Sum-metric model; variograms;air relative humidity

2014-11-05.

國家自然科學(xué)基金項目(41171313)；蘇州市科技計劃項目( SYG201319)；地理空間信息工程國家測繪地理信息局重點實驗室開放研究基金項目 (201329 )；湖北省自然科學(xué)基金項目 (2014CFB725).

1000-1190(2015)04-0596-07

P208< class="emphasis_bold">文獻標(biāo)識碼： A

A

*通訊聯(lián)系人. E-mail: shu_hong@whu.edu.cn.