郭家力，郭東淏，丁光旭，李 穎1，，李英海1，，楊 旭，張海榮

（1.水電工程施工與管理湖北省重點實驗室（三峽大學(xué)），湖北 宜昌 443002； 2.三峽大學(xué)水利與環(huán)境學(xué)院，湖北 宜昌 443002；3.智慧長江與水電科學(xué)湖北省重點實驗室（中國長江電力股份有限公司），湖北 宜昌 443133）

0 引 言

降水是最基本的氣象要素之一，準(zhǔn)確掌握流域降水的時空分布特征對實現(xiàn)數(shù)字孿生流域“四預(yù)”功能具有重要意義［1-3］。目前獲取降水資料的主要途徑主要有3種：地面雨量站、地基天氣雷達(dá)估測和氣象衛(wèi)星反演［4］，各自優(yōu)缺點十分明顯，但在一定程度上可以互補［5］。地面雨量站對降水量的觀測是獲得降水?dāng)?shù)據(jù)最直接也是最常用的方法，在整個降水觀測體系中起到基礎(chǔ)性的作用，可以獲得高精度的單點數(shù)據(jù)，但地面站點觀測無法覆蓋到無人區(qū)以及地形相對復(fù)雜的區(qū)域，從而限制了站點觀測數(shù)據(jù)的使用；氣象衛(wèi)星反演降水通過間接手段進行觀測，不受地理條件的限制，具有全天候、全覆蓋的特點，然而其精度會受到傳感器性能、反演算法以及云層性質(zhì)的影響，雖然可一定程度上反映降水的時空分布，但在單點上精度低于地面觀測數(shù)據(jù)［4，6］；雷達(dá)通過間接手段觀測降水，其優(yōu)點在于能夠跟蹤雨區(qū)范圍，在中小尺度地區(qū)的短期強降水監(jiān)測中具有較強的優(yōu)勢，但雷達(dá)對降水的觀測范圍有限，且易受地形等因素影響，估測的降水?dāng)?shù)據(jù)誤差較大［7］。受制于降水類型、地形地貌、反演算法等因素，單一降水產(chǎn)品難以準(zhǔn)確反映區(qū)域降水過程，導(dǎo)致各降水產(chǎn)品在不同區(qū)域有著不同的性能。

因此，如何通過多源降水?dāng)?shù)據(jù)融合的方式研制高質(zhì)量和高時空分辨率的降水產(chǎn)品，一直是水文氣象領(lǐng)域研究的熱點問題之一［8］。王筱譯［9］等將加權(quán)最小二乘估計法應(yīng)用于降水?dāng)?shù)據(jù)融合，引入SM2RAIN-CCI降水?dāng)?shù)據(jù)對IMERG數(shù)據(jù)進行校正，并以渭河流域2014年至2015年的地面實測降水?dāng)?shù)據(jù)進行效果檢驗，融合數(shù)據(jù)的各項統(tǒng)計誤差較原始IMERG數(shù)據(jù)均有較大的改善；譚偉偉［10］等提出了點面融合方法和站點偏差校正估計兩個多源數(shù)據(jù)融合方案，考慮了經(jīng)緯度、DEM、IMERG插值數(shù)據(jù)等輔助變量，融合了湖北省2016年7月19日的IMERG日尺度降水?dāng)?shù)據(jù)與氣象站點數(shù)據(jù)資料，有效提升了IMERG降水?dāng)?shù)據(jù)的精度；高真［11］等以地面雨量站降水?dāng)?shù)據(jù)作為參考，評估了NLDAS2、Stage IV、TMPARP、TMPART和IMERG降水?dāng)?shù)據(jù)的精度，并分別驅(qū)動分布式水文模型，發(fā)現(xiàn)多源遙感降水?dāng)?shù)據(jù)融合產(chǎn)品精度更高，更適合于水文模擬應(yīng)用；孟慶博［9］等基于地面觀測數(shù)據(jù)和TRMM、CMORPH、CHIRPS、PERSIAAN-CDR、GL‐DAS_Noah多源衛(wèi)星降水?dāng)?shù)據(jù)，使用了集合卡爾曼濾波算法進行數(shù)據(jù)融合同化，同化后的5種降水產(chǎn)品精度均有顯著提升，證明集合卡爾曼濾波算法在降水?dāng)?shù)據(jù)融合同化中有較大潛力。從上述研究可以看出，目前在多源降水融合研究方面，集合卡爾曼濾波被證明是一種具有較大潛力的方法，但目前的應(yīng)用研究卻開展極少。

