鄭 興 文，舒 紅，胡 泓 達(dá)

（1.貴州電力設(shè)計(jì)研究地理信息中心，貴州 貴陽(yáng) 550002；2.武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079）

0 引言

自然環(huán)境中大量的現(xiàn)象可被視為時(shí)空隨機(jī)場(chǎng)的實(shí)現(xiàn)［1，2］，可以建立基于隨機(jī)場(chǎng)的時(shí)空模型來(lái)分析這些現(xiàn)象或觀(guān)測(cè)數(shù)據(jù)。但是，目前地質(zhì)統(tǒng)計(jì)學(xué)中多是純空間模型，單純使用空間模型分析時(shí)空數(shù)據(jù)會(huì)造成時(shí)間維有用信息的丟失。變異函數(shù)（亦稱(chēng)為半變異函數(shù)或是變差函數(shù)）是地質(zhì)統(tǒng)計(jì)學(xué)中重要的組成部分，不僅因?yàn)樗窃S多地質(zhì)統(tǒng)計(jì)學(xué)計(jì)算（如估計(jì)方差、離散方差、正則化變量的變異函數(shù)計(jì)算等）的基礎(chǔ)，還由于它能反映和刻畫(huà)區(qū)域化變量的許多性質(zhì)［3］。實(shí)現(xiàn)空間變異函數(shù)模型到時(shí)空變異函數(shù)模型的擴(kuò)展，是將變異函數(shù)中的空間距離（hs，hs∈Rd，d是物理空間維）擴(kuò)展為時(shí)空距離（h，h∈Rd＋1，d＋1表示物理空間維加一個(gè)時(shí)間維）的過(guò)程。De Ilaoc等指出變異函數(shù)從空間到時(shí)空的擴(kuò)展并不是一個(gè)簡(jiǎn)單的過(guò)程［4］，要解決時(shí)空量綱不一致（或時(shí)間和空間“距離”單位的統(tǒng)一）問(wèn)題和確保時(shí)空變異函數(shù)模型（或時(shí)空協(xié)方差）的最優(yōu)擬合問(wèn)題等。針對(duì)這些問(wèn)題，一些學(xué)者對(duì)時(shí)空變異函數(shù)進(jìn)行了理論探討。根據(jù)時(shí)空的結(jié)合方式將時(shí)空模型分為時(shí)空可分離模型和時(shí)空不可分離模型，時(shí)空可分離模型［5，6］只是簡(jiǎn)單地將時(shí)間協(xié)方差和空間協(xié)方差相乘或相加作為時(shí)空協(xié)方差，這種模型實(shí)現(xiàn)簡(jiǎn)單，但損失了精細(xì)的時(shí)空結(jié)構(gòu)信息或丟失了時(shí)空交互信息，比如積模型中任意兩空間點(diǎn)對(duì)上的兩個(gè)時(shí)間序列交叉協(xié)方差函數(shù)相同，而實(shí)際上時(shí)空協(xié)方差函數(shù)應(yīng)該隨空間滯后距離呈正比：Cst（h1，ht）∝Cst（h2，ht）。通常，時(shí)空不可分離模型［7］使用Bochner理論和選擇合適的譜密度構(gòu)造，模型考慮了數(shù)據(jù)內(nèi)在的豐富的時(shí)空結(jié)構(gòu)信息，但構(gòu)造方法復(fù)雜，很難獲得模型。特別地，De Cesare等提出一種時(shí)空變異函數(shù)的積和模型［8，9］，其比一般不可分離模型更加靈活。

本文依據(jù)積和模型理論給出了時(shí)空變異函數(shù)的兩種實(shí)現(xiàn)方法，給出數(shù)據(jù)時(shí)空平穩(wěn)化的預(yù)處理方法。將空間和時(shí)間變異函數(shù)積和成時(shí)空變異函數(shù)并實(shí)現(xiàn)可視化效果，通過(guò)這種時(shí)空變異函數(shù)分析東北氣象數(shù)據(jù)的時(shí)空結(jié)構(gòu)特性。

1 時(shí)空變異函數(shù)

假設(shè)Z＝｛Z（s，t），s、t∈Rd＋1｝（d＋1表示歐式空間維加時(shí)間維且一般d≤3）表示為時(shí)空隨機(jī)場(chǎng)，設(shè)定時(shí)空隨機(jī)場(chǎng)中兩位置的時(shí)空距離h＝（hs，ht），hs為矢量，代表樣本空間距離同時(shí)也包含方向信息，ht為樣本時(shí)間“距離”。當(dāng)Z（s，t）滿(mǎn)足二階平穩(wěn)時(shí)可定義其協(xié)方差函數(shù)為：

