基于高斯過程算法的日尺度IMERG降水數(shù)據(jù)與站點數(shù)據(jù)的融合研究
——以湖北省為例

譚偉偉，曾 超，沈煥鋒,3,4，田禮喬*

(1.武漢大學測繪遙感信息國家重點實驗室，武漢 430079；2.武漢大學資源與環(huán)境科學學院，武漢 430079；3.地球空間信息技術(shù)協(xié)同創(chuàng)新中心，武漢 430079；4.地理信息系統(tǒng)教育部重點實驗室，武漢 430079)

作為全球水循環(huán)的主要驅(qū)動力，降水是氣象學、生態(tài)學和水文學的一個重要參量，在物質(zhì)交換和能量平衡中也有重要的作用[1-2].降水數(shù)據(jù)的精度和空間分辨率往往較低,作為水文模型的重要輸入?yún)⒘浚洳淮_定性會導致水文模型的輸出結(jié)果具有很大的不確定性，因此提高降水數(shù)據(jù)的精度和空間分辨率在水文、氣象、生態(tài)等諸多領(lǐng)域具有非常重要的意義[3].

降水數(shù)據(jù)可以從地面站點網(wǎng)、雷達降水估計、衛(wèi)星降水估計等重要的來源中獲取.地面觀測站的降水觀測值僅僅能代表一定范圍內(nèi)的降水特征，具有稀疏和空間分布不均勻的缺點.基于地面站點的降水數(shù)據(jù)插值以獲取更高空間分辨率的面狀降水數(shù)據(jù)，其誤差往往比較大，插值數(shù)據(jù)在地面站點外的區(qū)域的誤差更為顯著[4].天氣雷達能提供高時空分辨率的降水數(shù)據(jù)，但是雷達降水估計的精度受復雜山區(qū)的影響很大.過去幾十年，多個衛(wèi)星降水計劃得到實施，如GPCP[5]、PERSIANN[6]、TRMM[7]、GSMaP[8]、GPM[9]等.一系列降水衛(wèi)星的研發(fā)和相繼升空，使得全球區(qū)域的各個區(qū)域可以獲得降水數(shù)據(jù)源.遙感衛(wèi)星為獲取時空連續(xù)的柵格化降水數(shù)據(jù)提供了一種重要的來源[10].盡管衛(wèi)星降水數(shù)據(jù)豐富了降水數(shù)據(jù)源，但是衛(wèi)星降水數(shù)據(jù)依然有局限性，主要有:受諸多環(huán)境因素影響，衛(wèi)星降水數(shù)據(jù)的精度往往較低，特別是衛(wèi)星降水數(shù)據(jù)的時間分辨率越高，與氣象站點觀測數(shù)據(jù)的一致性往往越低；衛(wèi)星降水數(shù)據(jù)雖然有覆蓋范圍廣泛的優(yōu)勢，但是其空間分辨率較低，無法滿足小流域尺度的研究需求.因此，不同來源的降水數(shù)據(jù)資料都有其優(yōu)缺點，如何綜合多源降水數(shù)據(jù)的優(yōu)點以獲取高精度和高分辨率的柵格化降水數(shù)據(jù)具有重要的意義.

近十多年來，國內(nèi)外有大量關(guān)于降水數(shù)據(jù)融合的研究，降水數(shù)據(jù)融合方法為獲取高精度和高空間分辨率的柵格化降水數(shù)據(jù)提供了可靠的思路.國內(nèi)外提出了多個降水數(shù)據(jù)融合方法，主要包括最優(yōu)插值[11]、卡爾曼濾波[12]、貝葉斯估計[13]、概率密度匹配[14]、小波變換分析[15]、變分方法[16]等.目前關(guān)于最新一代IMERG衛(wèi)星日分辨率降水數(shù)據(jù)的降水融合研究較少.針對日尺度分辨率衛(wèi)星降水數(shù)據(jù)的精度和空間分辨率較低的問題，本文利用2016年7月19日湖北省IMERG日降水數(shù)據(jù)和站點觀測數(shù)據(jù)，采用自適應樣條多元回歸[17]、隨機森林[18]、高斯過程回歸[19]三種算法，提出點面融合估計和站點偏差校正估計兩個融合方案獲取融合的降水數(shù)據(jù)，并對比了兩個融合方案和每個算法的融合效果和精度.

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)概況

湖北省位于中國中部地區(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ū)和氣象站點分布圖.

圖1 湖北省區(qū)域及氣象站點分布Fig.1 Locations of rain gauge stations in Hubei province

1.2 研究數(shù)據(jù)

