最優(yōu)插值法在對(duì)中國(guó)自動(dòng)站降水量空間分析中的參數(shù)優(yōu)化

沈 艷, 潘 ,徐 賓, 宇婧婧

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

0 引言

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

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

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

1.1 中國(guó)自動(dòng)站小時(shí)觀測(cè)降水量

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

圖1 中國(guó)自動(dòng)站的分布及其區(qū)域劃分

1.2 方法

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

1.2.1 降水比率OI插值方法

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

1.2.2 精度評(píng)估方法

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

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

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

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

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

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

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

Cor=0.3*exp(-d/10.0)

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

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

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

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

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

1.3.4 初估場(chǎng)

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

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

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

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

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

進(jìn)一步統(tǒng)計(jì)了新舊參數(shù)交叉檢驗(yàn)統(tǒng)計(jì)指標(biāo)的頻率分布,數(shù)據(jù)統(tǒng)計(jì)時(shí)間段是2010年7月,圖5僅給出全國(guó)平均結(jié)果,不同區(qū)域的統(tǒng)計(jì)結(jié)果與全國(guó)類似(圖略)。從系統(tǒng)偏差(Bias)的頻率分布而言:新參數(shù)對(duì)站點(diǎn)數(shù)據(jù)的低估較舊參數(shù)比例大,而對(duì)站點(diǎn)高估的比例小于舊參數(shù)。尤其是系統(tǒng)偏差在[-0.005,0]時(shí),舊參數(shù)對(duì)應(yīng)的樣本百分率為40.6%,而新參數(shù)只有28.8%。從Rmse的頻率分布看:新、舊參數(shù)在Rmse＜1.5mm/hr所占的樣本百分率合計(jì)為63.3%和56.6%,而舊參數(shù)在Rmse≥1.5mm/hr時(shí)對(duì)應(yīng)的樣本百分率均大于新參數(shù)所占的比率,說(shuō)明新參數(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%,而對(duì)應(yīng)的舊參數(shù)下只占6.7%。

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

3 結(jié)論與討論

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

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

致謝:文章中使用的OI插值程序(FORTRAN語(yǔ)言)由NOAA氣候預(yù)測(cè)中心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方法對(duì)雨量計(jì)逐時(shí)資料的分析[J].應(yīng)用氣象學(xué)報(bào),2010,21(4):416-422.

[13] 任芝花,趙平,張強(qiáng),等.適用于全國(guó)自動(dòng)站小時(shí)降水資料的質(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),張洪政,等.我國(guó)逐日降水量格點(diǎn)化方法[J].應(yīng)用氣象學(xué)報(bào),2010,21(3):279-286.

[16] 徐賓.中國(guó)地面降水逐日氣候背景0.25°格點(diǎn)場(chǎng)數(shù)據(jù)集說(shuō)明文檔[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.