沈 艷, 潘 ,徐 賓, 宇婧婧

(國家氣象信息中心,北京100081)

0 引言

中國已建成自動氣象觀測站(包括國家級自動站和區(qū)域自動站)約32000多個,這些資料對于氣象災(zāi)害預(yù)警、決策服務(wù)、數(shù)值預(yù)報模式評估和精度檢驗等具有重要意義。但受目前資料交換政策的限制,用戶無法獲得自動站觀測的小時數(shù)據(jù),因此如何將該數(shù)據(jù)快速、準(zhǔn)確地轉(zhuǎn)換成可供決策服務(wù)部門和廣大科研用戶使用的高質(zhì)量產(chǎn)品,是綜合提高氣象防災(zāi)減災(zāi)能力、服務(wù)能力和科研水平的有效途徑。

利用特定的空間分析方法,已有一些研究提出了適宜于日、侯或月時間尺度的降水量空間插值方法[1-10]。chen等[9-11]比較了Cressman、Shepard和最優(yōu)插值(Optimal Interpolation,OI)3種方法插值逐日降水量時的精度表明:與只考慮距離權(quán)重的Cressman或Shepard插值不同,OI方法應(yīng)用到空間內(nèi)插時,除了距離還要考慮格點與站點以及站點與站點之間的空間相關(guān)性,其插值結(jié)果優(yōu)于Cressman和Shepard兩種方法。近年來,OI方法也應(yīng)用到了小時降水量的空間內(nèi)插中,如李春暉等[12]利用自動觀測站數(shù)據(jù),在中國站點分布相對均勻的廣東省和廣西省,進(jìn)行了區(qū)域小時降水量的插值試驗,認(rèn)為降水比率OI方法為最優(yōu)。作者同時指出,為了進(jìn)一步提高小時降水量的空間插值精度,需要對核心參數(shù)進(jìn)行優(yōu)化。實際上中國自動站的分布非常不均衡,“東密西疏”的站網(wǎng)布局非常明顯,約90%的自動站都建在了中國東部區(qū)域。在站點稠密區(qū)與某臺站距離最近的臺站間平均距離不到20km,但在站點稀疏區(qū)能達(dá)到100km以上。因此,針對現(xiàn)有的自動站站網(wǎng)布局,合理優(yōu)化插值參數(shù),對于提高中國小時降水量空間插值場精度是必要的。將中國分成3個區(qū)域,確定了OI插值中數(shù)據(jù)信噪比、相關(guān)隨距離的函數(shù)關(guān)系。討論影響臺站數(shù)和初估場對小時降水量精度地影響。

1 數(shù)據(jù)和方法

1.1 中國自動站小時觀測降水量

截止到2010年10月,中國已建成32742個自動觀測站,站點的空間分布圖見圖1?？梢?自動站在中國的分布極不均勻,約80%的站點分布在中國東南沿海地區(qū),而青藏高原等地形復(fù)雜區(qū)的站點非常稀疏。為此,以40°N和100°E將中國分成3個區(qū)域(圖1),分別對應(yīng)中國站點較為稀疏的新疆北部、內(nèi)蒙古大部和東三省(簡稱1區(qū))、站點較為稠密的中國東部區(qū)域(簡稱2區(qū))和青藏高原區(qū)域(簡稱3區(qū))。自動站觀測資料已實現(xiàn)了實時上傳,針對降水量的質(zhì)量控制方案也已于2009年7月實現(xiàn)了業(yè)務(wù)化[13],方案充分考慮了氣候?qū)W界限值、區(qū)域界限值、時間一致性、空間一致性等檢查,顯著提高了降水資料的可靠性。選取了2010年7月經(jīng)質(zhì)量控制后的自動站小時觀測降水量數(shù)據(jù)并將其格點化到0.1°×0.1°網(wǎng)格點上。

圖1 中國自動站的分布及其區(qū)域劃分

1.2 方法

Cressman、站點算術(shù)平均法、Shepard、Kriging等空間插值方法均是氣象上常用的空間分析方法。但Xie et al.[14]、沈艷等[15]、李春暉等表明:考慮降水比率的OI方法在降水量空間插值時具有優(yōu)勢。

1.2.1 降水比率OI插值方法

