基于深度學(xué)習(xí)的多源降水?dāng)?shù)據(jù)融合方法及其應(yīng)用

覃 曉 東,朱 仟,周 東 旸,丁 潔

(東南大學(xué) 土木工程學(xué)院,江蘇 南京 211189)

0 引 言

降水是水文循環(huán)中的重要組成部分,同時也是水文模擬中的主要輸入變量,降水估計的準確性直接影響水文模擬的精度[1]。以往的研究表明,水文模擬中70%～80%的不確定性源于降水?dāng)?shù)據(jù)[2-3],因此,亟需獲取高精度、高分辨率的實時降水?dāng)?shù)據(jù)為水文模擬和洪澇災(zāi)害預(yù)報提供支撐。

傳統(tǒng)的地面降水觀測數(shù)據(jù)是最為直接準確的數(shù)據(jù)來源,但是受到成本、地形等因素的限制,地面觀測站點的數(shù)量以及覆蓋范圍較為有限,分布不均,難以充分反映降水的空間分布。衛(wèi)星遙感經(jīng)過多年的發(fā)展,已能夠提供全天候、全覆蓋的降水?dāng)?shù)據(jù)產(chǎn)品。已有較多研究對遙感降水產(chǎn)品的精度及水文方面的應(yīng)用進行了評估[4-6],并表明遙感降水產(chǎn)品能夠比較準確地反映降水的時空分布特征,但在精度方面仍有較大的改進空間。

自20 世紀80 年代數(shù)據(jù)融合首次應(yīng)用于雷達-觀測降水?dāng)?shù)據(jù)融合以來[7],已逐漸成為改善降水?dāng)?shù)據(jù)分辨率和質(zhì)量的主流方法之一。Chao等[8]運用混合地理加權(quán)回歸的方法,對地面觀測降水?dāng)?shù)據(jù)和CMORPH遙感降水?dāng)?shù)據(jù)進行了融合,顯著提高了遙感降水?dāng)?shù)據(jù)的空間分辨率,并改善了其對降水空間異質(zhì)性的捕捉。Chen等[9]使用兩個神經(jīng)網(wǎng)絡(luò)模型,設(shè)計了一個降水?dāng)?shù)據(jù)融合系統(tǒng),其利用地面雷達降水估計作為橋梁,對地面觀測降水和遙感降水進行融合,有效提高了遙感降水產(chǎn)品的精度。Wu等[10]則通過深度學(xué)習(xí)方法,融合了地面觀測降水?dāng)?shù)據(jù)和TRMM遙感降水?dāng)?shù)據(jù),在中國范圍內(nèi)生成了一套準確度和空間分辨率更高的降雨數(shù)據(jù)集。以上研究都表明降水?dāng)?shù)據(jù)融合產(chǎn)品能夠充分結(jié)合地面觀測降水?dāng)?shù)據(jù)的精度以及遙感降水?dāng)?shù)據(jù)的高分辨率,而事實上其質(zhì)量很大程度上取決于融合算法的選擇和設(shè)計。本次研究采用基于長短期記憶神經(jīng)網(wǎng)絡(luò)(Long Short-Term Memory network,LSTM)的空間插值方法和時間動態(tài)貝葉斯模型平均(Dynamic Bayesian Model Averaging,DBMA)數(shù)據(jù)融合方法,對地面觀測降水?dāng)?shù)據(jù)和基于綜合多衛(wèi)星反演的全球降水測量(Integrated Multi-satellite Retrievals for Global Precipitation Measurement,IMERG)遙感降水?dāng)?shù)據(jù)進行融合,旨在獲得更高精度、高分辨率的降水產(chǎn)品,并提高其在水文模擬中的應(yīng)用效果。

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

1.1 研究區(qū)域

湘江是長江最大的支流之一,也是湖南省最長的河流,全長約856 km。湘江流域位于長江中下游地區(qū)(東經(jīng)110.50°～114.25°,北緯24.50°～28.25°),湘潭站以上流域面積約為82 375 km2(見圖1)。流域?qū)賮啛釒Ъ撅L(fēng)氣候,年均氣溫約為17 ℃,年均降水量約為1400～1700 mm。流域地形以山地丘陵為主,地勢上總體呈現(xiàn)南高北低的特點。流域內(nèi)的土地利用類型以林地為主,西部和南部分布著大量農(nóng)業(yè)用地[11],是中國重要的糧食基地。而流域下游城鎮(zhèn)密集人口眾多,是湖南省的經(jīng)濟發(fā)展中心。

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