顯然，協(xié)方差函數(shù)只與距離有關(guān)，與空間和時(shí)間位置無(wú)關(guān)。在隨機(jī)變量Z（s，t）的期望不變（時(shí)空平穩(wěn)性）且協(xié)方差矩陣［Cst］n×n正定的條件下，下列函數(shù)為有效變異函數(shù)且可表示為：

式中：σ2為Z（s，t）的方差，正定條件即任意ai∈?，i＝1，…，n，和任意正整數(shù)n，Cst必須符合：

采用一類(lèi)積和式變異函數(shù)擬合時(shí)空地理數(shù)據(jù)的時(shí)空變異結(jié)構(gòu)，協(xié)方差函數(shù)和變異函數(shù)如下：

其中，Cst為時(shí)空協(xié)方差，Cs為空間協(xié)方差，Ct為時(shí)間協(xié)方差，γst、γs、γt分別是對(duì)應(yīng)的時(shí)間、空間、時(shí)空變異函數(shù)，而Cst（0，0）、Cs（0）、Ct（0）分別是對(duì)應(yīng)的基臺(tái)值。這里假定γst（0，0）＝γs（0）＝γt（0）＝0，模型中的系數(shù)k1、k2、k3滿(mǎn)足k1＞0、k2≥0、k3≥0，并通過(guò)正定條件由下式?jīng)Q定，推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)［8］：

式（5）中γs（hs）、γt（ht）的估值通過(guò)γst（hs，0）、γst（0，ht）的估值獲得，兩者通過(guò)式（5）和上述假定條件存在如下關(guān)系：

其中時(shí)間變異函數(shù)的滯后向量為rt，空間容差為δs，距離集合的分組數(shù)為：

2 時(shí)空變異函數(shù)的計(jì)算

依據(jù)時(shí)空變異函數(shù)的理論模型運(yùn)用R語(yǔ)言設(shè)計(jì)時(shí)空變異函數(shù)，由純空間和純時(shí)間變異函數(shù)通過(guò)積和模式構(gòu)造時(shí)空變異函數(shù)，本文運(yùn)用兩種方案對(duì)時(shí)空數(shù)據(jù)進(jìn)行純時(shí)間和純空間變異函數(shù)構(gòu)造，并通過(guò)模擬k1值確定Cst（0，0）將純空間和純時(shí)間變異函數(shù)積和為時(shí)空模型。

時(shí)空數(shù)據(jù)純空間和純時(shí)間變異函數(shù)構(gòu)造方案一：純空間變異函數(shù)先對(duì)相同測(cè)站的時(shí)間序列觀(guān)測(cè)數(shù)據(jù)求平均，然后對(duì)不同測(cè)站的觀(guān)測(cè)平均值進(jìn)行變異函數(shù)計(jì)算。純時(shí)間變異函數(shù)采用相同時(shí)態(tài)不同測(cè)站觀(guān)測(cè)屬性數(shù)據(jù)求平均，然后對(duì)不同時(shí)態(tài)的空間平均屬性值進(jìn)行變異函數(shù)計(jì)算。運(yùn)用R中的變異函數(shù)自動(dòng)擬合獲取時(shí)空和空間變異函數(shù)公式。時(shí)空數(shù)據(jù)純空間和純時(shí)間變異函數(shù)構(gòu)造方案二：純空間變異函數(shù)采用對(duì)空間測(cè)站，不同時(shí)態(tài)的屬性值分別計(jì)算變異函數(shù)值，獲得變異函數(shù)塊金值、偏基臺(tái)值、變程等參數(shù)，然后對(duì)各時(shí)態(tài)計(jì)算的變異函數(shù)參數(shù)求平均，用塊金值、偏基臺(tái)值、變程的平均值生成最終空間變異函數(shù)。純時(shí)間變異函數(shù)采用對(duì)應(yīng)時(shí)態(tài)，不同空間測(cè)站時(shí)間序列屬性數(shù)據(jù)分別計(jì)算變異函數(shù)，獲得塊金值、偏基臺(tái)值、變程，然后對(duì)各測(cè)站點(diǎn)屬性數(shù)據(jù)的變異函數(shù)參數(shù)求平均，用塊金值、偏基臺(tái)值、變程平均值生成時(shí)間變異函數(shù)。