降水空間變率大導(dǎo)致直接對其插值的誤差也相對較大。因此,通常的插值思路是:先構(gòu)建一個初估場,進(jìn)一步借助初估場來定義一個新的要素,如降水量的差值[9]或比值[14]進(jìn)行空間插值,以減小因降水量空間分布不連續(xù)而帶來的插值誤差。Xie等的研究表明:定義比值能獲得更高的插值精度[14]。因而,針對小時降水的空間插值中也采用了通過初估場定義小時降水量比值的插值思路[12,15]。

1.2.2 精度評估方法

站點觀測數(shù)據(jù)是最能“真實”反映觀測點的數(shù)據(jù)源,因此在評估OI核心參數(shù)精度時,采用交叉檢驗方法,將站點數(shù)隨機(jī)分成10等份,利用90%的站點值反插到另外10%的臺站位置上,如此循環(huán)10次,使每一個臺站都有一個實測值和一個內(nèi)插值,然后利用平均偏差、均方根誤差和相關(guān)系數(shù)等指標(biāo)來定量評估調(diào)整OI方法中核心參數(shù)時的效果。

圖2 不同分區(qū)在ee取不同值時對應(yīng)的均方根誤差

1.3 OI方法參數(shù)設(shè)置

1.3.1 觀測誤差標(biāo)準(zhǔn)差和初估場誤差標(biāo)準(zhǔn)差的比率ee

這是OI方法中非常重要的參數(shù)之一,表示觀測數(shù)據(jù)誤差與初估場誤差的比率,可近似認(rèn)為是數(shù)據(jù)的信噪比。其中觀測數(shù)據(jù)誤差一是來源于觀測誤差,但通常認(rèn)為站點觀測數(shù)據(jù)的誤差可近似忽略,另外一個誤差來源是空間代表性誤差。而初估場誤差是時空分辨率和降水量的函數(shù),一般很難被量化[14]。因此在3個分區(qū)中,通過人為變換不同的ee值(從0.01變換到1.5),并將誤差最小時對應(yīng)的ee值作為最優(yōu)的觀測誤差標(biāo)準(zhǔn)差和初估場誤差標(biāo)準(zhǔn)差的比率。3個分區(qū)中,ee取不同值時對應(yīng)的均方根誤差見圖2,可見:不同區(qū)域下,ee取不同值時均方根誤差的變化趨勢比較一致,即隨著ee增大,各區(qū)域的均方根誤差均有小幅減小,之后隨著ee的增大,均方根誤差基本保持不變。各區(qū)域?qū)?yīng)的均方根誤差均在ee=0.3時達(dá)到最小,因此,在針對自動站小時觀測降水量空間插值時,比率確定為0.3。

1.3.2 相關(guān)隨距離的衰減關(guān)系

利用OI方法進(jìn)行數(shù)據(jù)融合或同化時,需要確定兩點之間誤差相關(guān)與距離的定量函數(shù)關(guān)系式。站點資料空間分析時,一般認(rèn)為站點觀測不存在誤差,或者認(rèn)為這種誤差可以忽略不計,因此可以簡化成是兩點間降水量時間相關(guān)系數(shù)與兩點間距離的關(guān)系式。選取2010年7月份逐時降水觀測資料(744個時次),分別計算了3個分區(qū)任意兩點間降水量相關(guān)系數(shù)與距離。進(jìn)一步按照距離大小分成多個等份,計算每個等份內(nèi)所有樣本的平均相關(guān)系數(shù),圖3顯示了3個分區(qū)降水量相關(guān)系數(shù)與距離之間的關(guān)系及函數(shù)擬和式?？梢?個分區(qū)中兩點間降水量相關(guān)系數(shù)均隨距離增加而減小,均可近似表達(dá)成隨距離增加而呈e指數(shù)遞減的關(guān)系,且3個區(qū)的擬和關(guān)系式基本一致。因此,在中國范圍可近似用下式來表達(dá)相關(guān)隨距離的衰減關(guān)系,即:

Cor=0.3*exp(-d/10.0)

其中Cor為兩點間降水量的時間相關(guān)系數(shù),d表示兩點間的距離。

圖3 不同分區(qū)兩點間降水量相關(guān)系數(shù)與距離之間的關(guān)系

