許家琦，李顏伶，舒 紅＊

（1.武漢大學(xué)測(cè)繪遙感信息工程國家重點(diǎn)實(shí)驗(yàn)室，武漢430079；2.青島勘察測(cè)繪研究院，山東青島266032）

在對(duì)空間數(shù)據(jù)的分析中總是假設(shè)空間協(xié)方差結(jié)構(gòu)平穩(wěn)，然而這一假設(shè)在實(shí)際中常常不成立.例如降水?dāng)?shù)據(jù)就顯示出空間非平穩(wěn)性［1－2］，月降水中硫酸濃度也呈現(xiàn)出了非平穩(wěn)結(jié)構(gòu)［3］.空間數(shù)據(jù)的平穩(wěn)性檢驗(yàn)有多種方法.地理加權(quán)回歸法將空間的概念引入到一般線性回歸的回歸參數(shù)中，回歸系數(shù)隨著空間位置的變化而變化［4］.地理加權(quán)回歸法檢驗(yàn)的是變量之間相關(guān)關(guān)系的空間異質(zhì)性，地理加權(quán)回歸法（GWR）在已知數(shù)據(jù)非平穩(wěn)的情況下分析這種非平穩(wěn)性，而不是用于數(shù)據(jù)檢驗(yàn)是否平穩(wěn).本文采用的3種方法（三維散點(diǎn)圖法，趨勢(shì)線法和鄰域圓法）都是用來檢驗(yàn)屬性數(shù)據(jù)本身的空間異質(zhì)性.譜方法用局部平穩(wěn)過程的卷積來表示整個(gè)非平穩(wěn)過程，并用這個(gè)非平穩(wěn)模型來進(jìn)行譜密度的參數(shù)和非參數(shù)估計(jì)［5－6］.由于每個(gè)被視為平穩(wěn)的局部區(qū)域中都要有足夠多的數(shù)據(jù)來計(jì)算空間變異函數(shù)然或者空間協(xié)方差函數(shù)［7］，因此這種方法并不適用于我國較為稀疏的氣象站點(diǎn)數(shù)據(jù).三維散點(diǎn)圖法擬合的是整個(gè)研究區(qū)域內(nèi)屬性數(shù)據(jù)的趨勢(shì)面；趨勢(shì)線法擬合的是投影面上所有屬性數(shù)據(jù)的趨勢(shì)；雖然鄰域圓法也是把空間分割為一些小的區(qū)域來研究，但鄰域圓法關(guān)注的是局部區(qū)域內(nèi)的二階矩.均值和方差的計(jì)算與經(jīng)驗(yàn)變異函數(shù)的計(jì)算不同，對(duì)于參與計(jì)算的氣象站點(diǎn)沒有太多數(shù)量上的要求.在非平穩(wěn)的空間協(xié)方差結(jié)構(gòu)分析中，把不平穩(wěn)各向異性的地理空間映射到平穩(wěn)各向同性的頻譜空間中進(jìn)行處理，然后將處理結(jié)果逆變換到地理空間中［8］.Sampson＆Guttorp法（SG模型）計(jì)算復(fù)雜，物理含義不夠直觀明確.而本文采用3種方法：（1）三維散點(diǎn)圖法有著簡(jiǎn)單直觀的優(yōu)點(diǎn)，通過查看氣象數(shù)據(jù)在三維空間中的數(shù)值大小來初步判斷數(shù)據(jù)是否平穩(wěn)；（2）在趨勢(shì)線法中，把氣象數(shù)據(jù)投影到東西和南北兩個(gè)面上，在這兩個(gè)面上分別用多項(xiàng)式法擬合出兩條趨勢(shì)線，從而判斷出空間數(shù)據(jù)在東西方向和南北方向上是否平穩(wěn)［9］.（3）由于隨機(jī)變量的一階矩和二階矩已經(jīng)足夠?yàn)榇蠖鄶?shù)實(shí)際問題提供可以接受的近似解決方法［10］，在鄰域圓法中通過統(tǒng)計(jì)圓內(nèi)樣本數(shù)據(jù)的均值和方差來檢驗(yàn)我國東北地區(qū)氣象數(shù)據(jù)的空間平穩(wěn)性［11］.

1 研究區(qū)域概況

