譚偉偉,曾 超,沈煥鋒,3,4,田禮喬*
(1.武漢大學(xué)測繪遙感信息國家重點(diǎn)實(shí)驗(yàn)室,武漢 430079;2.武漢大學(xué)資源與環(huán)境科學(xué)學(xué)院,武漢 430079;3.地球空間信息技術(shù)協(xié)同創(chuàng)新中心,武漢 430079;4.地理信息系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室,武漢 430079)
作為全球水循環(huán)的主要驅(qū)動力,降水是氣象學(xué)、生態(tài)學(xué)和水文學(xué)的一個(gè)重要參量,在物質(zhì)交換和能量平衡中也有重要的作用[1-2].降水?dāng)?shù)據(jù)的精度和空間分辨率往往較低,作為水文模型的重要輸入?yún)⒘?,其不確定性會導(dǎo)致水文模型的輸出結(jié)果具有很大的不確定性,因此提高降水?dāng)?shù)據(jù)的精度和空間分辨率在水文、氣象、生態(tài)等諸多領(lǐng)域具有非常重要的意義[3].
降水?dāng)?shù)據(jù)可以從地面站點(diǎn)網(wǎng)、雷達(dá)降水估計(jì)、衛(wèi)星降水估計(jì)等重要的來源中獲取.地面觀測站的降水觀測值僅僅能代表一定范圍內(nèi)的降水特征,具有稀疏和空間分布不均勻的缺點(diǎn).基于地面站點(diǎn)的降水?dāng)?shù)據(jù)插值以獲取更高空間分辨率的面狀降水?dāng)?shù)據(jù),其誤差往往比較大,插值數(shù)據(jù)在地面站點(diǎn)外的區(qū)域的誤差更為顯著[4].天氣雷達(dá)能提供高時(shí)空分辨率的降水?dāng)?shù)據(jù),但是雷達(dá)降水估計(jì)的精度受復(fù)雜山區(qū)的影響很大.過去幾十年,多個(gè)衛(wèi)星降水計(jì)劃得到實(shí)施,如GPCP[5]、PERSIANN[6]、TRMM[7]、GSMaP[8]、GPM[9]等.一系列降水衛(wèi)星的研發(fā)和相繼升空,使得全球區(qū)域的各個(gè)區(qū)域可以獲得降水?dāng)?shù)據(jù)源.遙感衛(wèi)星為獲取時(shí)空連續(xù)的柵格化降水?dāng)?shù)據(jù)提供了一種重要的來源[10].盡管衛(wèi)星降水?dāng)?shù)據(jù)豐富了降水?dāng)?shù)據(jù)源,但是衛(wèi)星降水?dāng)?shù)據(jù)依然有局限性,主要有:受諸多環(huán)境因素影響,衛(wèi)星降水?dāng)?shù)據(jù)的精度往往較低,特別是衛(wèi)星降水?dāng)?shù)據(jù)的時(shí)間分辨率越高,與氣象站點(diǎn)觀測數(shù)據(jù)的一致性往往越低;衛(wèi)星降水?dāng)?shù)據(jù)雖然有覆蓋范圍廣泛的優(yōu)勢,但是其空間分辨率較低,無法滿足小流域尺度的研究需求.因此,不同來源的降水?dāng)?shù)據(jù)資料都有其優(yōu)缺點(diǎn),如何綜合多源降水?dāng)?shù)據(jù)的優(yōu)點(diǎn)以獲取高精度和高分辨率的柵格化降水?dāng)?shù)據(jù)具有重要的意義.
近十多年來,國內(nèi)外有大量關(guān)于降水?dāng)?shù)據(jù)融合的研究,降水?dāng)?shù)據(jù)融合方法為獲取高精度和高空間分辨率的柵格化降水?dāng)?shù)據(jù)提供了可靠的思路.國內(nèi)外提出了多個(gè)降水?dāng)?shù)據(jù)融合方法,主要包括最優(yōu)插值[11]、卡爾曼濾波[12]、貝葉斯估計(jì)[13]、概率密度匹配[14]、小波變換分析[15]、變分方法[16]等.目前關(guān)于最新一代IMERG衛(wèi)星日分辨率降水?dāng)?shù)據(jù)的降水融合研究較少.針對日尺度分辨率衛(wèi)星降水?dāng)?shù)據(jù)的精度和空間分辨率較低的問題,本文利用2016年7月19日湖北省IMERG日降水?dāng)?shù)據(jù)和站點(diǎn)觀測數(shù)據(jù),采用自適應(yīng)樣條多元回歸[17]、隨機(jī)森林[18]、高斯過程回歸[19]三種算法,提出點(diǎn)面融合估計(jì)和站點(diǎn)偏差校正估計(jì)兩個(gè)融合方案獲取融合的降水?dāng)?shù)據(jù),并對比了兩個(gè)融合方案和每個(gè)算法的融合效果和精度.
湖北省位于中國中部地區(qū),介于29°01′53″~33°16′47″N,108°21′42″~116°07′50″E,總面積約18.59萬km2.區(qū)域處于地勢第二級階梯向第三級階梯過渡的地帶,地貌類型豐富,包括山地、丘陵、崗地和平原,山地和丘陵約占全省面積80%,地勢高低懸殊.湖北省的地形對降水影響較大,降水地域分布呈由南向北遞減趨勢,且有季節(jié)變化規(guī)律,一般是夏季最多,冬季最少.6月中旬至7月中旬是梅雨期,降水量最大.圖1為研究區(qū)和氣象站點(diǎn)分布圖.
圖1 湖北省區(qū)域及氣象站點(diǎn)分布Fig.1 Locations of rain gauge stations in Hubei province
研究采用降水量豐富的2016年7月19日 的IMERG數(shù)據(jù)、DEM數(shù)據(jù)、Bright temperature(BT)亮溫?cái)?shù)據(jù)和氣象站點(diǎn)實(shí)測數(shù)據(jù).IMERG日分辨率降水?dāng)?shù)據(jù)下載自NASA(https://pmm.nasa.gov/data-access/downloads/gpm),數(shù)據(jù)覆蓋范圍為90°S~90°N,空間分辨率為0.1°×0.1°(約10 km).DEM數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http://www.gscloud.cn/).BT數(shù)據(jù)來自NASA的MYDTBGA亮溫?cái)?shù)據(jù),空間分辨率為500 m×500 m.研究的范圍有55個(gè)氣象站點(diǎn),數(shù)據(jù)來源于中國氣象局國家氣象中心氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/).以上數(shù)據(jù)都使用像元聚合的方法重采樣到1 km×1 km的空間分辨率,該重采樣方法的特點(diǎn)是對采樣窗口內(nèi)部的像元等權(quán)加權(quán)計(jì)算重采樣值,能抑制高分辨率輔助數(shù)據(jù)的異常值或像元值的劇烈變化的影響.
高斯過程回歸(GPR)在機(jī)器學(xué)習(xí)領(lǐng)域應(yīng)用比較廣泛,具有嚴(yán)格的統(tǒng)計(jì)學(xué)習(xí)理論基礎(chǔ).它在貝葉斯線性回歸方法的基礎(chǔ)上,用核函數(shù)替代貝葉斯線性回歸的線性核,因此高斯過程回歸在對處理高維數(shù)、小樣本、非線性等復(fù)雜的問題具有很好的適應(yīng)性,具有泛化能力強(qiáng)的特點(diǎn).
高斯過程(GP)從函數(shù)空間角度上描述函數(shù)分布,其性質(zhì)由均值函數(shù)m(x)和協(xié)方差函數(shù)K(x,x′)決定,形式分別為:
m(x)=E[f(x)],
(1)
K(x,x′)=E[(f(x)-m(x))(f(x′)-m(x′))],
(2)
其中x、x′屬于Rd為任意隨機(jī)變量.GP定義為f(x)~GP(m(x),K(x,x′)).
對于高斯過程回歸,假設(shè)
(3)
其中,x為輸入向量,y為受噪聲污染的觀測值,f為GP預(yù)測的函數(shù)值,ε為噪聲.
觀測值y的先驗(yàn)分布為:
(4)
而觀測值y和預(yù)測值f*的聯(lián)合先驗(yàn)分布為:
(5)
由以上可以得出預(yù)測值f*的后驗(yàn)概率密度分布為:
(6)
其中
(7)
cov(f*)=K(x*,x*)-
(8)
高斯過程中選擇的最廣泛的協(xié)方差函數(shù)為平方指數(shù)協(xié)方差函數(shù):
(9)
文獻(xiàn)[3]提出了一種高時(shí)間分辨率的衛(wèi)星降水?dāng)?shù)據(jù)和氣象站點(diǎn)數(shù)據(jù)的融合方法,該方法的思路求出衛(wèi)星降水?dāng)?shù)據(jù)相對于站點(diǎn)數(shù)據(jù)的偏差,使用高分辨率的輔助數(shù)據(jù)模擬偏差場來求出校正的降水融合數(shù)據(jù).本文稱該融合方法為站點(diǎn)偏差校正估計(jì)方法,并提出點(diǎn)面融合估計(jì)方法融合IMERG數(shù)據(jù)和站點(diǎn)數(shù)據(jù).本實(shí)驗(yàn)選擇GPR算法和使用兩種融合方案估計(jì)研究區(qū)域內(nèi)的日降水量,核函數(shù)采用平方指數(shù)協(xié)方差函數(shù),并與其他兩個(gè)統(tǒng)計(jì)學(xué)習(xí)方法包括自適應(yīng)多元樣條回歸(MARS)和隨機(jī)森林(RF)的實(shí)驗(yàn)結(jié)果作對比.
假設(shè)pf是融合的降水?dāng)?shù)據(jù),ps是衛(wèi)星降水?dāng)?shù)據(jù),po是站點(diǎn)觀測數(shù)據(jù).點(diǎn)面融合方法估計(jì)的步驟如下:
1) 雙線性內(nèi)插0.1°×0.1°分辨率的IMERG降水?dāng)?shù)據(jù),將IMERG數(shù)據(jù)降尺度到1 km×1 km的分辨率,表示為IMERGBi;反距離加權(quán)法內(nèi)插點(diǎn)狀氣象站點(diǎn)數(shù)據(jù)為1 km×1 km分辨率的面狀數(shù)據(jù),表示為降水空間因子preps.
2) 訓(xùn)練站點(diǎn)處觀測值(po)與各個(gè)預(yù)測因子之間的模型,模型表示為
po=f(Lon,Lat,DEM,preps,IMERGBi),
(10)
其中,Lon、Lat、DEM、preps、IMERGBi分別代表氣象站點(diǎn)處的緯度、經(jīng)度、高程、插值的IMERGE降水值、降水空間因子值,f是以站點(diǎn)處的觀測值與各個(gè)因子訓(xùn)練得到的模型.
3)將第2)步訓(xùn)練的模型應(yīng)用到1 km×1 km分辨率的面狀數(shù)據(jù),得到融合的降水?dāng)?shù)據(jù),表示為
pf=f(Lon,Lat,DEM,preps,IMERGBi,BT),
(11)
其中,Lon、Lat、DEM、preps、IMERGEBi分別代表相應(yīng)的1 km×1 km分辨率的面狀數(shù)據(jù),pf為最終融合的降水?dāng)?shù)據(jù).圖2為該方案1的流程圖.
圖2 基于點(diǎn)面融合估計(jì)的降水融合方法流程圖Fig.2 Flow chart of precipitation merging method based on point-to-area estimation
偏差校正估計(jì)的步驟如下:
1)同2.2的(1),插值得到IMERGBi和preps.
2)假設(shè)融合數(shù)據(jù)pf,衛(wèi)星降水?dāng)?shù)據(jù)ps和站點(diǎn)觀測數(shù)據(jù)具有如下關(guān)系:
pf=ps+f(po-ps),
(12)
式中,po-ps即站點(diǎn)觀測值與衛(wèi)星降水?dāng)?shù)據(jù)的差,代表衛(wèi)星降水的偏差.f(po-ps)表示為:
f(po-ps)=f(Lon,Lat,DEM,
preps,IMERGBi,BT),
(13)
其中,Lon、Lat、DEM、preps、IMERGBi同2.2的第2)步,f是以站點(diǎn)處的衛(wèi)星降水偏差值與各個(gè)因子訓(xùn)練得到的模型.
3)將第(2)步訓(xùn)練的模型應(yīng)用到1 km×1 km分辨率的面狀數(shù)據(jù),估計(jì)整個(gè)偏差場f′(po-ps),表示為
f′(po-ps)=f(Lon,Lat,DEM,
preps,IMERGBi,BT),
(14)
其中,Lon、Lat、DEM、preps、IMERGBi同2.2的第3)步.
4)將插值的IMERGBi(ps)與偏差場f′(po-ps)相加得到偏差校正結(jié)果pf,表示為:
pf=ps+f′(po-ps).
(15)
式中,pf為最終融合的偏差校正降水?dāng)?shù)據(jù).圖3為該方案2的流程圖.
圖3 基于偏差校正的降水融合方法流程圖Fig.3 Flow chart of precipitation merging method based Bias correction
本文選取2016年7月19日的IMERG降水?dāng)?shù)據(jù)作為例子,保持1 km×1 km高分辨的亮溫?cái)?shù)據(jù)與衛(wèi)星降水?dāng)?shù)據(jù)的時(shí)間同步.分別使用點(diǎn)面融合估計(jì)方法(方案1)和站點(diǎn)偏差校正估計(jì)方法(方案2)對IMERG日降水?dāng)?shù)據(jù)與站點(diǎn)觀測數(shù)據(jù)融合,并分別使用GPR、MARS、RF三個(gè)統(tǒng)計(jì)學(xué)習(xí)方法估計(jì)IMERG日降水?dāng)?shù)據(jù)與站點(diǎn)觀測數(shù)據(jù)的融合結(jié)果.方案1和方案2的融合結(jié)果分別如圖4和圖5所示.
圖4(a)表明原始的IMERG降水圖像(圖4(a))空間分辨率較低,氣象站點(diǎn)數(shù)據(jù)直接插值的結(jié)果(圖4(b))相比原始數(shù)據(jù)雖然分辨率得到提升,但是圖像明顯很模糊,缺少細(xì)節(jié)變化信息.利用MARS方法通過方案1融合IMERG數(shù)據(jù)和站點(diǎn)數(shù)據(jù)的結(jié)果(圖4(c))的分辨率得到提升,但是圖像細(xì)節(jié)信息不豐富;RF融合結(jié)果(圖4(d))有足夠的細(xì)節(jié)紋理信息,但是出現(xiàn)了明顯的塊狀,如圖4(d)的最北端的矩形塊狀和圖像中心的圓形塊狀,這明顯不符合降水的分布模式;圖4(e)是基于GPR算法點(diǎn)面融合的降水圖像,融合結(jié)果的圖像分辨率和細(xì)節(jié)信息均有較大的提升,空間變化合理,且改善了反距離加權(quán)插值帶來的“牛眼效應(yīng)”.
使用方案2的流程對IMERG降水?dāng)?shù)據(jù)和站點(diǎn)數(shù)據(jù)融合的結(jié)果(圖4(c)~(e))均能提升IMERG降水?dāng)?shù)據(jù)的空間分辨率.通過比較兩個(gè)方案的融合結(jié)果,圖5(c)相比于圖4(c),其空間細(xì)節(jié)信息得到改善;圖5(d)仍然出現(xiàn)圖4(d)的塊狀現(xiàn)象,采用偏差校正方法不能避免塊狀問題;圖5(e) 明顯優(yōu)于圖5(c)~(d)的融合效果,且相比于圖4(e)變化不大,但是圖5(e)保留了IMERGBi的幾處插值痕跡.
氣象站點(diǎn)數(shù)據(jù)用于融合結(jié)果的精度驗(yàn)證,精度驗(yàn)證的最常用的方法有保留交叉驗(yàn)證、K折交叉驗(yàn)證和留一交叉驗(yàn)證法.研究區(qū)域內(nèi)的站點(diǎn)數(shù)較少,本文選擇留一交叉驗(yàn)證法,即假設(shè)站點(diǎn)數(shù)為N,每次使用N-1個(gè)站點(diǎn)建模并預(yù)測剩下的站點(diǎn)處的融合值,該過程循環(huán)N次以保證所有的站點(diǎn)都參與精度評價(jià)過程.精度評價(jià)選擇的量化指標(biāo)有相關(guān)系數(shù)R,均方根誤差RMSE和偏差Bias.R表示融合的降水?dāng)?shù)據(jù)與站點(diǎn)觀測數(shù)據(jù)的一致性,相關(guān)性越高則R越接近1;RMSE表示評價(jià)融合數(shù)據(jù)的整體誤差;Bias表示數(shù)據(jù)的系統(tǒng)誤差,評價(jià)融合結(jié)果的偏離程度.三個(gè)指標(biāo)的計(jì)算公式如下:
(16)
(17)
圖4 IMERG數(shù)據(jù)(a),站點(diǎn)插值數(shù)據(jù)preps(b),基于MARS(c),RF(d),GPR(e)的點(diǎn)面融合方法的融合數(shù)據(jù)Fig.4 Comparison of (a) IMERGE precipitation,(b) preps,the merged results by point-to-surface fusion method using (c) MARS,(d) RF,and (e) GPR
(18)
本文用氣象站點(diǎn)數(shù)據(jù)評估IMERG降水?dāng)?shù)據(jù)的原始精度,作為與融合數(shù)據(jù)的精度對比的參照,圖6為精度評價(jià)結(jié)果.從圖6可得,該日尺度的精度較低(R=0.43,Bias=0.07%),無法滿足氣象、水文等領(lǐng)域?qū)Ω呔刃l(wèi)星降水?dāng)?shù)據(jù)的需求.
本文以站點(diǎn)降水?dāng)?shù)據(jù)為真值,用留一交叉驗(yàn)證法定量評估兩種融合方案的效果,站點(diǎn)驗(yàn)證結(jié)果如表1所示,其中GPR_1表示使用方案1的GPR融合結(jié)果,GPR_2表示使用方案2的GPR融合結(jié)果,其余類似.
從模型精度結(jié)果來看,兩種方案的模型訓(xùn)練的精度(R)都很高,站點(diǎn)處的擬合效果都比較好,因此比較模型的泛化能力更重要.留一交叉驗(yàn)證結(jié)果表明,對于兩種融合方案,GPR的精度最高,都優(yōu)于RF的融合結(jié)果,而MARS方法表現(xiàn)最差.對比方案1和方案2,方案1的三種算法的融合結(jié)果也都優(yōu)于方案2.因此偏差校正融合方案相比點(diǎn)面融合估計(jì)沒有任何優(yōu)勢.對比最佳的降水融合結(jié)果(GPR_1)和原始的IMERG降水?dāng)?shù)據(jù),融合結(jié)果的精度得到較大的提升,R從0.43提升到0.79,Bias從0.07%降低到-0.02%.
圖5 IMERG數(shù)據(jù)(a),站點(diǎn)插值數(shù)據(jù)preps(b),基于MARS(c),RF(d),GPR(e)的偏差校正方法的融合數(shù)據(jù)Fig.5 Comparison of (a) IMERG precipitation,(b) preps,the merged result by bias-correction fusion method using (c)MARS,(d) RF,and (e) GPR
表1 兩種融合方案的模型精度和降水融合結(jié)果的驗(yàn)證精度Tab.1 Model performances and validation performances of the above two merging schemes
圖6 2016年7月19日IMERG降水?dāng)?shù)據(jù)精度評估結(jié)果圖Fig.6 The accuracy evaluation result of IMERG precipitation data on July 19,2016
方案2估計(jì)衛(wèi)星降水?dāng)?shù)據(jù)和氣象站點(diǎn)數(shù)據(jù)的偏差場,本質(zhì)上也是點(diǎn)面融合.方案2相比方案1,增加了偏差計(jì)算的步驟,估計(jì)的偏差場與插值的IMERGBi數(shù)據(jù)相加,反而使最終的融合結(jié)果受到IMERGBi插值痕跡的影響比較大,最終影響融合結(jié)果的精度.
隨機(jī)森林方法本質(zhì)上是一種回歸樹,傳統(tǒng)的回歸樹算法采用節(jié)點(diǎn)分裂機(jī)制,即對訓(xùn)練數(shù)據(jù)集進(jìn)行劃分,該機(jī)制會導(dǎo)致融合結(jié)果在空間上有被劃分的塊狀痕跡.自適應(yīng)多元樣條回歸采用前向選擇和后向刪除的機(jī)制,對訓(xùn)練數(shù)據(jù)集尤其是對于小樣本數(shù)據(jù)集非常敏感,在樣本稀疏的情況下容易造成非常平滑且有極端值的結(jié)果,因此該方法不適用于融合稀疏分布的站點(diǎn)數(shù)據(jù)與衛(wèi)星降水?dāng)?shù)據(jù).
本文選擇日分辨率的IMERG衛(wèi)星降水?dāng)?shù)據(jù)和站點(diǎn)觀測數(shù)據(jù)兩種降水?dāng)?shù)據(jù)資料,并考慮經(jīng)緯度、DEM、亮溫?cái)?shù)據(jù)、插值的IMERG數(shù)據(jù)等作為高分辨率的輔助變量,利用高斯過程回歸算法結(jié)合點(diǎn)面融合估計(jì)和站點(diǎn)偏差校正估計(jì)兩種方案融合兩種降水?dāng)?shù)據(jù)資料,并與隨機(jī)森林、自適應(yīng)樣條多元回歸的結(jié)果作對比.首先對IMERG衛(wèi)星降水?dāng)?shù)據(jù)做精度評估,評估結(jié)果表明IMERG數(shù)據(jù)與站點(diǎn)觀測數(shù)據(jù)的不一致較高(R=0.43,Bias=0.07%),然后從融合圖像的結(jié)果定性分析每個(gè)方案的融合效果,最后使用站點(diǎn)觀測數(shù)據(jù)作為真值定量評價(jià)融合數(shù)據(jù)的精度.實(shí)驗(yàn)結(jié)果表明:
1) RF融合結(jié)果會出現(xiàn)塊狀痕跡,MARS方法的融合結(jié)果很模糊,而提倡的GPR方法的融合結(jié)果變化合理,且細(xì)節(jié)信息優(yōu)于站點(diǎn)數(shù)據(jù)直接插值的結(jié)果.
2) 三種算法的點(diǎn)面融合的精度都優(yōu)于基于偏差校正融合結(jié)果的精度,且避免了偏差校正的融合結(jié)果保留的IMERG數(shù)據(jù)的一些插值痕跡.
3) 使用GPR方法對IMERG衛(wèi)星降水?dāng)?shù)據(jù)和站點(diǎn)觀測數(shù)據(jù)做點(diǎn)面融合,融合結(jié)果的精度(R=0.79,Bias=-0.02%)相比原始的衛(wèi)星降水?dāng)?shù)據(jù)有較大的提升,該方案可對日分辨率降水?dāng)?shù)據(jù)的融合提供一個(gè)有價(jià)值的參考.
盡管本文提出的基于GPR的點(diǎn)面融合方法能改善IMERG數(shù)據(jù)的質(zhì)量,但是該方法的融合結(jié)果仍然不能完全消除IMERG插值數(shù)據(jù)的痕跡,還需要更深入地研究如何完全消除插值數(shù)據(jù)的影響.