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

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

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

作為全球水循環的主要驅動力，降水是氣象學、生態學和水文學的一個重要參量，在物質交換和能量平衡中也有重要的作用[1-2].降水數據的精度和空間分辨率往往較低,作為水文模型的重要輸入參量，其不確定性會導致水文模型的輸出結果具有很大的不確定性，因此提高降水數據的精度和空間分辨率在水文、氣象、生態等諸多領域具有非常重要的意義[3].

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

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

1 研究區域與數據

1.1 研究區概況

湖北省位于中國中部地區，介于29°01′53″～33°16′47″N，108°21′42″～116°07′50″E，總面積約18.59萬km2.區域處于地勢第二級階梯向第三級階梯過渡的地帶，地貌類型豐富，包括山地、丘陵、崗地和平原，山地和丘陵約占全省面積80%，地勢高低懸殊.湖北省的地形對降水影響較大，降水地域分布呈由南向北遞減趨勢，且有季節變化規律，一般是夏季最多，冬季最少.6月中旬至7月中旬是梅雨期，降水量最大.圖1為研究區和氣象站點分布圖.

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

1.2 研究數據

研究采用降水量豐富的2016年7月19日 的IMERG數據、DEM數據、Bright temperature(BT)亮溫數據和氣象站點實測數據.IMERG日分辨率降水數據下載自NASA(https://pmm.nasa.gov/data-access/downloads/gpm)，數據覆蓋范圍為90°S～90°N，空間分辨率為0.1°×0.1°(約10 km).DEM數據來源于地理空間數據云(http://www.gscloud.cn/).BT數據來自NASA的MYDTBGA亮溫數據，空間分辨率為500 m×500 m.研究的范圍有55個氣象站點，數據來源于中國氣象局國家氣象中心氣象數據網(http://data.cma.cn/).以上數據都使用像元聚合的方法重采樣到1 km×1 km的空間分辨率，該重采樣方法的特點是對采樣窗口內部的像元等權加權計算重采樣值，能抑制高分辨率輔助數據的異常值或像元值的劇烈變化的影響.

2 降水數據融合

2.1 高斯過程回歸模型

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

高斯過程(GP)從函數空間角度上描述函數分布，其性質由均值函數m(x)和協方差函數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′)).

對于高斯過程回歸，假設

(3)

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

觀測值y的先驗分布為：

(4)

而觀測值y和預測值f*的聯合先驗分布為：

(5)

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

(6)

其中

(7)

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

(8)

高斯過程中選擇的最廣泛的協方差函數為平方指數協方差函數：

(9)

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

2.2 點面融合估計方法

假設pf是融合的降水數據，ps是衛星降水數據，po是站點觀測數據.點面融合方法估計的步驟如下：

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

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

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

(10)

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

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

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

(11)

其中，Lon、Lat、DEM、preps、IMERGEBi分別代表相應的1 km×1 km分辨率的面狀數據，pf為最終融合的降水數據.圖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)假設融合數據pf，衛星降水數據ps和站點觀測數據具有如下關系：

pf=ps+f(po-ps)，

(12)

式中，po-ps即站點觀測值與衛星降水數據的差，代表衛星降水的偏差.f(po-ps)表示為：

f(po-ps)=f(Lon，Lat，DEM，

preps，IMERGBi，BT)，

(13)

其中，Lon、Lat、DEM、preps、IMERGBi同2.2的第2)步，f是以站點處的衛星降水偏差值與各個因子訓練得到的模型.

3)將第(2)步訓練的模型應用到1 km×1 km分辨率的面狀數據，估計整個偏差場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)相加得到偏差校正結果pf，表示為：

pf=ps+f′(po-ps).

(15)

式中，pf為最終融合的偏差校正降水數據.圖3為該方案2的流程圖.

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

3 實驗結果和討論

3.1 日尺度降水融合結果

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

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

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

3.2 精度驗證

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

(16)

(17)

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

本文用氣象站點數據評估IMERG降水數據的原始精度，作為與融合數據的精度對比的參照，圖6為精度評價結果.從圖6可得，該日尺度的精度較低(R=0.43，Bias=0.07%)，無法滿足氣象、水文等領域對高精度衛星降水數據的需求.

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

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

圖5 IMERG數據(a),站點插值數據preps(b),基于MARS(c),RF(d),GPR(e)的偏差校正方法的融合數據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 兩種融合方案的模型精度和降水融合結果的驗證精度Tab.1 Model performances and validation performances of the above two merging schemes

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

4 討論和結論

4.1 實驗結果討論

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

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

4.2 結論

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

1) RF融合結果會出現塊狀痕跡，MARS方法的融合結果很模糊，而提倡的GPR方法的融合結果變化合理，且細節信息優于站點數據直接插值的結果.

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

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

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