本次研究選用美國國家航空航天局NASA (National Aeronautics and Space Administration,NASA) 發(fā)布的IMERG V05B作為遙感降水?dāng)?shù)據(jù)來源[12]。IMERG由多種降雨反演算法組成,并結(jié)合了從GPM傳感器獲得的多個遙感數(shù)據(jù)源,能夠以0.1°的空間分辨率和30 min的時間分辨率提供降水估算。本次研究選用IMERG系統(tǒng)產(chǎn)生的滯后時間為4 h的IMERG Early Run(以下簡稱IMERG-E)作為降水?dāng)?shù)據(jù)融合的遙感降水?dāng)?shù)據(jù)來源。

從中國國家氣象信息中心(https:∥data.cma.cn)獲取湘江流域27個氣象站的日降水?dāng)?shù)據(jù)(以下簡稱GAUGE)作為降水?dāng)?shù)據(jù)融合的地面站點觀測降水?dāng)?shù)據(jù)來源。同時,應(yīng)用27個氣象站的其他氣象數(shù)據(jù)(最高和最低溫度、相對濕度、風(fēng)速和太陽輻射)獲得用于驅(qū)動水文模型的氣象變量。

為了對融合降水?dāng)?shù)據(jù)進行評估,本次研究選用中國國家氣象局發(fā)布的中國自動站與CMORPH降水產(chǎn)品融合的逐時降水量網(wǎng)格數(shù)據(jù)集1.0版(以下簡稱CMA)作為參考降水?dāng)?shù)據(jù)集[13]。該產(chǎn)品具有較小的偏差,可以捕捉到強降水期間的降水時空變化特征,在中國區(qū)域優(yōu)于國際同類型產(chǎn)品的精度[14]。研究中使用的湘潭水文站2014～2017年間每日觀測徑流數(shù)據(jù)來自湖南省水文局,使用的2015年歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)數(shù)據(jù)來自中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心(http:∥www.resdc.cn/)。

2 研究方法

本次研究選用湘江流域27個地面站點處的GAUGE降水?dāng)?shù)據(jù)和IMERG-E遙感降水?dāng)?shù)據(jù),在0.1°的空間分辨率下,對2014年4月1日至2017年12月31日的日降水?dāng)?shù)據(jù)進行數(shù)據(jù)融合。首先采用基于LSTM的空間插值方法,對27個地面站點處的GAUGE降水?dāng)?shù)據(jù)進行空間插值,獲得0.1°空間分辨率下的降水估計(以下簡稱GAUGE-grid)。然后通過DBMA數(shù)據(jù)融合方法,對GAUGE-grid和IMERG-E在時間序列上進行數(shù)據(jù)融合,以獲得整個流域上基于遙感降水的高精度降水估計(以下簡稱CMA-E)。最后對融合降水產(chǎn)品的精度及水文應(yīng)用方面進行評估。圖2展示了上述的時空動態(tài)數(shù)據(jù)融合方法。

圖2 時空動態(tài)數(shù)據(jù)融合方法Fig.2 Spatio-temporal dynamic data fusion framework

2.1 基于長短期記憶(LSTM)神經(jīng)網(wǎng)絡(luò)的空間插值方法

本次研究基于LSTM對降水?dāng)?shù)據(jù)進行空間插值,該過程主要包括建立時間序列樣本和模型訓(xùn)練兩部分。假設(shè)在研究流域有N個地面觀測站點,這N個站點的時間序列樣本則為

P=[P1,P2,…,PN]∈R1×N

(1)

Lat=[lat1,lat2,…,latN]∈R1×N

(2)

Lon=[lon1,lon2,…,lonN]∈R1×N

(3)

H=[H1,H2,…,HN]∈R1×N

(4)

NDVI=[NDVI1,NDVI2,…,NDVIN]∈R1×N

(5)

式中:P為降水?dāng)?shù)據(jù);lat為緯度;lon為經(jīng)度;H為高程;NDVI為歸一化植被指數(shù)。

由于在插值過程中較近的觀測站點對目標點的影響較大,而較遠的觀測站點對目標點的影響較小,本次研究選用距目標點歐氏距離最近的8個地面觀測站點在目標點進行插值。此外,直接計算已知點與未知點之間的距離可能會丟失方向信息,為了解決這個問題,本次研究考慮了點與點之間的相對位置,用相對緯度和相對經(jīng)度作為輸入變量,其計算公式如下:

RLat=(lato-Lat)∈R1×N

(6)

RLon=(lono-Lon)∈R1×N

(7)

式中:lato和lono為目標點的經(jīng)緯度。在添加空間信息后,模型的輸入和輸出可表示為

{[P,RLat,RLon,H,NDVI]→Po}

(8)

模型訓(xùn)練選擇均方誤差MSE作為LSTM的目標函數(shù),并采用學(xué)習(xí)率LR為0.000 1的Adam算法訓(xùn)練模型。神經(jīng)網(wǎng)絡(luò)通常將數(shù)據(jù)集分為訓(xùn)練數(shù)據(jù)集與驗證數(shù)據(jù)集,前者用于訓(xùn)練擬合模型,后者用于調(diào)整模型的超參數(shù)并在模型迭代訓(xùn)練時驗證模型的泛化能力。在本次研究中,將27個地面觀測站點數(shù)據(jù)的時間序列樣本以2∶1的比例分為訓(xùn)練集和驗證集,并根據(jù)訓(xùn)練期和驗證期的模型性能綜合選取最佳的模型參數(shù)。確定最佳模型參數(shù)后將湘江流域748個0.1°格點的高程、歸一化植被指數(shù)和相對經(jīng)緯度以及地面觀測站點的降水?dāng)?shù)據(jù)的時間序列樣本作為輸入,將地面觀測降水?dāng)?shù)據(jù)在湘江流域插值至0.1°。

2.2 動態(tài)貝葉斯模型平均(DBMA)數(shù)據(jù)融合方法

在本次研究中,使用DBMA數(shù)據(jù)融合方法[15],其由不同模型的概率密度函數(shù)(Probability Density Function,PDF)分配不同的權(quán)重來改善降水估計,以獲得與參考降水?dāng)?shù)據(jù)之間良好的適配度。與傳統(tǒng)的為不同模型分配固定權(quán)重的貝葉斯模型平均法相比[16-17],DBMA數(shù)據(jù)融合方法為不同模型分配動態(tài)權(quán)重,該方法的動態(tài)性體現(xiàn)在權(quán)重會隨時間變化,以使其適應(yīng)時間變化引起的天氣模式變化。DBMA數(shù)據(jù)融合結(jié)果的PDF為

(9)

p[y(t)|O(t)]=

(10)

(11)

式中:ai和bi是通過參考降雨數(shù)據(jù)以及不同模型降雨數(shù)據(jù)簡單線性回歸得到,其可以視為簡單的偏差矯正過程。從本質(zhì)上來說,DBMA融合的降水估計是在給定參考降水?dāng)?shù)據(jù)的基礎(chǔ)上,基于不同降水產(chǎn)品的似然值進行加權(quán)平均獲得的。由于直接求解DBMA參數(shù)解析解的復(fù)雜性,通常會采用迭代的方法,本研究選用期望最大化(expectation-maximization,EM)算法求解DBMA參數(shù)[18]。

DBMA融合降水估計y的后驗平均值E(y|O)和方差Var(y|O)為

(12)

(13)

該步驟中訓(xùn)練期是一個滑動窗口,對于時間點t來說,在其之前40 d內(nèi)的各個模型的降雨估計都是其訓(xùn)練數(shù)據(jù)集,隨后利用線性回歸對訓(xùn)練期各個降雨估計進行偏差矯正。最后使用EM算法計算擬合該時間點下各個模型降雨產(chǎn)品的權(quán)重系數(shù),最終實現(xiàn)權(quán)重系數(shù)隨時間動態(tài)變化。

2.3 DHSVM模型

DHSVM模型是由美國太平洋西北國家實驗室(Pacific Northwest National Laboratory,PNNL)和華盛頓大學(xué)(University of Washington,UW)開發(fā)的分布式水文模型[19]。DHSVM利用近地表氣象數(shù)據(jù)(氣溫、風(fēng)速、濕度、降水量以及短波和長波輻射)作為水文氣象輸入來計算能量和水量平衡問題。本次研究使用的模型版本為DHSVM 3.1.2,設(shè)置的網(wǎng)格分辨率為3 000 m。在DHSVM模型中將湘江流域被劃分為6種土壤類型和8種植被類型,其空間分布如圖3所示。

圖3 DHSVM模型中湘江流域Fig.3 Division of Xiangjiang River Basin in DHSVM model