研究采用降水量豐富的2016年7月19日 的IMERG數(shù)據(jù)、DEM數(shù)據(jù)、Bright temperature(BT)亮溫數(shù)據(jù)和氣象站點實測數(shù)據(jù).IMERG日分辨率降水數(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亮溫數(shù)據(jù)，空間分辨率為500 m×500 m.研究的范圍有55個氣象站點，數(shù)據(jù)來源于中國氣象局國家氣象中心氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/).以上數(shù)據(jù)都使用像元聚合的方法重采樣到1 km×1 km的空間分辨率，該重采樣方法的特點是對采樣窗口內(nèi)部的像元等權(quán)加權(quán)計算重采樣值，能抑制高分辨率輔助數(shù)據(jù)的異常值或像元值的劇烈變化的影響.

2 降水數(shù)據(jù)融合

2.1 高斯過程回歸模型

高斯過程回歸(GPR)在機器學習領(lǐng)域應用比較廣泛，具有嚴格的統(tǒng)計學習理論基礎(chǔ).它在貝葉斯線性回歸方法的基礎(chǔ)上，用核函數(shù)替代貝葉斯線性回歸的線性核，因此高斯過程回歸在對處理高維數(shù)、小樣本、非線性等復雜的問題具有很好的適應性，具有泛化能力強的特點.

高斯過程(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為任意隨機變量.GP定義為f(x)～GP(m(x)，K(x，x′)).

對于高斯過程回歸，假設(shè)

(3)

其中，x為輸入向量，y為受噪聲污染的觀測值，f為GP預測的函數(shù)值，ε為噪聲.

觀測值y的先驗分布為：

(4)

而觀測值y和預測值f*的聯(lián)合先驗分布為：

(5)

由以上可以得出預測值f*的后驗概率密度分布為:

(6)

其中

(7)

cov(f*)=K(x*，x*)-

(8)

高斯過程中選擇的最廣泛的協(xié)方差函數(shù)為平方指數(shù)協(xié)方差函數(shù)：

(9)

文獻[3]提出了一種高時間分辨率的衛(wèi)星降水數(shù)據(jù)和氣象站點數(shù)據(jù)的融合方法，該方法的思路求出衛(wèi)星降水數(shù)據(jù)相對于站點數(shù)據(jù)的偏差，使用高分辨率的輔助數(shù)據(jù)模擬偏差場來求出校正的降水融合數(shù)據(jù).本文稱該融合方法為站點偏差校正估計方法，并提出點面融合估計方法融合IMERG數(shù)據(jù)和站點數(shù)據(jù).本實驗選擇GPR算法和使用兩種融合方案估計研究區(qū)域內(nèi)的日降水量，核函數(shù)采用平方指數(shù)協(xié)方差函數(shù)，并與其他兩個統(tǒng)計學習方法包括自適應多元樣條回歸(MARS)和隨機森林(RF)的實驗結(jié)果作對比.

2.2 點面融合估計方法

假設(shè)pf是融合的降水數(shù)據(jù)，ps是衛(wèi)星降水數(shù)據(jù)，po是站點觀測數(shù)據(jù).點面融合方法估計的步驟如下：

1) 雙線性內(nèi)插0.1°×0.1°分辨率的IMERG降水數(shù)據(jù)，將IMERG數(shù)據(jù)降尺度到1 km×1 km的分辨率，表示為IMERGBi；反距離加權(quán)法內(nèi)插點狀氣象站點數(shù)據(jù)為1 km×1 km分辨率的面狀數(shù)據(jù)，表示為降水空間因子preps.

2) 訓練站點處觀測值(po)與各個預測因子之間的模型，模型表示為

po=f(Lon，Lat，DEM，preps，IMERGBi)，

(10)

其中，Lon、Lat、DEM、preps、IMERGBi分別代表氣象站點處的緯度、經(jīng)度、高程、插值的IMERGE降水值、降水空間因子值，f是以站點處的觀測值與各個因子訓練得到的模型.

3)將第2)步訓練的模型應用到1 km×1 km分辨率的面狀數(shù)據(jù)，得到融合的降水數(shù)據(jù)，表示為

pf=f(Lon，Lat，DEM，preps，IMERGBi，BT)，

(11)

其中，Lon、Lat、DEM、preps、IMERGEBi分別代表相應的1 km×1 km分辨率的面狀數(shù)據(jù)，pf為最終融合的降水數(shù)據(jù).圖2為該方案1的流程圖.

圖2 基于點面融合估計的降水融合方法流程圖Fig.2 Flow chart of precipitation merging method based on point-to-area estimation

2.3 站點偏差校正估計

偏差校正估計的步驟如下：

1)同2.2的(1)，插值得到IMERGBi和preps.

2)假設(shè)融合數(shù)據(jù)pf，衛(wèi)星降水數(shù)據(jù)ps和站點觀測數(shù)據(jù)具有如下關(guān)系：