基于此，本文將基于集合卡爾曼濾波算法融合多源降水資料，生成一套高精度的降水融合產(chǎn)品。首先，收集建立融合模型所需的多源數(shù)據(jù)，包括地面站點地形數(shù)據(jù)、實測數(shù)據(jù)和衛(wèi)星降水?dāng)?shù)據(jù)，然后對原始數(shù)據(jù)進行預(yù)處理，將所有網(wǎng)格降水?dāng)?shù)據(jù)重采樣至統(tǒng)一的空間分辨率，處理成融合模型所需要的數(shù)據(jù)，最后使用集合卡爾曼濾波算法對多源降水?dāng)?shù)據(jù)進行融合，得到一套0.05°×0.05°的融合降水?dāng)?shù)據(jù)。

1 研究區(qū)域概況與數(shù)據(jù)獲取

1.1 流域概況

清江屬于長江的一級支流，位于東經(jīng)108°58'～111°21'，北緯29°41'～30°87'之間，屬長江上游地區(qū)，流域位于湖北省境內(nèi)。清江干流發(fā)源于恩施州利川市境內(nèi)，全長423 km，平均比降3.2%，流域面積1.67 萬km2。清江流域內(nèi)部群山聳立、河流深切、高差懸殊，高山地區(qū)（海拔>1 200 m）占流域總面積的28%；中高山（海拔800～1 200 m）面積為1.2萬km2，占52%；低山丘陵區(qū)（<800 m）為4 652 km2，占20%，流域最大相對高差2 217 m。清江流域雨量豐沛，水汽主要來自太平洋季風(fēng)和印度洋季風(fēng)，年平均降水量為1 561 mm，最大年降雨量達(dá)1 700多 mm，最少年降水量僅900 mm左右，其中80%以上的年降水量集中在汛期5-10月，各月降水差異明顯，屬于亞熱帶季風(fēng)氣候。清江流域地理位置如圖1所示。

1.2 數(shù)據(jù)來源

1.1.1 地形高程數(shù)據(jù)

地形高程數(shù)據(jù)選擇SRTM DEM 90 m分辨率原始高程數(shù)據(jù)，數(shù)據(jù)從地理空間數(shù)據(jù)云下載（http：//www.gscloud.cn/）。

1.1.2 地面站點數(shù)據(jù)

本文收集到了流域內(nèi)部28個及流域周邊5個共33個地面觀測站點的逐日降水?dāng)?shù)據(jù)，站點分布如圖1所示，所選數(shù)據(jù)通過了嚴(yán)格的質(zhì)量控制，資料完整，數(shù)據(jù)可靠，最終得到1998年1月1日-1999年12月31日的逐日降水資料。在日尺度進行多源降水?dāng)?shù)據(jù)融合，數(shù)據(jù)序列長度730 d，基本滿足融合結(jié)果的精度檢驗要求。

1.1.3 衛(wèi)星降水?dāng)?shù)據(jù)

本文所用的衛(wèi)星降水?dāng)?shù)據(jù)包括TRMM（Tropical Rainfall Measuring Mission）和CMORPH（Climate Prediction Center morph‐ing technique）數(shù)據(jù)，其空間分辨率均為0.25°，時間分辨率為1 d。TRMM衛(wèi)星由美國國家航空航天局和日本宇宙航空研究開發(fā)機構(gòu)聯(lián)合研制，是專門用于定量測量熱帶、亞熱帶降雨的氣象衛(wèi)星，現(xiàn)已擴展到近似覆蓋全球范圍［13］，研究采用美國國家航空航天局發(fā)布的TMPA（TRMM Multi-satellite Precipitation Analysis）第7版日降水?dāng)?shù)據(jù)。CMORPH數(shù)據(jù)是美國國家海洋和大氣管理局發(fā)布的衛(wèi)星降水?dāng)?shù)據(jù)集［14］，本文采用的是CMORPH_V1.0_BLD日降水?dāng)?shù)據(jù)。

1.1.4 再分析降水?dāng)?shù)據(jù)