本次研究選用DHSVM模型進行水文模擬,將研究時期分為3個階段:預(yù)熱期(2014年4～12月)、率定期(2015年1月至2016年12月)和驗證期(2017年1～12月)。以徑流模擬的納什效率系數(shù)NSE為目標函數(shù),采用基于ε支配的第二代非支配排序遺傳算法(ε-dominance Non-dominated Sorted Genetic Algorithm Ⅱ,ε-NSGAII)的自動率定模塊進行參數(shù)率定[20]。

2.4 評估指標

本次研究選用相關(guān)系數(shù)(Correlation Coefficient,CC)、均方根誤差 (Root Mean Square Error,RMSE)、相對偏差(Relative Bias,BIAS)3個指標定量評估遙感降水在湘江流域的精度,并采用NSE和BIAS-P評估遙感降水在徑流和洪水事件模擬中的表現(xiàn),具體指標計算公式如下:

(14)

(15)

(16)

(17)

(18)

3 結(jié)果與討論

3.1 融合降水評估

圖4顯示了湘江流域2014年4月1日至2017年12月31日,CMA,IMERG,GAUGE-grid,CMA-E降水產(chǎn)品的平均日降水量的空間分布。根據(jù)CMA的降水空間分布,湘江流域的降水具有較大的空間異質(zhì)性,從東北向西南呈顯著下降趨勢。IMERG-E在流域中西部降水相對較少,其他地區(qū)降水分布相對均勻。GAUGE-grid在流域北部降水較多,可以顯示出較多降水空間分布的細節(jié)。CMA-E則與CMA的空間分布比較吻合,經(jīng)過對IMERG-E和GAUGE-grid的融合,其降低了IMERG-E在流域南部對降水的高估,而在流域中部的降水估計則有所提升,呈現(xiàn)出北多南少,從東北向西南遞減的降水空間分布。

研究統(tǒng)計了以GAUGE為參考降水?dāng)?shù)據(jù),CMA,IMERG-E和CMA-E降水產(chǎn)品在地面站點處的CC,RMSE和BIAS,如表1和圖5所示。根據(jù)評估指標CC,CMA-E降水產(chǎn)品展現(xiàn)出了與地面觀測降水?dāng)?shù)據(jù)較強的相關(guān)性,在大部分站點處其CC要大于IMERG-E和CMA。在部分地面站點,如編號為57 669,57 679,57 682和57 780的站點,CMA-E降水產(chǎn)品的CC高于0.65,表現(xiàn)出了較高的相關(guān)性。在RMSE方面,CMA-E降水產(chǎn)品與IMERG-E,CMA降水產(chǎn)品相比均有所改善,其均值約為11.23 mm,略小于IMERG-E的13.3 mm以及CMA的12.03 mm。通過比較27個地面觀測站點處的BIAS,IMERG-E是3種降水產(chǎn)品中表現(xiàn)最好的,其BIAS絕對值的平均值約為10%,而CMA和CMA-E則為20%左右。CMA對其中22個地面觀測站點的降水存在明顯低估,而CMA-E表現(xiàn)與CMA類似。

注:箱線圖顯示25,50和75百分位,其中正方形表示平均值。圖5 CMA,IMERG-E和CMA-E降水產(chǎn)品在地面觀測站點處的CC,RMSE和BIAS箱線圖Fig.5 CC,RMSE and BIAS of CMA,IMERG-E and CMA-E at ground observation stations

基于以上結(jié)果,本文提出的方法所得到的降水融合產(chǎn)品空間分布與CMA基本一致,降低了IMERG-E在流域西南區(qū)域上的高估值并修正了一些誤報。同時根據(jù)表1結(jié)果,最終所得到的降水與地面觀測結(jié)果有著更強的相關(guān)性,其CC和RMSE相較于IMERG-E有較大的改善。其原因是降水融合過程中采用了LSTM的空間插值方法,該方法融合了研究流域內(nèi)的高程以及NDVI數(shù)據(jù),將27個地面站點處的GAUGE降水?dāng)?shù)據(jù)進行空間插值。但是其BIAS較低,存在對降水的低估現(xiàn)象,這是由于神經(jīng)網(wǎng)絡(luò)方法的局限性,其原理機制存在一定的黑箱效應(yīng),在使用神經(jīng)網(wǎng)絡(luò)方法進行降水融合時也存在類似低估的現(xiàn)象[21]。

總體而言,融合降水產(chǎn)品 CMA-E能夠很好地捕捉湘江流域降水的空間分布,在基于地面站點的評估中,其CC和RMSE相較于IMERG-E有較大的改善,與地面觀測降水有著更高的相關(guān)性。