1.3.3 影響臺站數(shù)和平均搜索距離

與只用距離作為權(quán)重的空間插值方法不同,OI方法在考慮搜索半徑的同時還要考慮搜索半徑內(nèi)臺站個數(shù)的影響,而且影響臺站數(shù)的判斷條件優(yōu)先于搜索半徑。選取2010年7月1日00時～2010年7月10日23時的逐時觀測資料,當(dāng)影響臺站數(shù)從1變化到15時,逐時次交叉檢驗計算均方根誤差(rmse)和平均搜索距離在3個區(qū)域的變化如圖4,可見3個區(qū)域的插值誤差均隨著影響臺站個數(shù)的增加而減小,之后維持著相對穩(wěn)定的誤差水平。對于1區(qū),當(dāng)影響臺站個數(shù)達(dá)到4個及以上時,rmse就相對穩(wěn)定,而對于2區(qū)和3區(qū),該數(shù)值分別是5個和7個。因此,在針對逐小時降水量的OI插值方案中,影響臺站數(shù)應(yīng)至少在4個以上,以5～7個之間為宜。

圖4 研究區(qū)均方根誤差和平均搜索距離隨影響臺站個數(shù)的交叉檢驗結(jié)果

1.3.4 初估場

OI方法中需要提前建立初估場,在針對逐日、逐月站點數(shù)據(jù)的空間分析中,將日或月的氣候標(biāo)準(zhǔn)值數(shù)據(jù)作為了初估場[11,14-15]。但由于小時降水的時空變率大,建立穩(wěn)定可靠的小時降水量背景場需要更多的氣象臺站觀測資料,而現(xiàn)有的資料質(zhì)量和密度尚無法建立該初估場。因此在針對小時降水量的空間分析中仍采用了逐日降水量氣候標(biāo)準(zhǔn)值作為初估場[12]。設(shè)計了3種初估場方案,即(1)NOAA氣候預(yù)測中心(CPC)研制的考慮地形海拔高度影響的逐日降水量背景場[14];(2)國家氣象信息中心制作的綜合考慮了海拔、坡度、坡向與風(fēng)向夾角的逐日降水量背景場[16];(3)用Shepard空間插值法[3]插值生成同時刻的降水分析場。采用3種初估場的交叉檢驗結(jié)果見表1。由于生成初估場1和初估場2使用的均是國家級2400多臺站的日降水量標(biāo)準(zhǔn)值,盡管采用的插值算法略有不同,但對2套數(shù)據(jù)造成的插值結(jié)果對比發(fā)現(xiàn):利用初估場1和2得到的不同區(qū)域插值誤差的交叉檢驗結(jié)果基本一致,區(qū)域 2的相對偏差(R Bias)在10.0%以內(nèi),區(qū)域 1和 3的相對偏差分別超過-20.0%和-30.0%。由于初估場2進(jìn)一步考慮了坡度、坡向與風(fēng)向夾角對降水量的影響,因而能更準(zhǔn)確得再現(xiàn)降水的分布特征,交叉檢驗的系統(tǒng)偏差均小于初估場1對應(yīng)的結(jié)果。采用初估場3雖然顯著降低了區(qū)域1和區(qū)域3的相對偏差,如區(qū)域1的相對偏差由-21.1%降低到了-5.6%,區(qū)域3的相對偏差由-34.2%降低到了-9.1%,但各區(qū)域的均方根誤差(Rmse)均有所增大,其中區(qū)域2的Rmse由0.956增加到了1.250mm/hr。Rmse增大的原因可以從插值數(shù)據(jù)的時間尺度上得到一定解釋:初估場3是由小時數(shù)據(jù)內(nèi)插得到的,與用30年氣候標(biāo)準(zhǔn)值內(nèi)插生成的初估場1和2相比,必然包含更大的隨機(jī)性,其結(jié)果將導(dǎo)致隨機(jī)誤差即Rmse變大。

表1 采用3種初估場得到的區(qū)域插值誤差的交叉檢驗結(jié)果

2 新舊參數(shù)下的交叉檢驗結(jié)果