為了進一步驗證融合數(shù)據(jù)集的可靠性，選用CMFD（China Meteorological Forcing Dataset，CMFD）、ERA5（ECMWF Reanaly‐sis v5）和MSWEP（Multi-Source Weighted-Ensemble Precipita‐tion）三種再分析數(shù)據(jù)集與其進行對比。CMFD是中國科學(xué)院青藏高原研究所開發(fā)的近地面氣象與環(huán)境要素再分析數(shù)據(jù)，該數(shù)據(jù)以國際上已有的GLDAS（Global Land Data Assimilation Sys‐tem）等再分析數(shù)據(jù)和TRMM數(shù)據(jù)為背景數(shù)據(jù)，融合中國氣象局常規(guī)氣象觀測資料制作而成［15］。ERA5是歐洲天氣預(yù)報中心ECMWF（European Centre for Medium-Range Weather Forecasts）對全球氣候的第五代大氣再分析數(shù)據(jù)集，它使用了ECMWF綜合預(yù)報系統(tǒng)中的四維變分同化系統(tǒng)和模式預(yù)報系統(tǒng)，降水部分融合了TRMM等數(shù)據(jù)集、雷達(dá)降水?dāng)?shù)據(jù)集和地面氣象觀測資料，時間分辨率為6 h，空間分辨率為0.25°×0.25°［16］。MSWEP是美國普林斯頓大學(xué)研制的融合降水?dāng)?shù)據(jù)集，該產(chǎn)品基于Budyko框架和全球徑流觀測對地形復(fù)雜區(qū)域降水進行校正，較大程度上減小了估算誤差，并進一步融合了CMORPH等衛(wèi)星降水?dāng)?shù)據(jù)、ERA-Interim（ECMWF Reanalysis Interim）等再分析數(shù)據(jù)和地面氣象觀測資料［17］。以上多源降水?dāng)?shù)據(jù)總結(jié)于表1。

表1 地面站點、衛(wèi)星和再分析降水?dāng)?shù)據(jù)基本信息Tab.1 Basic information on ground stations，satellites and reanalysis of precipitation data

2 研究方法

2.1 空間插值方法

如表1所示，考慮到地面站點和衛(wèi)星降水?dāng)?shù)據(jù)的空間分辨率不一致，在數(shù)據(jù)融合之前需要在網(wǎng)格上統(tǒng)一。首先通過薄板光滑樣條模型（ANUSPLIN）插值方法將分布不規(guī)則的站點觀測數(shù)據(jù)插值為0.05°×0.05°分辨率的規(guī)則網(wǎng)格數(shù)據(jù)。ANUSPLIN是一種使用薄盤光滑樣條函數(shù)進行空間插值和分析的工具，在氣象變量插值中取得了優(yōu)異的效果，廣泛應(yīng)用于氣象數(shù)據(jù)集的構(gòu)建。除了獨立的樣條變量外，ANUSPLIN允許引入?yún)?shù)化的線性子模型或線性協(xié)變量子模型［18］。

清江流域地形復(fù)雜，海拔差異大，而地形對于降水的動力、熱力和微物理過程具有重要的影響，是導(dǎo)致區(qū)域降水量變化的一個主要因素［19］。因此選取經(jīng)度和緯度作為獨立變量，海拔作為協(xié)變量，對站點觀測數(shù)據(jù)進行插值，得到的結(jié)果將其規(guī)定為插值數(shù)據(jù)SI。

2.2 融合方法

集合卡爾曼濾波算法是集合預(yù)報與卡爾曼濾波方法的結(jié)合，它是貝葉斯信息估計理論的一種直接實現(xiàn)［20］。集合卡爾曼濾波算法通過蒙特卡羅方法計算狀態(tài)的預(yù)報誤差協(xié)方差，用集合的思想解決實際應(yīng)用中背景誤差協(xié)方差矩陣估計和預(yù)報困難的問題，有效降低了計算量。與標(biāo)準(zhǔn)的卡爾曼濾波相似，集合卡爾曼濾波算法主要分為預(yù)測和更新兩個步驟［21］。

在k=0時刻對狀態(tài)變量集合進行初始化（本文僅有降水這一個狀態(tài)變量），將k時刻的狀態(tài)變量代入模型，得到狀態(tài)量k+1時刻的預(yù)測值。