3.2 基于融合降水的水文模擬

分別以3個降水產(chǎn)品(CMA,IMERG-E,CMA-E)作為DHSVM的降水輸入進行水文模擬,并選用NSE和BIAS-P對模擬性能進行評估,其降雨-徑流過程如圖6所示。由圖6可知,CMA-E徑流模擬的NSE基本與CMA相當(dāng),在率定期分別為0.79和 0.80,在驗證期分別為0.87和0.86,二者均遠高于IMERG-E的0.55與0.75。將CMA和CMA-E作為降水輸入得到的模擬結(jié)果良好,證明了DHSVM模型在湘江流域進行徑流模擬的適用性。

注:CAL表示率定期,VAL表示驗證期。圖6 湘江流域CMA,IMERG-E,CMA-E驅(qū)動下的日尺度DHSVM模型模擬徑流Fig.6 Daily streamflow simulated with DHSVM model drivenby CMA,IMERG-E and CMA-E precipitation in the Xiangjiang River Basin

為了評估遙感降水產(chǎn)品在洪水模擬中的表現(xiàn),本次研究選取了研究時間范圍內(nèi)的11場歷史洪水事件進行評估,選擇的依據(jù)為洪峰流量超過8 600 m3/s(相當(dāng)于97分位數(shù))。表2顯示了將CMA,IMERG-E和CMA-E作為降水輸入的洪水事件模擬的性能。相較于IMERG-E,CMA-E有效改進了洪水事件模擬的性能,率定期的平均NSE從0.26提高到了0.61,驗證期的平均NSE從0.41提高至0.84,對比CMA的0.80也有一定的改進。且IMERG-E所有洪水事件模擬的NSE均大于零,有效改善了IMERG-E無法準確捕捉到的洪水事件模擬。BIAS-P方面,其結(jié)果與NSE相似,CMA-E模擬的洪水事件較IMERG-E在BIAS-P均值上亦有所降低,率定期的均值從21%下降至14%,驗證期的均值由40%下降至21%。

表2 CMA、IMERG-E和CMA-E驅(qū)動下的DHSVM模型對洪水事件模擬的NSE和BIAS-PTab.2 NSE,BIAS-P of flood events simulation using DHSVM driven by CMA,IMERG-E and CMA-E

基于以上結(jié)果可以看出,IMERG-E作為降水輸入時所得到的徑流模擬效果較差,在洪水事件的表現(xiàn)上有著最大的不確定性,無法很好地捕捉洪峰,對多個洪峰存在顯著的高估,低流量時期也存在類似的高估現(xiàn)象。這是因為作為降水輸入的IMERG-E本身存在誤差,而降水估計的誤差會在水文模擬的過程中傳播,導(dǎo)致最后徑流模擬效果較差[22- 23]。相對而言,CMA-E作為降水輸入時所得到的徑流模擬效果有著更好的表現(xiàn),其結(jié)果有著更高的NSE,同時能模擬出更為平滑穩(wěn)定的徑流,更符合實際觀測所得的徑流過程,這證明了本研究所采用的數(shù)據(jù)融合方法能夠得到更高精度的降雨數(shù)據(jù),進而改進其在徑流模擬中的性能。

4 結(jié) 論

本文提出了一種時空動態(tài)數(shù)據(jù)融合方法,并通過評估融合降水?dāng)?shù)據(jù)的精度及在水文模擬中的應(yīng)用,有效證明了所提出的時空動態(tài)數(shù)據(jù)融合方法在改善遙感降水的精度及其在水文模擬中應(yīng)用的適用性。具體結(jié)論如下:

(1) 利用本文所提出的基于LSTM和DBMA的時空動態(tài)數(shù)據(jù)融合方法對湘江流域的地面觀測降水和IMERG-E遙感降水?dāng)?shù)據(jù)進行融合,能夠獲得高精度、高分辨率的CMA-E降水產(chǎn)品。其能夠很好地捕捉湘江流域的降水空間分布,與地面觀測降水具有較高的相關(guān)性。

(2) 基于分布式水文模型DHSVM對CMA-E降水產(chǎn)品在水文模擬中的應(yīng)用進行評估,發(fā)現(xiàn)與CMA,IMERG-E相比,CMA-E能夠更準確地模擬徑流,并且在洪水事件期間,NSE與BISA-P也均有一定程度上的改善。