氣溫數據的時空變異結構分析

鄭 興 文，舒 紅，胡 泓 達

（1.貴州電力設計研究地理信息中心，貴州 貴陽 550002；2.武漢大學測繪遙感信息工程國家重點實驗室，湖北 武漢 430079）

0 引言

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

本文依據積和模型理論給出了時空變異函數的兩種實現方法，給出數據時空平穩化的預處理方法。將空間和時間變異函數積和成時空變異函數并實現可視化效果，通過這種時空變異函數分析東北氣象數據的時空結構特性。

1 時空變異函數

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

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

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

采用一類積和式變異函數擬合時空地理數據的時空變異結構，協方差函數和變異函數如下：

其中，Cst為時空協方差，Cs為空間協方差，Ct為時間協方差，γst、γs、γt分別是對應的時間、空間、時空變異函數，而Cst（0，0）、Cs（0）、Ct（0）分別是對應的基臺值。這里假定γst（0，0）＝γs（0）＝γt（0）＝0，模型中的系數k1、k2、k3滿足k1＞0、k2≥0、k3≥0，并通過正定條件由下式決定，推導過程見文獻［8］：

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

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

2 時空變異函數的計算

依據時空變異函數的理論模型運用R語言設計時空變異函數，由純空間和純時間變異函數通過積和模式構造時空變異函數，本文運用兩種方案對時空數據進行純時間和純空間變異函數構造，并通過模擬k1值確定Cst（0，0）將純空間和純時間變異函數積和為時空模型。

時空數據純空間和純時間變異函數構造方案一：純空間變異函數先對相同測站的時間序列觀測數據求平均，然后對不同測站的觀測平均值進行變異函數計算。純時間變異函數采用相同時態不同測站觀測屬性數據求平均，然后對不同時態的空間平均屬性值進行變異函數計算。運用R中的變異函數自動擬合獲取時空和空間變異函數公式。時空數據純空間和純時間變異函數構造方案二：純空間變異函數采用對空間測站，不同時態的屬性值分別計算變異函數值，獲得變異函數塊金值、偏基臺值、變程等參數，然后對各時態計算的變異函數參數求平均，用塊金值、偏基臺值、變程的平均值生成最終空間變異函數。純時間變異函數采用對應時態，不同空間測站時間序列屬性數據分別計算變異函數，獲得塊金值、偏基臺值、變程，然后對各測站點屬性數據的變異函數參數求平均，用塊金值、偏基臺值、變程平均值生成時間變異函數。

3 實驗分析

3.1 區域及數據介紹

東北三省地處119.70°～132.97°E，38.90°～52.97°N，北靠西伯利亞東部，深受寒冷干燥的冬季風影響，屬寒帶大陸性季風氣候。年氣溫一般在－5～10℃，冬季較長。周圍環山，中部以平原為主，山脈減弱了鄂霍次克海氣流的直接侵入，朝鮮半島和日本列島阻擋了其與太平洋的直接接觸，降低了海洋調節氣候的作用。在全國綜合自然區劃中，東北區屬于寒溫帶、溫帶、暖溫帶的濕潤、半濕潤地區［10］。

本實驗數據通過中國氣象科學數據共享服務網獲取，共有觀測站點87個，由于黑龍江肇州、吉林羅子河等6個站點數據嚴重缺失，實驗只采用81個觀測站1960年1月－2011年12月的月平均氣溫數據。

3.2 數據預處理

時空變異函數構造的重要前提是假設時空變量滿足二階平穩。時空氣溫數據可以看作為空間站點上的時間序列，時間序列一般包括：周期項、趨勢項、隨機項［11］。因此氣溫變量的時間序列分解為：

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

為保持數據連續性和整體變化趨勢，時間序列分解后的趨勢項和隨機項保留，趨勢項不明顯，對數據平穩性影響不大。針對上述氣溫時間序列去除周期項后進行自相關法平穩性檢測（圖2），可見自回歸系數隨時間延遲階數增大呈類似周期性余弦衰減，自相關系數很快衰減為0，說明時間序列去周期項后平穩［12］。

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

圖5 月平均氣溫站點數據去空間趨勢后擬合面Fig.5 Fitting surface for monthly average temperature data removed spatial trending

先對各站時間序列進行平穩性處理，再對整個空間數據進行平穩處理，最終達到構造時空變異函數的前提條件，即隨機變量滿足時空上的二階平穩。