計算卡爾曼增益K，并對k+1時刻的觀測值對集合的狀態(tài)進行更新得到分析值。

將地面站點插值數(shù)據(jù)作為融合算法觀測值，將衛(wèi)星降水?dāng)?shù)據(jù)作為融合算法背景值，統(tǒng)一在0.05°×0.05°分辨率下的網(wǎng)格進行數(shù)據(jù)融合。

2.3 精度評定方法和指標(biāo)

為了評估融合模型的精度，并對比主流的再分析數(shù)據(jù)，采用留一交叉驗證法進行驗證。留一交叉驗證（Leave-One-Out Cross-Validation，LOO-CV）是評估模型泛化能力的方法，具體做法為：每次使用5個流域周邊氣象站點和清江流域內(nèi)28個氣象站點中的27個站點作為模型輸入，剩下的1個氣象站點留做精度評定，如此逐一循環(huán)28次，使得流域內(nèi)每個氣象站都參與了精度評定，得到插值數(shù)據(jù)在流域內(nèi)每個站點的交叉驗證值。衛(wèi)星數(shù)據(jù)、融合數(shù)據(jù)以及再分析數(shù)據(jù)則直接提取28個站點所在網(wǎng)格的數(shù)據(jù)進行精度評定。使用相關(guān)性系數(shù)R，平均絕對誤差MAE和均方根誤差RMSE作為具體的評價指標(biāo)，其表達(dá)式如公式（8）～（10）所示。

式中：Gi表示第i天的站點實測降水；Si表示第i天的衛(wèi)星降水（或插值降水、融合降水和再分析降水），mm；N表示序列長度，d。

3 結(jié)果與分析

3.1 融合模型精度評價

通過采用留一交叉驗證的方法，得出清江流域內(nèi)28個站點日尺度下各評價指標(biāo)。圖2為背景數(shù)據(jù)（TRMM、CMORPH和SI）、再分析數(shù)據(jù)（CMFD、ERA5和MSWEP）與MSAP融合數(shù)據(jù)評價指標(biāo)箱線圖。其中，紅色代表3種背景數(shù)據(jù)，綠色代表MSAP數(shù)據(jù)，藍(lán)色代表3種再分析數(shù)據(jù)。

圖2 日尺度融合數(shù)據(jù)相關(guān)性系數(shù)、平均絕對誤差、均方根誤差的箱線圖Fig.2 Boxplots of R， MAE， and RMSE of daily-scale fusion data

首先，對比背景數(shù)據(jù)和MSAP數(shù)據(jù)的精度。圖2（a）為相關(guān)性系數(shù)箱線圖：TRMM、CMORPH和SI的中位數(shù)分別為0.67、0.81和0.88，MSAP中位數(shù)為0.89，相較于TRMM、CMORPH和SI分別提高了32.1%、9.3%和1.1%；數(shù)值分布上，MSAP在各站點均大于0.8，而SI有3個站點小于0.8，可見MSAP略優(yōu)于SI。圖2（b）為平均絕對誤差箱線圖：TRMM中位數(shù)為3.7 mm/d，最小值為3.1 mm/d；CMORPH中位數(shù)為2.5 mm/d，數(shù)據(jù)分布于1.4～3.4 mm/d之間；SI中位數(shù)為1.9 mm/d，最大值為3.6 mm/d；MSAP中位數(shù)為1.9 mm/d，與SI持平，但優(yōu)于TRMM和CMORPH；數(shù)據(jù)分布上，MSAP最大值為2.9 mm/d，低于SI，數(shù)據(jù)分布均勻，可見MSAP在平均絕對誤差指標(biāo)上能基本保持與三種背景數(shù)據(jù)中的最優(yōu)者SI持平。圖2（c）為均方根誤差箱線圖：TRMM數(shù)據(jù)均在7 mm/d以上，誤差較大；CMORPH分布于4.0～10.9 mm/d之間，中位數(shù)為6.9 mm/d；SI分布于2.8～12.74 mm/d之間，中位數(shù)為5.5 mm/d，多個站點值大于8 mm/d；MSAP相較于TRMM和CMORPH分別提高了41.2%、24.3%，與SI持平；MSAP的均方根誤差數(shù)值分布上極差最小，集中分布于3.7～9.5 mm/d之間，最大值為9.5 mm/d，均小于3種背景數(shù)據(jù)的最大值。通過對3個評價指標(biāo)對比結(jié)果的定量分析可以得出結(jié)論，EnKF融合算法融合了各背景數(shù)據(jù)的優(yōu)勢，MSAP數(shù)據(jù)的精度優(yōu)于3種背景數(shù)據(jù)。