pf=ps+f(po-ps)，

(12)

式中，po-ps即站點觀測值與衛(wèi)星降水數(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是以站點處的衛(wèi)星降水偏差值與各個因子訓練得到的模型.

3)將第(2)步訓練的模型應用到1 km×1 km分辨率的面狀數(shù)據(jù)，估計整個偏差場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為最終融合的偏差校正降水數(shù)據(jù).圖3為該方案2的流程圖.

圖3 基于偏差校正的降水融合方法流程圖Fig.3 Flow chart of precipitation merging method based Bias correction

3 實驗結(jié)果和討論

3.1 日尺度降水融合結(jié)果

本文選取2016年7月19日的IMERG降水數(shù)據(jù)作為例子，保持1 km×1 km高分辨的亮溫數(shù)據(jù)與衛(wèi)星降水數(shù)據(jù)的時間同步.分別使用點面融合估計方法(方案1)和站點偏差校正估計方法(方案2)對IMERG日降水數(shù)據(jù)與站點觀測數(shù)據(jù)融合，并分別使用GPR、MARS、RF三個統(tǒng)計學習方法估計IMERG日降水數(shù)據(jù)與站點觀測數(shù)據(jù)的融合結(jié)果.方案1和方案2的融合結(jié)果分別如圖4和圖5所示.

圖4(a)表明原始的IMERG降水圖像(圖4(a))空間分辨率較低，氣象站點數(shù)據(jù)直接插值的結(jié)果(圖4(b))相比原始數(shù)據(jù)雖然分辨率得到提升，但是圖像明顯很模糊，缺少細節(jié)變化信息.利用MARS方法通過方案1融合IMERG數(shù)據(jù)和站點數(shù)據(jù)的結(jié)果(圖4(c))的分辨率得到提升，但是圖像細節(jié)信息不豐富；RF融合結(jié)果(圖4(d))有足夠的細節(jié)紋理信息，但是出現(xiàn)了明顯的塊狀，如圖4(d)的最北端的矩形塊狀和圖像中心的圓形塊狀，這明顯不符合降水的分布模式；圖4(e)是基于GPR算法點面融合的降水圖像，融合結(jié)果的圖像分辨率和細節(jié)信息均有較大的提升，空間變化合理，且改善了反距離加權(quán)插值帶來的“牛眼效應”.

使用方案2的流程對IMERG降水數(shù)據(jù)和站點數(shù)據(jù)融合的結(jié)果(圖4(c)～(e))均能提升IMERG降水數(shù)據(jù)的空間分辨率.通過比較兩個方案的融合結(jié)果，圖5(c)相比于圖4(c)，其空間細節(jié)信息得到改善；圖5(d)仍然出現(xiàn)圖4(d)的塊狀現(xiàn)象，采用偏差校正方法不能避免塊狀問題；圖5(e) 明顯優(yōu)于圖5(c)～(d)的融合效果，且相比于圖4(e)變化不大，但是圖5(e)保留了IMERGBi的幾處插值痕跡.

3.2 精度驗證

氣象站點數(shù)據(jù)用于融合結(jié)果的精度驗證，精度驗證的最常用的方法有保留交叉驗證、K折交叉驗證和留一交叉驗證法.研究區(qū)域內(nèi)的站點數(shù)較少，本文選擇留一交叉驗證法，即假設(shè)站點數(shù)為N，每次使用N-1個站點建模并預測剩下的站點處的融合值，該過程循環(huán)N次以保證所有的站點都參與精度評價過程.精度評價選擇的量化指標有相關(guān)系數(shù)R，均方根誤差RMSE和偏差Bias.R表示融合的降水數(shù)據(jù)與站點觀測數(shù)據(jù)的一致性，相關(guān)性越高則R越接近1；RMSE表示評價融合數(shù)據(jù)的整體誤差；Bias表示數(shù)據(jù)的系統(tǒng)誤差，評價融合結(jié)果的偏離程度.三個指標的計算公式如下：

(16)

(17)

圖4 IMERG數(shù)據(jù)(a),站點插值數(shù)據(jù)preps(b),基于MARS(c),RF(d),GPR(e)的點面融合方法的融合數(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)

本文用氣象站點數(shù)據(jù)評估IMERG降水數(shù)據(jù)的原始精度，作為與融合數(shù)據(jù)的精度對比的參照，圖6為精度評價結(jié)果.從圖6可得，該日尺度的精度較低(R=0.43，Bias=0.07%)，無法滿足氣象、水文等領(lǐng)域?qū)Ω呔刃l(wèi)星降水數(shù)據(jù)的需求.