本文采用我國東北地區(qū)的氣象觀測(cè)站點(diǎn)數(shù)據(jù)，包括2006年1月1日的日平均氣溫?cái)?shù)據(jù)，2006年1月1日的日平均風(fēng)速數(shù)據(jù)以及2006年5月的月降水?dāng)?shù)據(jù)作為實(shí)驗(yàn)數(shù)據(jù).氣象站點(diǎn)分布如圖1所示.

冬季，受極地大陸氣團(tuán)控制，東北區(qū)氣候嚴(yán)寒.1月是全年最冷月，全區(qū)的月平均氣溫均在－4℃以下.溫度梯度之大，為全年之冠.大小興安嶺北段的平均氣溫在－30℃左右，是我國冬季最冷的地方.冬季（以1月為代表）在蒙古冷高壓影響下，空氣穩(wěn)定，風(fēng)速較小.特別是松花江流域以北地區(qū)，因位置最北，常在高壓內(nèi)部控制之下，冬季為全年風(fēng)速最小季節(jié).1月份全區(qū)平均風(fēng)速在1～5 m／s之間；山區(qū)13 m／s；平原多在3～4 m／s.本區(qū)的降水主要集中于夏季.各地降水一般從5月份明顯增加，9月份又陡降.降水量自東南向西北有遞減的趨勢(shì).降水日多于20 d的區(qū)域由東北部的山區(qū)延伸到了松遼平原的東部，降水日數(shù)的高值區(qū)仍是小興安嶺南部和長白山一帶［12］.

圖1 中國東北地區(qū)氣象觀測(cè)站點(diǎn)分布圖Fig.1 Location map of meteorological stations in the Northeast area of China

2 研究方法

設(shè)一隨機(jī)函數(shù)Z，其空間分布律不因平移而改變，即若對(duì)任一向量h，兩個(gè)k維向量的隨機(jī)變量｛Z（x1），Z（x2），…，Z（xk）｝和｛Z（x1＋h），Z（x2＋h），…，Z（xk＋h）｝.

有相同的分布律.這種平穩(wěn)假設(shè)至少要求Z（x）的各階矩均存在，且平穩(wěn)，而在實(shí)際工作中卻很難滿足.在線性地質(zhì)統(tǒng)計(jì)學(xué)研究中，只需假設(shè)其1、2階矩存在且平穩(wěn)就夠了，因而提出二階平穩(wěn)假設(shè).即

假定：1）在整個(gè)研究區(qū)內(nèi)，區(qū)域化變量Z（x）的期望存在且等于常數(shù)：

2）在整個(gè)研究區(qū)內(nèi)，區(qū)域化變量的空間協(xié)方差函數(shù)存在且平穩(wěn)［13］：

2.1 三維散點(diǎn)圖法

空間現(xiàn)象可以分解為兩部分：區(qū)域性趨勢(shì)（平均值）和殘差（變化性）.二階平穩(wěn)意味著區(qū)域范圍內(nèi)隨機(jī)變量的平均值和方差均沒有趨勢(shì)變化.平穩(wěn)性就是觀測(cè)值圍繞固定常數(shù)上下波動(dòng)，且離開均值上下波動(dòng)的幅度大體相似（方差為常量）.反之，如果趨勢(shì)面不是一個(gè)水平面或者數(shù)據(jù)圍繞趨勢(shì)面波動(dòng)的幅度不均衡就是非平穩(wěn).所以可以對(duì)觀測(cè)值作散點(diǎn)圖，同時(shí)顯示趨勢(shì)面，通過人眼初步觀察一下是否滿足平穩(wěn)性.

2.2 趨勢(shì)線法

將三維散點(diǎn)圖分別向東西方向和南北方向上投影，繪制氣象數(shù)據(jù)與x軸，y軸的二維散點(diǎn)圖.用一條紅色的趨勢(shì)線（既可以是直線也可以使曲線）來擬合這些數(shù)據(jù)，通過兩個(gè)方向上的剖面圖可以清楚看到數(shù)據(jù)趨勢(shì)在哪個(gè)方向上.

2.3 鄰域圓法

