最優插值法在對中國自動站降水量空間分析中的參數優化

沈 艷, 潘 ,徐 賓, 宇婧婧

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

0 引言

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

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

1 數據和方法

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

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

圖1 中國自動站的分布及其區域劃分

1.2 方法

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

1.2.1 降水比率OI插值方法

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

1.2.2 精度評估方法

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

圖2 不同分區在ee取不同值時對應的均方根誤差

1.3 OI方法參數設置

1.3.1 觀測誤差標準差和初估場誤差標準差的比率ee

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

1.3.2 相關隨距離的衰減關系

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

Cor=0.3*exp(-d/10.0)

其中Cor為兩點間降水量的時間相關系數,d表示兩點間的距離。

圖3 不同分區兩點間降水量相關系數與距離之間的關系

1.3.3 影響臺站數和平均搜索距離

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

圖4 研究區均方根誤差和平均搜索距離隨影響臺站個數的交叉檢驗結果

1.3.4 初估場

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

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

2 新舊參數下的交叉檢驗結果

利用優化后的OI核心參數,進一步分析了新舊參數下的交叉檢驗結果,數據統計時間段是2010年7月,新舊參數在不同區域的交叉檢驗結果見表2。可見不同區域采用新參數后的均方根誤差均小于舊參數下的相應值,如區域2的Rmse由舊參數下的1.409mm/hr減小到了新參數下的1.328mm/hr。說明采用新參數獲得的降水量值與站點觀測更為接近。新參數在不同區域對應的相關系數也得以提高,如全國平均的相關系數從舊參數下的0.257提高到了新參數下的0.294,說明采用新參數得到的降水空間分布型與站點更為接近。

表2 新舊參數下得到的不同區域交叉檢驗統計結果

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

圖5 新舊參數交叉檢驗統計指標(Bias,Rmse,Cor Coe)的頻率分布

3 結論與討論

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

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

致謝:文章中使用的OI插值程序(FORTRAN語言)由NOAA氣候預測中心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].應用氣象學報,2010,21(4):416-422.

[13] 任芝花,趙平,張強,等.適用于全國自動站小時降水資料的質量控制方法[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] 沈艷,馮明農,張洪政,等.我國逐日降水量格點化方法[J].應用氣象學報,2010,21(3):279-286.

[16] 徐賓.中國地面降水逐日氣候背景0.25°格點場數據集說明文檔[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.