本文以站點降水數(shù)據(jù)為真值，用留一交叉驗證法定量評估兩種融合方案的效果，站點驗證結(jié)果如表1所示，其中GPR_1表示使用方案1的GPR融合結(jié)果，GPR_2表示使用方案2的GPR融合結(jié)果，其余類似.

從模型精度結(jié)果來看，兩種方案的模型訓練的精度(R)都很高，站點處的擬合效果都比較好，因此比較模型的泛化能力更重要.留一交叉驗證結(jié)果表明，對于兩種融合方案，GPR的精度最高，都優(yōu)于RF的融合結(jié)果，而MARS方法表現(xiàn)最差.對比方案1和方案2，方案1的三種算法的融合結(jié)果也都優(yōu)于方案2.因此偏差校正融合方案相比點面融合估計沒有任何優(yōu)勢.對比最佳的降水融合結(jié)果(GPR_1)和原始的IMERG降水數(shù)據(jù)，融合結(jié)果的精度得到較大的提升，R從0.43提升到0.79，Bias從0.07%降低到-0.02%.

圖5 IMERG數(shù)據(jù)(a),站點插值數(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é)果的驗證精度Tab.1 Model performances and validation performances of the above two merging schemes

圖6 2016年7月19日IMERG降水數(shù)據(jù)精度評估結(jié)果圖Fig.6 The accuracy evaluation result of IMERG precipitation data on July 19,2016

4 討論和結(jié)論

4.1 實驗結(jié)果討論

方案2估計衛(wèi)星降水數(shù)據(jù)和氣象站點數(shù)據(jù)的偏差場，本質(zhì)上也是點面融合.方案2相比方案1，增加了偏差計算的步驟，估計的偏差場與插值的IMERGBi數(shù)據(jù)相加，反而使最終的融合結(jié)果受到IMERGBi插值痕跡的影響比較大，最終影響融合結(jié)果的精度.

隨機森林方法本質(zhì)上是一種回歸樹，傳統(tǒng)的回歸樹算法采用節(jié)點分裂機制，即對訓練數(shù)據(jù)集進行劃分，該機制會導致融合結(jié)果在空間上有被劃分的塊狀痕跡.自適應多元樣條回歸采用前向選擇和后向刪除的機制，對訓練數(shù)據(jù)集尤其是對于小樣本數(shù)據(jù)集非常敏感，在樣本稀疏的情況下容易造成非常平滑且有極端值的結(jié)果，因此該方法不適用于融合稀疏分布的站點數(shù)據(jù)與衛(wèi)星降水數(shù)據(jù).

4.2 結(jié)論

本文選擇日分辨率的IMERG衛(wèi)星降水數(shù)據(jù)和站點觀測數(shù)據(jù)兩種降水數(shù)據(jù)資料，并考慮經(jīng)緯度、DEM、亮溫數(shù)據(jù)、插值的IMERG數(shù)據(jù)等作為高分辨率的輔助變量，利用高斯過程回歸算法結(jié)合點面融合估計和站點偏差校正估計兩種方案融合兩種降水數(shù)據(jù)資料，并與隨機森林、自適應樣條多元回歸的結(jié)果作對比.首先對IMERG衛(wèi)星降水數(shù)據(jù)做精度評估，評估結(jié)果表明IMERG數(shù)據(jù)與站點觀測數(shù)據(jù)的不一致較高(R=0.43，Bias=0.07%)，然后從融合圖像的結(jié)果定性分析每個方案的融合效果，最后使用站點觀測數(shù)據(jù)作為真值定量評價融合數(shù)據(jù)的精度.實驗結(jié)果表明：

1) RF融合結(jié)果會出現(xiàn)塊狀痕跡，MARS方法的融合結(jié)果很模糊，而提倡的GPR方法的融合結(jié)果變化合理，且細節(jié)信息優(yōu)于站點數(shù)據(jù)直接插值的結(jié)果.

2) 三種算法的點面融合的精度都優(yōu)于基于偏差校正融合結(jié)果的精度，且避免了偏差校正的融合結(jié)果保留的IMERG數(shù)據(jù)的一些插值痕跡.

3) 使用GPR方法對IMERG衛(wèi)星降水數(shù)據(jù)和站點觀測數(shù)據(jù)做點面融合，融合結(jié)果的精度(R=0.79，Bias=-0.02%)相比原始的衛(wèi)星降水數(shù)據(jù)有較大的提升，該方案可對日分辨率降水數(shù)據(jù)的融合提供一個有價值的參考.

盡管本文提出的基于GPR的點面融合方法能改善IMERG數(shù)據(jù)的質(zhì)量，但是該方法的融合結(jié)果仍然不能完全消除IMERG插值數(shù)據(jù)的痕跡，還需要更深入地研究如何完全消除插值數(shù)據(jù)的影響.