利用優(yōu)化后的OI核心參數(shù),進(jìn)一步分析了新舊參數(shù)下的交叉檢驗結(jié)果,數(shù)據(jù)統(tǒng)計時間段是2010年7月,新舊參數(shù)在不同區(qū)域的交叉檢驗結(jié)果見表2?？梢姴煌瑓^(qū)域采用新參數(shù)后的均方根誤差均小于舊參數(shù)下的相應(yīng)值,如區(qū)域2的Rmse由舊參數(shù)下的1.409mm/hr減小到了新參數(shù)下的1.328mm/hr。說明采用新參數(shù)獲得的降水量值與站點觀測更為接近。新參數(shù)在不同區(qū)域?qū)?yīng)的相關(guān)系數(shù)也得以提高,如全國平均的相關(guān)系數(shù)從舊參數(shù)下的0.257提高到了新參數(shù)下的0.294,說明采用新參數(shù)得到的降水空間分布型與站點更為接近。

表2 新舊參數(shù)下得到的不同區(qū)域交叉檢驗統(tǒng)計結(jié)果

進(jìn)一步統(tǒng)計了新舊參數(shù)交叉檢驗統(tǒng)計指標(biāo)的頻率分布,數(shù)據(jù)統(tǒng)計時間段是2010年7月,圖5僅給出全國平均結(jié)果,不同區(qū)域的統(tǒng)計結(jié)果與全國類似(圖略)。從系統(tǒng)偏差(Bias)的頻率分布而言:新參數(shù)對站點數(shù)據(jù)的低估較舊參數(shù)比例大,而對站點高估的比例小于舊參數(shù)。尤其是系統(tǒng)偏差在[-0.005,0]時,舊參數(shù)對應(yīng)的樣本百分率為40.6%,而新參數(shù)只有28.8%。從Rmse的頻率分布看:新、舊參數(shù)在Rmse＜1.5mm/hr所占的樣本百分率合計為63.3%和56.6%,而舊參數(shù)在Rmse≥1.5mm/hr時對應(yīng)的樣本百分率均大于新參數(shù)所占的比率,說明新參數(shù)能獲得更小的均方根誤差。新參數(shù)得到的相關(guān)系數(shù)也較舊參數(shù)有所提高,如在相關(guān)系數(shù)比較低的[0,0.2]區(qū)間,舊參數(shù)出現(xiàn)頻率為28.3%,而新參數(shù)只占17.3%。隨著相關(guān)系數(shù)的提高,新參數(shù)出現(xiàn)高相關(guān)系數(shù)的頻率均比舊參數(shù)有所提高,如新參數(shù)在相關(guān)系數(shù)≥0.4的樣本百分率為14.1%,而對應(yīng)的舊參數(shù)下只占6.7%。

圖5 新舊參數(shù)交叉檢驗統(tǒng)計指標(biāo)(Bias,Rmse,Cor Coe)的頻率分布

3 結(jié)論與討論

OI插值的核心參數(shù)包括:數(shù)據(jù)信噪比、相關(guān)隨距離的變化關(guān)系、搜索半徑(或表達(dá)成影響臺站數(shù))以及初估場等。進(jìn)一步考慮自動站的空間分布將中國分成了三大區(qū)域,探討了利用OI對全國3萬多個自動站小時降水?dāng)?shù)據(jù)進(jìn)行空間分析時核心參數(shù)的優(yōu)化和確定方法。研究表明:“優(yōu)化的”數(shù)據(jù)信噪比約為0.3,相關(guān)隨距離增加呈e指數(shù)遞減關(guān)系;當(dāng)站網(wǎng)密度差異較大時,通過確定影響臺站數(shù)來進(jìn)行空間插值,一方面可以提高站點稠密區(qū)的插值效率,另一方面可以提高站點相對稀疏區(qū)的插值精度。另外,初估場對插值結(jié)果也有影響,認(rèn)為采用氣候標(biāo)準(zhǔn)值生成的初估場能顯著降低插值結(jié)果的隨機(jī)誤差,結(jié)果比采用同時刻的降水分析場作為初估場更穩(wěn)定。進(jìn)一步采用站點交叉檢驗的方法對新舊參數(shù)下的插值結(jié)果(以2010年7月數(shù)據(jù))進(jìn)行了對比分析:新參數(shù)對應(yīng)的均方根誤差均比舊方案小,而空間相關(guān)系數(shù)有一定提高,說明采用新參數(shù)能獲得更加可靠真實的1小時、0.1°分辨率的降水量分析場。