所謂鄰域圓法就是在一個(gè)觀測(cè)點(diǎn)的周圍以該點(diǎn)為圓心畫一個(gè)圓，根據(jù)圓內(nèi)的數(shù)據(jù)計(jì)算均值和相對(duì)標(biāo)準(zhǔn)差，然后對(duì)每個(gè)觀測(cè)點(diǎn)都做同樣的計(jì)算，通過圓內(nèi)均值和相對(duì)標(biāo)準(zhǔn)差在上是否平穩(wěn).由于這種方法可以根據(jù)氣象站點(diǎn)之間的距離來計(jì)算出合適的鄰域圓半徑并且可以準(zhǔn)確地遍歷整個(gè)研究區(qū)域，比起移動(dòng)窗口法［11］，這種方法更適合我國稀疏的氣象站點(diǎn)數(shù)據(jù).

標(biāo)準(zhǔn)差描述的是區(qū)域化變量觀測(cè)值間的離散情況，因此計(jì)算其標(biāo)準(zhǔn)差就相當(dāng)于求得離散情況，但是由于均值不同時(shí)變異程度就不同，即區(qū)域內(nèi)的標(biāo)準(zhǔn)差依賴于區(qū)域內(nèi)的均值，均值不同不能直接采用標(biāo)準(zhǔn)差，因此這里采用相對(duì)標(biāo)準(zhǔn)差，也稱變異系數(shù)，是一個(gè)無量綱的統(tǒng)計(jì)量，能反映單位均值上的離散程度，可以消除單位或平均數(shù)不同對(duì)兩個(gè)或多個(gè)資料變異程度比較的影響［14］.用公式表示為：

式中s表示圓內(nèi)樣點(diǎn)標(biāo)準(zhǔn)差，m表示圓內(nèi)樣點(diǎn)均值.

如果數(shù)據(jù)平穩(wěn)，圓內(nèi)的均值應(yīng)該圍繞某個(gè)固定常數(shù)上下隨機(jī)波動(dòng)，沒有明顯的趨勢(shì)，相對(duì)標(biāo)準(zhǔn)差應(yīng)該基本在0值附近波動(dòng).

3 結(jié)果分析

本文用R語言實(shí)現(xiàn)了對(duì)我國東北地區(qū)氣象要素（氣溫，風(fēng)速，降水）的空間平穩(wěn)性檢驗(yàn).R語言有著優(yōu)秀的統(tǒng)計(jì)分析能力和繪圖環(huán)境.R里面自帶的包中包含了統(tǒng)計(jì)學(xué)家們開發(fā)的前沿算法.通過直接調(diào)用這些包里面的函數(shù)（如persp3d（））可以避免通過大量編程實(shí)現(xiàn)統(tǒng)計(jì)算法.

3.1 三維散點(diǎn)圖法

利用rgl包中的persp3d（）函數(shù)繪制三維散點(diǎn)圖［15］，如圖2所示.

圖2 中國東北地區(qū)氣象觀測(cè)數(shù)據(jù)散點(diǎn)圖（氣溫，風(fēng)速，降水）Fig.2 Scatterplots of meteorological data in the Northeast area of China（temperature，wind speed and precipitation）

從圖2中可以看出，日平均氣溫有一個(gè)由北向南線性增加的趨勢(shì)，且增幅明顯.東北地區(qū)2006年1月1日的日平均氣溫從高緯地區(qū)的－25℃到低緯地區(qū)的－5℃.但是氣溫變量的隨機(jī)波動(dòng)很小，都是圍繞趨勢(shì)面小幅度的上下波動(dòng).風(fēng)速有一個(gè)自西向東減少的趨勢(shì)，趨勢(shì)不明顯.但是風(fēng)速變量的隨機(jī)波動(dòng)很大，即數(shù)據(jù)偏離趨勢(shì)面的幅度很大.降水?dāng)?shù)據(jù)的趨勢(shì)面是個(gè)明顯的曲面，趨勢(shì)面不平直，變量的隨機(jī)性一般.

3.2 趨勢(shì)線法

利用r中的l m（）和predict（）函數(shù)繪制趨勢(shì)線，如圖3所示.

從圖3中可以看出，氣溫?cái)?shù)據(jù)有一個(gè)極為明顯的由北向南增加的趨勢(shì)，由西向東方向則微弱減少.風(fēng)速數(shù)據(jù)基本上可以看成是一個(gè)平穩(wěn)的數(shù)據(jù)集，不論是在南北還是東西方向上都圍繞平均風(fēng)速2 m／s上下隨機(jī)波動(dòng).而5月的降水?dāng)?shù)據(jù)由西向東增加，由南向北減少，有著明顯的趨勢(shì)且隨機(jī)波動(dòng)較大.