3.3 時空變異函數實現過程

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

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

對比發現方案一和方案二空間變異函數的參數值具有較大差異，時間變異函數基本相同，分析原因發現第二種方案在確定空間變異函數中對應不同時間擬合變異函數（一共擬合360個）時有8個存在異?，F象，異常值嚴重影響了平均值，因此先去除異常值再進行平均，獲得相應參數見表1，結果與方案一基本相同。但方案二所用時間是方案一所用時間的幾十倍，而且方案一針對一組數據進行擬合可自動選擇最優擬合模型，方案二獲得均值參數后要根據過程擬合中的模型選擇最后變異函數模型。

空間和時間變異函數構造完成后，根據積和模型構造時空變異函數和時空協方差函數。時空變異函數：r1＋r2－k1＊r1＊r2（r1：純空間變異函數，r2：純時間變異函數）；時空協方差函數：k1＊（Snugget＋Spsill）＊（Tnugget＋Tpsill）＋k2＊（Snugget＋Spsill）＋k3＊（Tnugget＋Tpsill）－（r1＋r2－k1＊r1＊r2）（Snugget、Spsill為空間變異函數塊金值和偏基臺值，Tnugget、Tpsill為時間變異函數塊金值和偏基臺值）。k1、k2、k3的值滿足式（6）要求，k1的值決定了時空變異函數基臺值的大小。假設：

比較兩種方法構建的時空變異函數：方法二中塊金值比方法一小，方法二基臺值、變程值比方法一大，說明方法二刻畫的時空結構變異性要大，表達的時空結構性強，這與下文中變異程度結果相吻合。

4 結果分析

根據兩種生成純空間和純時間變異函數方案分別構造時空變異函數和協方差函數（圖7、圖8），并獲取基臺值、塊金值等相關參數（表2）。

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

空間變異程度是由隨機性因素引起的空間異質性占系統總變異的比例，當變異程度在區間25%～75%時認為變量為中等相關性，變異程度越小表示空間結構性越好?？臻g變異程度推廣到時空變異函數表示為隨機性因素引起的時空異質性占系統總變異的比例，當變異程度小于25%認為時空變量具有強時空相關性，說明具有很好的時空結構，變異程度在25%～75%之間認為時空變量具有中等時空相關性，大于75%時空相關性弱。從表3可知方案二變異程度較方案一小但都在25%～75%之間，表示東北月平均氣溫具有中等的時空相關性，時空結構較好。用方案二擬合的變異函數表述的東北月平均氣溫時空相關性和時空結構性略強于第一種方案。

可從圖6和表1獲知方案一和方案二中純空間和純時間變異函數的變異程度分別是：方案一，純空間變異函數變異程度為0，純時間變異函數變異程度為57%；方案二，純空間變異函數變異程度為0，純時間變異函數變異程度為33%?？梢姈|北三省月平均氣溫在空間上具有較強相關性，空間結構性強，而時間上相關性為中等，時間結構性較好。但時空積和模型的變異程度方案一方案二分別為69.7%和63%，說明時空結構性為中等。時空結構性下降主要可能是時間結構性不強造成的。

無論空間變量還是時間變量都認為在變程內是相關的，超出變程變量不相關。由圖7、圖8的時空變異函數可知空間變程不隨時間延遲而變化，類似的時間變程也不隨空間滯后而變化，這反映了時空具有相對穩定的結構特性。通過空間變程和時間變程可以獲得東北三省月平均氣溫認為在空間48 km內、時間3個月內具有一定的相關性，超出此范圍則認為東北月平均氣溫不存在相關性。

5 結論

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

本文討論了時空隨機場基礎上的時空變異函數的積和模型，給出了時空隨機場變異函數的計算方法，實驗分析了所給出時空變異函數計算方法的可行性，為基于時空變異函數的各種時空估計模型的實現提供了計算基礎。

［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］ 王仁鐸，胡光道.線性地質統計學［M］.北京：地質出版社，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］ 李莎，舒紅，徐正全.東北三省月降水量的時空克里金插值研究［J］.水文，2011，31（3）：31－35.

［12］ 李莎，舒紅，徐正全.利用時空Kriging進行氣溫插值研究［J］.武漢大學學報（信息科學版），2012，37（2）：237－241.

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