為了進(jìn)一步提高中國站點稀疏區(qū)的小時降水資料質(zhì)量,已將衛(wèi)星反演的降水產(chǎn)品(CMORPH和FY)融合到了自動站降水量空間分析場中,獲得了質(zhì)量和時空分辨率均更高的1小時、0.1°分辨率的降水量融合資料,產(chǎn)品已在中國氣象科學(xué)數(shù)據(jù)共享網(wǎng)實時發(fā)布(cdc.cma.gov.cn)。

致謝:文章中使用的OI插值程序(FORTRAN語言)由NOAA氣候預(yù)測中心Xie Pingping博士提供,在此深表感謝。

[1] Cressman G P.An operational objective analysis system[J].Mon.Wea.Rev.,1959,87:367-374.

[2] Gandin L S.Objective Analysis of meteorological fields[J].Israel Program for Scientific Translations,1965,242.

[3] Shepard D.A two dimensional interpolation function for irregularly spaced data[C].Prod.23rd National Conf.of the Association for Computing Machinery,Princeton,NJ,ACM,1968:517-524.

[4] New,M,M Hulme,P Jones.Representing twentieth century space-time climate variability.Part I:Development of a 1961-90 mean monthly terrestrial climatology[J].J.Climate,1999,12:829-856.

[5] New M,G M Hulme,P D Jones.Representing Twentieth-Century Space Time Climate Variability.Part II:Development of 1901-96 Monthly Grids of Terrestrial Surface Climate[J].J.Climate,2000,13:2217-2238.

[6] Rudolf B.Management and analysis of precipitation data on a routine basis[C].Proc.Int.Symp.on Precipitation and Evaporation,Bratislava,Slovakia,WMO,1993:69-76.

[7] Schneider U.The GPCC quality-control system for gauge measured precipitation data[C].Proc.Analysis Methods for Precipitation on a Global Scale:Report of a GEWEX Workshop,WCRP-81,WMO/TD-588,Koblenz,Germany,1993.

[8] Hutchinson M F.Interpolating mean rainfall using thin plate smoothing splines[J].Int.J.Geogr.Inf.Syst.,1995,9:385-403.

[9] Chen M,Xie P,Janowiak J E.Global land precipitation:A 50-yr monthly analysis based on gauge observations[J].J.Hydrometeorol,2002,(3):249-266.

[10] Rudolf B,U Schneider.Calculation of gridded precipitation data for the global land-surface using in-situ gauge observations[C].Proc.the 2nd Workshop of the International Precipitation Working Group IPWG,Monterey,2004,231-247.

[11] Chen M,W Shi,P Xie.2008:Assessing objective techniques for gauge-based analyses of global daily precipitation[J].J.Geophys.Res.,2008,113.

[12] 李春暉,梁建茵.基于Shepard和OI方法對雨量計逐時資料的分析[J].應(yīng)用氣象學(xué)報,2010,21(4):416-422.

[13] 任芝花,趙平,張強(qiáng),等.適用于全國自動站小時降水資料的質(zhì)量控制方法[J].氣象,2010,36(7):123-132.

[14] Xie P P,Yatagai A,Chen M,et al.A gauge-based analysis of daily precipitation over East Asia[J].J.Hydrometeor,2007,8:607-626.

[15] 沈艷,馮明農(nóng),張洪政,等.我國逐日降水量格點化方法[J].應(yīng)用氣象學(xué)報,2010,21(3):279-286.

[16] 徐賓.中國地面降水逐日氣候背景0.25°格點場數(shù)據(jù)集說明文檔[J].2009.

[18] Arkin P A,Ardanuy P E.Estimating climatic-scale precipitation from space:A review[J].J Climate,1989,(2):1229-1238.

[19] Barrett E C,Martin D W.The user of satellite data in rainfall monitoring[M].Academic Press,1981.

[20] Franke R.Smooth interpolation of scattered data by local thin plate splines[J].Comput.Math.Appl.,1982,8:273-281.

[21] Huffman G J.Estimates of Root-Mean-Square Random Error for Finite Samples of Estimated Precipitation[J].Journal of Applied Meteorology,1997,36(9):1191-1201.