3 實(shí)驗(yàn)分析

3.1 區(qū)域及數(shù)據(jù)介紹

東北三省地處119.70°～132.97°E，38.90°～52.97°N，北靠西伯利亞?wèn)|部，深受寒冷干燥的冬季風(fēng)影響，屬寒帶大陸性季風(fēng)氣候。年氣溫一般在－5～10℃，冬季較長(zhǎng)。周?chē)h(huán)山，中部以平原為主，山脈減弱了鄂霍次克海氣流的直接侵入，朝鮮半島和日本列島阻擋了其與太平洋的直接接觸，降低了海洋調(diào)節(jié)氣候的作用。在全國(guó)綜合自然區(qū)劃中，東北區(qū)屬于寒溫帶、溫帶、暖溫帶的濕潤(rùn)、半濕潤(rùn)地區(qū)［10］。

本實(shí)驗(yàn)數(shù)據(jù)通過(guò)中國(guó)氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)獲取，共有觀(guān)測(cè)站點(diǎn)87個(gè)，由于黑龍江肇州、吉林羅子河等6個(gè)站點(diǎn)數(shù)據(jù)嚴(yán)重缺失，實(shí)驗(yàn)只采用81個(gè)觀(guān)測(cè)站1960年1月－2011年12月的月平均氣溫?cái)?shù)據(jù)。

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

時(shí)空變異函數(shù)構(gòu)造的重要前提是假設(shè)時(shí)空變量滿(mǎn)足二階平穩(wěn)。時(shí)空氣溫?cái)?shù)據(jù)可以看作為空間站點(diǎn)上的時(shí)間序列，時(shí)間序列一般包括：周期項(xiàng)、趨勢(shì)項(xiàng)、隨機(jī)項(xiàng)［11］。因此氣溫變量的時(shí)間序列分解為：

式中隨機(jī)項(xiàng)R（t）滿(mǎn)足二階平穩(wěn)，以大連站氣溫觀(guān)測(cè)時(shí)間序列分解為例（圖1），將觀(guān)測(cè)氣溫時(shí)間序列（Observed）分解為趨勢(shì)項(xiàng)（Trend）、季節(jié)項(xiàng)（Seasonal）、隨機(jī)項(xiàng)（Random），可以看出自1960年1月到1989年12月30年的月平均氣溫值存在明顯的周期性，周期為一年，即冬季氣溫低，夏季氣溫高，月平均氣溫年較差在30℃左右，具有明顯的季節(jié)效應(yīng)。時(shí)間序列存在一個(gè)總體升高的趨勢(shì)，這可能與全球溫室效應(yīng)氣溫升高有關(guān)。

為保持?jǐn)?shù)據(jù)連續(xù)性和整體變化趨勢(shì)，時(shí)間序列分解后的趨勢(shì)項(xiàng)和隨機(jī)項(xiàng)保留，趨勢(shì)項(xiàng)不明顯，對(duì)數(shù)據(jù)平穩(wěn)性影響不大。針對(duì)上述氣溫時(shí)間序列去除周期項(xiàng)后進(jìn)行自相關(guān)法平穩(wěn)性檢測(cè)（圖2），可見(jiàn)自回歸系數(shù)隨時(shí)間延遲階數(shù)增大呈類(lèi)似周期性余弦衰減，自相關(guān)系數(shù)很快衰減為0，說(shuō)明時(shí)間序列去周期項(xiàng)后平穩(wěn)［12］。

變量在空間上也要滿(mǎn)足二階平穩(wěn)，將每一測(cè)站的時(shí)間序列經(jīng)過(guò)去季節(jié)效應(yīng)后計(jì)算移動(dòng)平均值，據(jù)此探索空間氣溫變化趨勢(shì)。用空間剖面圖法分析空間站點(diǎn)在東西和南北方向的分布趨勢(shì)（圖3），月平均氣溫在x軸（東西）方向呈中間略低兩側(cè)略高，在y軸（南北）方向呈南高北低趨勢(shì)，隨維度變化明顯，南北最大較差為15℃左右。用空間散點(diǎn)圖展示空間分布趨勢(shì)（圖4），采用3次多項(xiàng)式value＝ax＋by＋c（x＊y）＋d（x2）＋e（y2）＋f（x3）＋g（y3）進(jìn)行擬合（x，y為坐標(biāo)信息，value為月均值溫度），趨勢(shì)面擬合適應(yīng)度R2檢測(cè)值為0.73，具有較好的擬合效果，去除趨勢(shì)后空間數(shù)據(jù)呈平穩(wěn)分布（圖5）。