其次，對比MSAP與現(xiàn)有再分析降水?dāng)?shù)據(jù)的精度。從圖2綠色、藍(lán)色箱線圖對比可以看出，無論是相關(guān)性系數(shù)、平均絕對誤差還是均方根誤差，MSAP與3種再分析降水?dāng)?shù)據(jù)中的CMFD最為接近。MSAP相關(guān)性系數(shù)中位數(shù)為0.89，而CMFD為0.78，整體相關(guān)性系數(shù)優(yōu)于CMFD；盡管CMFD的平均絕對誤差的中位數(shù)小于MASP，但其分布范圍更寬，說明CMFD空間異質(zhì)性更大，而均方根誤差方面則是MASP明顯優(yōu)于CMFD。MSWEP的分布在3個指標(biāo)上略遜于CMFD，且無論是中位數(shù)還是分布范圍二者均類似。ERA5在3個指標(biāo)中表現(xiàn)均為最差，精度遠(yuǎn)低于CMFD、MSWEP和MSAP。總體上，MSAP在各個指標(biāo)上全面優(yōu)于3種再分析數(shù)據(jù)，精度排序為MSAP > CMFD > MSWEP >ERA5，進一步說明了EnKF算法在降水?dāng)?shù)據(jù)融合方面的有效性。

3.2 降水精度評價指標(biāo)空間分布

箱線圖僅能表達(dá)降水精度評價指標(biāo)誤差的范圍，但并不能直觀展示評價指標(biāo)的空間分布狀況。為此，圖3展示了背景數(shù)據(jù)、再分析數(shù)據(jù)和融合數(shù)據(jù)的精度評價指標(biāo)空間分布圖。由圖3背景數(shù)據(jù)（a）～（c）和融合數(shù)據(jù)圖3（g）可知，在3種背景數(shù)據(jù)各有所長。SI數(shù)據(jù)整體精度高，但空間異質(zhì)性較大，特別是在咸豐、汪營、五峰等流域邊界地區(qū)精度較低，這是因為融合產(chǎn)品質(zhì)量與站網(wǎng)密度呈正相關(guān)關(guān)系［22］；CMORPH各站點精度空間分布較為均勻，但在流域上游地區(qū)的精度低于其他地區(qū)，且在當(dāng)陽壩地區(qū)均方根誤差偏大；TRMM在總體精度和空間分布上均不如前面兩種數(shù)據(jù)，但從平均誤差指標(biāo)來看，流域內(nèi)的部分站點如西流水、茅田等精度高于SI和CMORPH。MSAP綜合了衛(wèi)星數(shù)據(jù)和高精度插值數(shù)據(jù)，改善了背景數(shù)據(jù)流域邊界地區(qū)部分區(qū)域精度低的缺點：在當(dāng)陽壩等站點，MSAP精度高于衛(wèi)星數(shù)據(jù)，與SI精度持平；在五峰和咸豐等站點，MSAP精度高于SI，略低于CMORPH；在汪營等站點，MSAP精度則優(yōu)于3種背景數(shù)據(jù)，證明了EnKF算法的有效性。在圖3（d）～（f）的再分析數(shù)據(jù)中，CMFD與MSWEP各站點精度評價指標(biāo)的空間異質(zhì)性高于MSAP數(shù)據(jù)，兩種數(shù)據(jù)在流域下游的精度高于流域上游。CMFD在所有站點的精度評價指標(biāo)都低于MSAP，MSAP融合了清江流域自動站觀測資料，而CMFD融合了中國氣象局常規(guī)氣象觀測數(shù)據(jù)，MSAP的背景數(shù)據(jù)在站網(wǎng)密度上高于CMFD，這是MSAP精度高于CMFD的重要原因［22］。

3.3 場次強降水量空間分布分析

基于更加直觀地驗證本文建立的融合模型所得到的降水與3種再分析降水產(chǎn)品的空間分布特征，根據(jù)融合后的降水對清江流域1998-1999年間的不同歷時的場次強降水進行分析，以適應(yīng)不同場景下降水產(chǎn)品的應(yīng)用需求。