3.3 鄰域圓檢驗(yàn)

圖3 東北地區(qū)氣象觀測(cè)數(shù)據(jù)趨勢(shì)線圖（氣溫，風(fēng)速，降水）Fig.3 Trend line maps of meteorological data in the Northeast area of China（temperature，wind speed and precipitation）

用R實(shí)現(xiàn)鄰域圓法的思路為：定義圓半徑值的大小，計(jì)算第一個(gè)點(diǎn)與所有點(diǎn)（包括第一個(gè)點(diǎn)本身）之間的距離，存儲(chǔ)距離小于半徑值的點(diǎn)并計(jì)算其均值和相對(duì)標(biāo)準(zhǔn)差，得到第一個(gè)點(diǎn)鄰域圓內(nèi)的均值和相對(duì)標(biāo)準(zhǔn)差.同理可計(jì)算得所有點(diǎn)鄰域圓內(nèi)的均值和相對(duì)標(biāo)準(zhǔn)差，并作圖顯示以判斷該區(qū)域化變量在該空間上是否平穩(wěn).鄰域圓按照從西向東、從北向南的順序遍歷整個(gè)研究區(qū)域.其結(jié)果如圖4～圖6所示.

圖4 東北地區(qū)日平均氣溫?cái)?shù)據(jù)鄰域圓法檢驗(yàn)圖Fig.4 Neighborhood circle plots of daily average temperature data in the Northeast area of China

圖5 東北地區(qū)日平均風(fēng)速數(shù)據(jù)鄰域圓法檢驗(yàn)圖Fig.5 Neighborhood circle plots of daily average wind speed data in the Northeast area of China

圖6 東北地區(qū)月降水?dāng)?shù)據(jù)鄰域圓法檢驗(yàn)圖Fig.6 Neighborhood circle plots of monthly precipitation data in the Northeast area of China

從圖4～圖6中可以看到2006年1月1日的風(fēng)速數(shù)據(jù)可以看成是平穩(wěn)的數(shù)據(jù)，而氣溫?cái)?shù)據(jù)和降水?dāng)?shù)據(jù)均不平穩(wěn).鄰域圓內(nèi)的氣溫均值逐步上升，由于氣溫圍繞均值波動(dòng)的幅度始終較小，因此相對(duì)標(biāo)準(zhǔn)差隨著均值的增加而逐步減小.鄰域圓內(nèi)的降水均值呈螺旋式上升的趨勢(shì)，相對(duì)標(biāo)準(zhǔn)差從樣本點(diǎn)號(hào)1到7是平穩(wěn)的，從9號(hào)樣本點(diǎn)到81也是平穩(wěn)的.但是從樣本號(hào)7到9之間降水?dāng)?shù)據(jù)的相對(duì)標(biāo)準(zhǔn)差急劇下降.樣本點(diǎn)號(hào)為7，8，9的氣象站如圖所示，由一個(gè)三角形圍住.從圖中可以看出，由于1～7號(hào)樣本點(diǎn)處于東北地區(qū)的西北角上，分布稀疏且呈一個(gè)狹長的形狀，因此落在每個(gè)鄰域圓的樣本數(shù)相對(duì)較少，標(biāo)準(zhǔn)差較大.而1～7號(hào)樣本點(diǎn)周圍鄰域圓內(nèi)的樣本均值卻比較小.所以1～7號(hào)樣本點(diǎn)周圍鄰域圓內(nèi)的相對(duì)標(biāo)準(zhǔn)差較大.

4 結(jié)論與展望