圖5 月平均氣溫站點(diǎn)數(shù)據(jù)去空間趨勢(shì)后擬合面Fig.5 Fitting surface for monthly average temperature data removed spatial trending

先對(duì)各站時(shí)間序列進(jìn)行平穩(wěn)性處理，再對(duì)整個(gè)空間數(shù)據(jù)進(jìn)行平穩(wěn)處理，最終達(dá)到構(gòu)造時(shí)空變異函數(shù)的前提條件，即隨機(jī)變量滿(mǎn)足時(shí)空上的二階平穩(wěn)。

3.3 時(shí)空變異函數(shù)實(shí)現(xiàn)過(guò)程

據(jù)時(shí)空變異函數(shù)的實(shí)現(xiàn)方案一對(duì)時(shí)間和空間數(shù)據(jù)求均值，根據(jù)均值計(jì)算純空間和純時(shí)間變異函數(shù)（圖6），擬合變異函數(shù)采用R庫(kù)函數(shù)包automap中的autofitVariogram（formula，input＿data，mode＝c（"Sph"，"Exp"，"Gau"，"Ste"），kappa＝c（0.05，seq（0.2，2，0.1），5，10））來(lái)實(shí)現(xiàn)，formula本文采用"z～1"形式是表示為構(gòu)造普通Kriging法的變異函數(shù)，z代表空間或時(shí)間月溫度均值，其他形式表示運(yùn)用其他種Kriging（如"z～x＋y"表示泛Kriging等），mode是可供選擇的擬合模型，常見(jiàn)的有橢球模型（"Sph"）、指數(shù)模型（"Exp"）、高斯模型（"Gau"）和Ste Matern模型（"Ste"）等［13］，kappa是擬合平滑度參數(shù)。空間變異函數(shù)描述空間相關(guān)性隨空間滯后的變化情況，隨空間距離變大變異函數(shù)值升高后趨于平穩(wěn)表示變量之間相關(guān)性變小直到不相關(guān)。從圖6A中可見(jiàn)變程47 km內(nèi)溫度是相關(guān)的，變程以外溫度不相關(guān)。時(shí)間變異函數(shù)描述時(shí)間變量隨時(shí)間延遲的相關(guān)程度，時(shí)間延遲是一維的，這與空間變異函數(shù)不同，從圖6B可見(jiàn)變程2.8個(gè)月內(nèi)認(rèn)為時(shí)間變量相關(guān)，塊金值0.97反映了時(shí)間變量隨機(jī)性的大小。

空間和時(shí)間變異函數(shù)的構(gòu)造方案二是“先擬合，再對(duì)參數(shù)求均值”，空間和時(shí)間變異函數(shù)采用參數(shù)的均值構(gòu)造。最終獲得的空間變異函數(shù)和時(shí)間變異函數(shù)的偏基臺(tái)值（Psill）分別為30.37、1.60，塊金值（Nugget）分 別 為 0、0.80，變 程 （Range）分別 為2 035.15、3.25。

對(duì)比發(fā)現(xiàn)方案一和方案二空間變異函數(shù)的參數(shù)值具有較大差異，時(shí)間變異函數(shù)基本相同，分析原因發(fā)現(xiàn)第二種方案在確定空間變異函數(shù)中對(duì)應(yīng)不同時(shí)間擬合變異函數(shù)（一共擬合360個(gè)）時(shí)有8個(gè)存在異常現(xiàn)象，異常值嚴(yán)重影響了平均值，因此先去除異常值再進(jìn)行平均，獲得相應(yīng)參數(shù)見(jiàn)表1，結(jié)果與方案一基本相同。但方案二所用時(shí)間是方案一所用時(shí)間的幾十倍，而且方案一針對(duì)一組數(shù)據(jù)進(jìn)行擬合可自動(dòng)選擇最優(yōu)擬合模型，方案二獲得均值參數(shù)后要根據(jù)過(guò)程擬合中的模型選擇最后變異函數(shù)模型。