1998年長江上游從6月中旬開始到8月中旬長達(dá)兩個多月的強對流云團的頻繁發(fā)生是1998年長江上游洪峰多，中游高水位長時間維持的重要原因［23］。這一時期長江上游地區(qū)暴雨洪水頻發(fā)，清江流域自5月中旬以來，也連續(xù)發(fā)生了8次洪峰，其中第四、五、七、八次洪峰與長江干流第一、二、四、六次洪峰不同程度地遭遇［24］。選取四次遭遇中洪峰最大的第四次和第八次進行降水空間分析，其歷時分別為5日和2日，對應(yīng)的流域平均面雨量為186和86 mm。

1998年6月28日-1998年7月2日清江流域發(fā)生了1998年最大的一次強降雨，此次降雨持續(xù)5天，暴雨中心位于流域北部的茅田地區(qū)，單站降水量達(dá)到330 mm，流域平均面雨量達(dá)到了186 mm，此次強降水也引發(fā)了清江流域1998年最大的一次洪水（洪峰流量12 300 m3/s，重現(xiàn)期20年一遇）［24］。圖4給出了各數(shù)據(jù)集對此次降雨的捕捉情況?？梢钥吹礁鲾?shù)據(jù)集都捕捉到了這個全流域強降水的特性，5日累計降水量都在100 mm以上，但在降水空間分布上有幾點差異：CMFD、MSWEP和MSAP能反映出暴雨中心，但是顯示暴雨中心位于流域北部，并向南部遞減，在流域南部上游地區(qū)達(dá)到最低；ERA5趨勢與CMFD趨勢一致，降水量最??；MSWEP與其他數(shù)據(jù)趨勢一致，降水量級為最大；MSAP顯示暴雨中心位于流域西北部，并向東南遞減，在量級上與CMFD保持一致。

清江流域?qū)儆趤啛釒Ъ撅L(fēng)氣候，夏季易受到東南季風(fēng)的影響。1998年8月15日-1998年8月16日清江流域出現(xiàn)了一次強降雨，暴雨中心位于流域東南地區(qū)的五峰站，48 h降水量為146 mm，流域平均面雨量為86 mm。圖5給出了各數(shù)據(jù)集對此次降雨的捕捉情況。從圖5中可以看出，各數(shù)據(jù)集均觀測到了此次降雨的空間分布：ERA5與其他3種數(shù)據(jù)有較大差異，其降雨中心位于流域中上游；CMFD、MSWEP和MSAP的降雨中心均位于流域東南部，并向流域上游遞減，其中CMFD與MSAP的空間分布更相似。

圖5 1998年8月15日-8月16日典型日降水分析Fig.5 Analysis of typical daily precipitation from August 15 to August 16， 1998

4 結(jié) 論

在日尺度上通過集合卡爾曼濾波法（EnKF）融合算法對清江流域1998～1999年的地面插值數(shù)據(jù)及2種衛(wèi)星降水產(chǎn)品進行數(shù)據(jù)融合，得到了清江流域的0.05°×0.05°融合降水序列，利用留一交叉驗證的方式，驗證了融合產(chǎn)品的精度。同時采用了國際上主流的再分析數(shù)據(jù)與本文融合降水?dāng)?shù)據(jù)進行對比，并對清江流域內(nèi)的場次強降水空間分布進行了分析。

（1）交叉驗證結(jié)果顯示，EnKF算法改善了衛(wèi)星降水整體精度低和插值數(shù)據(jù)流域邊界精度低的缺點，使得MSAP數(shù)據(jù)普遍優(yōu)于融合前的單一降水產(chǎn)品；融合數(shù)據(jù)同樣優(yōu)于3種再分析數(shù)據(jù)，CMFD數(shù)據(jù)與融合數(shù)據(jù)最為接近，MSWEP數(shù)據(jù)次之，ERA5數(shù)據(jù)精度最差。

（2）融合產(chǎn)品能準(zhǔn)確捕捉清江流域場次強降水的空間分布特征，MSAP融合數(shù)據(jù)與我國主流融合降水CMFD數(shù)據(jù)呈現(xiàn)相似的空間分布，且強降雨歷時越長，空間分布特征越接近。