本文介紹的3種平穩(wěn)性檢驗(yàn)方法各有優(yōu)缺點(diǎn)，三維散點(diǎn)圖法簡(jiǎn)單直觀，可以迅速看出數(shù)據(jù)的均值是否是一個(gè)水平的平面.如果是一個(gè)平直的面那么數(shù)據(jù)均值平穩(wěn)，反之不平穩(wěn).三維散點(diǎn)的缺點(diǎn)在于不便于發(fā)現(xiàn)均值不平穩(wěn)的方向性，以及二階矩是否平穩(wěn).趨勢(shì)線法把空間數(shù)據(jù)投影到東西和南北兩個(gè)方向上，通過兩個(gè)方向上的剖面圖可以清晰地看出空間不平穩(wěn)的方向性.鄰域圓法同時(shí)顯示了局部區(qū)域內(nèi)樣本點(diǎn)的均值圖和相對(duì)標(biāo)準(zhǔn)差圖，它可以同時(shí)清楚地展現(xiàn)一階矩和二階矩是否平穩(wěn).3種方法同時(shí)使用，對(duì)同樣3組數(shù)據(jù)（我國東北地區(qū)日氣溫，日風(fēng)速，月降水）進(jìn)行空間平穩(wěn)性分析，得到相同的結(jié)論.即日氣溫?cái)?shù)據(jù)有一個(gè)由北向南迅速遞增的線性趨勢(shì)，這與我國大部分地區(qū)處于北半球的溫帶，氣溫南高北低的氣候規(guī)律相吻合.東北地區(qū)日風(fēng)速數(shù)據(jù)在空間上比較平穩(wěn)，始終在2 m／s的平均值上下波動(dòng)，但是隨機(jī)波動(dòng)幅度較大，從0.5 m／s到4 m／s不等.月降水量有著東高西低，南高北低的趨勢(shì).我國東北地區(qū)南面臨近渤海、黃海，東面臨近日本海，臨海地區(qū)空氣濕潤，雨量相對(duì)充沛.降水的空間趨勢(shì)檢驗(yàn)結(jié)果也與我國東北地區(qū)的地理環(huán)境相吻合.

本文有效檢驗(yàn)了中國東北地區(qū)氣象要素的空間平穩(wěn)性.在統(tǒng)計(jì)建模中如何處理不平穩(wěn)的空間數(shù)據(jù)是一個(gè)有待進(jìn)一步研究的問題.

［1］ Berndtsson R.Temporal variability in spatial correlation of daily rainfal［J］.Water Resources Research，1988（24）：1511－1517.

［2］ Obled C，Creutin J D.Some developments in the use of empirical orthogonal functions for mapping meteorological Field［J］.Journal of Climate and Applied Meteorolog，1986（25）：1189－1204.

［3］ Loader C，Switzer P.Spatial covariance estimation for monitoring data［R］.Stanford：Department of statistics，Stanford University，1989.

［4］ Brunsdon C，F(xiàn)otheringham S，Charlton M.Geographically weighted regression－modeling spatial non－stationarity［J］.The Statistician，1998（47）：431－443.

［5］ Fuentes M.Spectral method for nonstationary spatial processes［J］.Biometrika，2002，1（89）：197－210.

［6］ Fuentes M.A formal test for nonstationarity of spatial stochastic processes［J］.Journal of Multivariate Analysis，2005（96）：30－54.

［7］ Goovaerts P.Geostatistics for natural resources evaluation［M］.New York，Oxford：Oxford University Press，1997.

［8］ Sampson P D，Guttorp P.Nonparametric estimation of nonstationary spatial covariance structure［J］.Journal of the A－merican Statistical Association，1992，87（417）：108－119.

［9］ 湯國安，楊 昕.ArcGIS地理信息系統(tǒng)空間分析實(shí)驗(yàn)教程［M］.北京：科學(xué)出版社，2006.

［10］ Journel A G，Huijbregts C J.Mining Geostatistics［M］.London，New York，San Francisco：Academic Press，1978：30－39.

［11］ Hass T C.Lognormal and moving window methods of estimating acid deposition［J］.Journal of the American Statistical Association，1990，85（412）：950－963.

［12］ 周 琳.東北氣候［M］.北京：氣象出版社，1991：25－28.

［13］ 嚴(yán) 冰，陽正熙，吳塹虹，等.地學(xué)數(shù)據(jù)分析教程［M］.北京：科學(xué)出版社，2008：91－96.

［14］ 常文淵，戴新剛，陳洪武.地質(zhì)統(tǒng)計(jì)學(xué)在氣象要素場(chǎng)插值的實(shí)例研究［J］.地球物理學(xué)報(bào).2004，6（47）：982－989.

［15］ 湯銀才.R語言與統(tǒng)計(jì)分析［M］.北京：高等教育出版社，2008.