空間和時(shí)間變異函數(shù)構(gòu)造完成后，根據(jù)積和模型構(gòu)造時(shí)空變異函數(shù)和時(shí)空協(xié)方差函數(shù)。時(shí)空變異函數(shù)：r1＋r2－k1＊r1＊r2（r1：純空間變異函數(shù)，r2：純時(shí)間變異函數(shù)）；時(shí)空協(xié)方差函數(shù)：k1＊（Snugget＋Spsill）＊（Tnugget＋Tpsill）＋k2＊（Snugget＋Spsill）＋k3＊（Tnugget＋Tpsill）－（r1＋r2－k1＊r1＊r2）（Snugget、Spsill為空間變異函數(shù)塊金值和偏基臺(tái)值，Tnugget、Tpsill為時(shí)間變異函數(shù)塊金值和偏基臺(tái)值）。k1、k2、k3的值滿(mǎn)足式（6）要求，k1的值決定了時(shí)空變異函數(shù)基臺(tái)值的大小。假設(shè)：

比較兩種方法構(gòu)建的時(shí)空變異函數(shù)：方法二中塊金值比方法一小，方法二基臺(tái)值、變程值比方法一大，說(shuō)明方法二刻畫(huà)的時(shí)空結(jié)構(gòu)變異性要大，表達(dá)的時(shí)空結(jié)構(gòu)性強(qiáng)，這與下文中變異程度結(jié)果相吻合。

4 結(jié)果分析

根據(jù)兩種生成純空間和純時(shí)間變異函數(shù)方案分別構(gòu)造時(shí)空變異函數(shù)和協(xié)方差函數(shù)（圖7、圖8），并獲取基臺(tái)值、塊金值等相關(guān)參數(shù)（表2）。

表2 兩種方案時(shí)空變異函數(shù)參數(shù)Table 2 Parameters of spatial-temporal variation function computed by two methods

空間變異程度是由隨機(jī)性因素引起的空間異質(zhì)性占系統(tǒng)總變異的比例，當(dāng)變異程度在區(qū)間25%～75%時(shí)認(rèn)為變量為中等相關(guān)性，變異程度越小表示空間結(jié)構(gòu)性越好?？臻g變異程度推廣到時(shí)空變異函數(shù)表示為隨機(jī)性因素引起的時(shí)空異質(zhì)性占系統(tǒng)總變異的比例，當(dāng)變異程度小于25%認(rèn)為時(shí)空變量具有強(qiáng)時(shí)空相關(guān)性，說(shuō)明具有很好的時(shí)空結(jié)構(gòu)，變異程度在25%～75%之間認(rèn)為時(shí)空變量具有中等時(shí)空相關(guān)性，大于75%時(shí)空相關(guān)性弱。從表3可知方案二變異程度較方案一小但都在25%～75%之間，表示東北月平均氣溫具有中等的時(shí)空相關(guān)性，時(shí)空結(jié)構(gòu)較好。用方案二擬合的變異函數(shù)表述的東北月平均氣溫時(shí)空相關(guān)性和時(shí)空結(jié)構(gòu)性略強(qiáng)于第一種方案。

可從圖6和表1獲知方案一和方案二中純空間和純時(shí)間變異函數(shù)的變異程度分別是：方案一，純空間變異函數(shù)變異程度為0，純時(shí)間變異函數(shù)變異程度為57%；方案二，純空間變異函數(shù)變異程度為0，純時(shí)間變異函數(shù)變異程度為33%。可見(jiàn)東北三省月平均氣溫在空間上具有較強(qiáng)相關(guān)性，空間結(jié)構(gòu)性強(qiáng)，而時(shí)間上相關(guān)性為中等，時(shí)間結(jié)構(gòu)性較好。但時(shí)空積和模型的變異程度方案一方案二分別為69.7%和63%，說(shuō)明時(shí)空結(jié)構(gòu)性為中等。時(shí)空結(jié)構(gòu)性下降主要可能是時(shí)間結(jié)構(gòu)性不強(qiáng)造成的。

無(wú)論空間變量還是時(shí)間變量都認(rèn)為在變程內(nèi)是相關(guān)的，超出變程變量不相關(guān)。由圖7、圖8的時(shí)空變異函數(shù)可知空間變程不隨時(shí)間延遲而變化，類(lèi)似的時(shí)間變程也不隨空間滯后而變化，這反映了時(shí)空具有相對(duì)穩(wěn)定的結(jié)構(gòu)特性。通過(guò)空間變程和時(shí)間變程可以獲得東北三省月平均氣溫認(rèn)為在空間48 km內(nèi)、時(shí)間3個(gè)月內(nèi)具有一定的相關(guān)性，超出此范圍則認(rèn)為東北月平均氣溫不存在相關(guān)性。

5 結(jié)論

本文基于一種積和模型的時(shí)空變異函數(shù)理論研究其實(shí)現(xiàn)過(guò)程，通過(guò)兩種構(gòu)造時(shí)空數(shù)據(jù)純空間和純時(shí)間變異函數(shù)方案分析其對(duì)時(shí)空變異函數(shù)構(gòu)造的影響，發(fā)現(xiàn)第二種方案描述的東北月平均氣溫時(shí)空結(jié)構(gòu)性更好，而方案一實(shí)現(xiàn)更為簡(jiǎn)單。提出一種使時(shí)空數(shù)據(jù)平穩(wěn)的方法，思路是先對(duì)測(cè)站點(diǎn)上時(shí)間序列去周期性，再去除空間趨勢(shì)面，以保證時(shí)空數(shù)據(jù)的平穩(wěn)性。通過(guò)預(yù)處理的時(shí)空數(shù)據(jù)構(gòu)造空間和時(shí)間變異函數(shù)，最終依據(jù)積和模型擬合獲取時(shí)空變異函數(shù)并可視化。分析發(fā)現(xiàn)東北月平均氣溫的時(shí)空結(jié)構(gòu)性較好，空間結(jié)構(gòu)性好于時(shí)間結(jié)構(gòu)性，并且在空間48 km內(nèi)、時(shí)間3個(gè)月內(nèi)（不考慮周期因素）東北月平均氣溫是相關(guān)的，超出此范圍不相關(guān)。

本文討論了時(shí)空隨機(jī)場(chǎng)基礎(chǔ)上的時(shí)空變異函數(shù)的積和模型，給出了時(shí)空隨機(jī)場(chǎng)變異函數(shù)的計(jì)算方法，實(shí)驗(yàn)分析了所給出時(shí)空變異函數(shù)計(jì)算方法的可行性，為基于時(shí)空變異函數(shù)的各種時(shí)空估計(jì)模型的實(shí)現(xiàn)提供了計(jì)算基礎(chǔ)。

［1］ EYNON B P，SWITZER P.The variability of rainfall acidity［J］.Canadian Journal of Statistics，1983，11（1）：11－23.

［2］ NHU LE D，PETKAU A J.The variability of rainfall acidity revisited［J］.Canadian Journal of Statistics，1988，16（1）：15－38.

［3］ 王仁鐸，胡光道.線(xiàn)性地質(zhì)統(tǒng)計(jì)學(xué)［M］.北京：地質(zhì)出版社，1988：38－40

［4］ DE IACO S，MYERS D E，POSA D.Nonseparable space-time covariance models：Some parametric families［J］.Mathematical Geology，2002，34（1）：23－42.

［5］ DE CESARE L，MYERS D，POSA D.Spatial-Temproal Modeling of SO2in Milan District［M］.Geostatistic Wollongong′96.Dordrecht，1996.

［6］ ROUHANI S，HALL T J.Space-time kriging of groundwater data［A］.Geostatistics［C］.Springer Netherlands，1989.639－650.

［7］ CRESSIE N，HUANG H C.Classes of nonseparable，spatiotemporal stationary covariance functions［J］.Journal of the A-merican Statistical Association，1999，94（448）：1330－1339.

［8］ DE CESARE L，MYERS D，POSA D.Estimating and modeling space-time correlation structures［J］.Statistics ＆Probability Letters，2001，51（1）：9－14.

［9］ DE CESARE L，MYERS D E，POSA D.Product-sum covariance for space-time modeling：An environmental application［J］.Environmetrics，2001，12（1）：11－23.

［10］ 周琳.東北氣候［M］.北京：氣象出版社，1991.1－3.

［11］ 李莎，舒紅，徐正全.東北三省月降水量的時(shí)空克里金插值研究［J］.水文，2011，31（3）：31－35.

［12］ 李莎，舒紅，徐正全.利用時(shí)空Kriging進(jìn)行氣溫插值研究［J］.武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2012，37（2）：237－241.

［13］ 張仁鐸.空間變異理論及應(yīng)用［M］.北京：科學(xué)出版社，2005. g g g