石羽佳,王忠靜,2,索 瀅

(1. 清華大學水利系,北京 100084；2. 寧夏大學西北土地退化與生態(tài)恢復省部共建國家重點實驗室培育基地,寧夏 銀川 750021)

海河流域位于中國北方,112°E—120°E、35°N—43°N,東臨渤海,南瀕黃河流域,流域面積為32.06萬km2,高程范圍為-10～2 864 m,總地勢是西北高東南低,大致分山地及平原2種地貌類型,是中國三大經(jīng)濟區(qū)之一和重要的工農(nóng)業(yè)基地,戰(zhàn)略地位十分重要。由于受氣候自然條件變化和人類社會經(jīng)濟活動影響,海河流域出現(xiàn)了水資源嚴重不足、土地鹽堿化、水污染等一列生態(tài)環(huán)境問題,流域內(nèi)水資源與水環(huán)境問題突出。降水資源作為流域水資源的最主要部分,對其正確地分析與準確地評價有助于理解流域水資源數(shù)量特征和空間分布特點,因此,亟需明晰海河流域近年來的降水特征及規(guī)律以輔助水資源與水環(huán)境綜合治理。

目前主要有4種估算降水的方法:地基測量觀測、地基雷達遙感、衛(wèi)星遙感和大氣再分析模型[1]。一般認為地基測量觀測值為降水真值,這也是以往研究分析海河流域多年平均年降水量及其空間分布規(guī)律的主要依據(jù)[2- 5]。但當站點空間分布不均、不密時,觀測數(shù)據(jù)只能代表一定范圍內(nèi)的降水值,難以描繪流域整體降水的空間分布。近年來,衛(wèi)星降水觀測因其大范圍面尺度優(yōu)勢成為各類研究與規(guī)劃中的重要參考[6],也為準確評價流域降水資源及其空間分布特征提供了可能,在海河流域也有所評價和應用[7- 10]。然而,衛(wèi)星遙感數(shù)據(jù)并非直接觀測數(shù)據(jù),是由反演得出的,存在一定的不確定性[11],需加以融合校正,以表達流域真實的降水空間分布特征,進而估算更為可靠的降水資源值。

為得到特定研究區(qū)準確的降水空間分布數(shù)據(jù),通常的做法是將衛(wèi)星反演的面降水數(shù)據(jù)與地基觀測的點降水數(shù)據(jù)融合,校正得到特定區(qū)域的降水融合產(chǎn)品,增強衛(wèi)星降水數(shù)據(jù)在應用中的可靠性[12- 15],為區(qū)域水資源管理提供數(shù)據(jù)支撐。然而,以往涉及海河流域的遙感降水成果存在2個主要問題：一是所使用的地基觀測校正站點分布較為稀疏[16- 17]；二是在校正時忽視了與降水密切相關的協(xié)變因子[16- 17],得到的結(jié)果難以再現(xiàn)海河流域降水實際空間分布。因此,有必要引入與降水密切相關的協(xié)變因子,將更密集的地基觀測數(shù)據(jù)與衛(wèi)星遙感相結(jié)合,提供海河流域高質(zhì)量、高分辨率的降水數(shù)據(jù)來更準確地評價海河流域降水資源量及其空間分布。

本文提出多要素融合的降水數(shù)據(jù)校正方法框架,擴展了一致性系數(shù)內(nèi)涵；基于地面站點觀測數(shù)據(jù)、衛(wèi)星降水數(shù)據(jù)、地理時空數(shù)據(jù)等進行多源數(shù)據(jù)融合校正,形成海河流域2001—2019年融合降水數(shù)據(jù)集；利用融合降水數(shù)據(jù)集對海河流域降水資源時空分布進行全面評價。研究成果可為海河流域水資源管理提供更可靠的降水資源數(shù)據(jù)基礎。

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

1.1 地面觀測降水資料

本文將海河流域分為區(qū)域Ⅰ(平原區(qū))、區(qū)域Ⅱ(山區(qū)迎風坡)和區(qū)域Ⅲ(山區(qū)背風坡)3個片區(qū)融合分析。使用的地面觀測降水數(shù)據(jù)來源于中國氣象局國家氣象中心氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/),包括海河流域244個氣象站點2001—2019年逐月實測降水數(shù)據(jù),區(qū)域及站點空間分布如圖1所示。

圖1 海河流域氣象站分布Fig.1 Meteorological stations in Haihe River basin

1.2 衛(wèi)星降水產(chǎn)品

本文使用了2種衛(wèi)星降水產(chǎn)品,分別是Integrated Multi- satellitE Retrievals for the Global Precipitation Measurement Final Run V06(IMERG- F)[18]和Global Satellite Mapping of Precipitation Guage- calibrated Standard V6[19](GSMaP- G)。時間范圍為2001—2019年,時間分辨率為月,空間分辨率為0.1°。IMERG- F來自NASA官網(wǎng)(https:∥pmm.nasa.gov/),GSMaP- G來自JAXA官網(wǎng)(https:∥sharaku.eorc.jaxa.jp/GSMaP/)。

1.3 其他數(shù)據(jù)來源

本文其他數(shù)據(jù)包括表征地面高程信息的ASTER GDEM數(shù)據(jù)(簡稱為DEM)和ERA5- Land數(shù)據(jù)。DEM數(shù)據(jù)來源于地理國情監(jiān)測云平臺(http:∥www.dsac.cn/),ERA5- Land數(shù)據(jù)來源于歐洲中期天氣預報中心(ECMWF，https:∥cds.climate.copernicus.eu/)。ERA5系列產(chǎn)品(ERA5和ERA5- Land)是ECMWF最新和最先進的全球再分析產(chǎn)品[20],已有諸多研究論證了該系列產(chǎn)品在溫度、風速、植被分析等方面與實際觀測具有良好的相關性和較低的誤差[21- 22],可以被作為環(huán)境信息的來源和創(chuàng)建高分辨率氣候數(shù)據(jù)集的基礎[20]。相較ERA5而言,ERA5- Land具有更高的空間分辨率,可以達到0.1°,與本研究選用的衛(wèi)星降水產(chǎn)品一致,因此本文選取ERA5- Land數(shù)據(jù)作為協(xié)變量數(shù)據(jù)來源之一。風場、溫度、氣壓、植被等[22- 25]都與降水有很緊密的反饋與聯(lián)系,本研究選取ERA5- Land數(shù)據(jù)集中的8個要素作為協(xié)變量,以反映10 m風場(10 m 風向u分量、10 m 風向v分量)、溫度與空氣距離飽和的程度(2 m露點溫度、空氣溫度)、氣壓(表面氣壓)、冬季降水(降雪量)和植被縱向分布情況(高大植被葉面積指數(shù)、矮小植被葉面積指數(shù))。

2 研究方法

2.1 融合校正方法

衛(wèi)星降水融合校正方法包括偏差校正法、插值展布法、多元回歸法和機器學習法等[13],效果各不相同[11,12,26- 29]。與傳統(tǒng)的插值或回歸方法相比,機器學習方法具有較強的學習和泛化能力,在建模更復雜的數(shù)據(jù)模式和處理大量數(shù)據(jù)方面具有優(yōu)勢[13,30]。同時,海河流域雨量站分布密度約為1 302 km2/站,研究表明在該雨量站分布密度下更適宜采用機器學習方法進行校正[28]。人工神經(jīng)網(wǎng)絡(Artificial Neural Network,ANN)是典型的機器學習方法,是基于模仿大腦神經(jīng)網(wǎng)絡結(jié)構和功能而建立的一種信息處理系統(tǒng),具有獨特的分布并行處理、非線性映射和自適應學習能力等優(yōu)點,在氣象領域的降水、徑流、蒸散發(fā)等研究中發(fā)揮了重要作用[29- 30],諸多研究表明ANN在數(shù)據(jù)關系擬合、挖掘等方面表現(xiàn)效果較優(yōu)[31- 32]。因此本文以人工神經(jīng)網(wǎng)絡為基礎,以地理要素、氣象要素、植被因子等作為協(xié)變量,構建了海河流域衛(wèi)星降水產(chǎn)品融合校正框架,如圖2所示。

圖2 海河流域衛(wèi)星降水產(chǎn)品融合校正框架Fig.2 Diagram and framework of satellite precipitation product fusion and correction in Haihe River basin

該融合校正框架是以站點降水數(shù)據(jù)為目標值,輸入相應位置柵格中的衛(wèi)星降水值,以經(jīng)緯度、高程及ERA5- Land要素作為協(xié)變量,利用ANN進行關系挖掘。測量值觀測數(shù)據(jù)在空間上隨機分為兩部分,一部分占比70%作為訓練數(shù)據(jù)集,另一部分占比30%作為測試數(shù)據(jù)集。在本研究中,采用了3層結(jié)構(即輸入層、隱含層和輸出層)的多層感知器(multilayer perceptron,MLP),這是神經(jīng)網(wǎng)絡算法中廣泛使用的一種形式,通過擇優(yōu)比較后,將隱含層設定為2層,節(jié)點數(shù)設置為20。

2.2 評價檢驗方法

2.2.1 站點處精度評價檢驗方法

利用氣象站點坐標提取衛(wèi)星降水對應格點數(shù)據(jù),以氣象站點降水觀測數(shù)據(jù)為基準進行站點處精度對比評價。評價指標分別為相關系數(shù)(r)、平均誤差(EM)、相對偏差(SBIA)和均方根誤差(ERMS)。這4個指標是常用的統(tǒng)計參數(shù),此處不再贅述各指標的計算方法。

2.2.2 無站點處精度評價檢驗方法

Xia等[26]提出了一種基于降水- 高程連續(xù)性的一致性假定,用于評價無地表測站柵格的遙感降水產(chǎn)品精度。該假定建立降水與高程因子之間的關系,用一致性系數(shù)(RC)量化指標,并對衛(wèi)星降水在中國8個山區(qū)和柴達木盆地的表現(xiàn)進行了一致性評價[26- 27]。本研究將其影響因子擴展到一般情況,形成廣義一致性系數(shù)。擴展后RC計算如下:

(1)

(2)

(3)

式中:Nsum為衛(wèi)星降水落入掩膜內(nèi)的柵格總數(shù)；N為評價區(qū)柵格總數(shù)；Ci為每個柵格的判定值,i=1,2,,N；Pti為每個柵格的衛(wèi)星降水值，i=1，2，，N；D為降水- 主要因子掩膜,由站點處主要因子數(shù)值與實測降水關系繪制得出,其方法參見文獻[26]。若無地面降水觀測站柵格上的衛(wèi)星降水值落入掩膜內(nèi),則認為其與地面測降水值有一致性,亦即可靠度較高。

3 結(jié)果及分析

3.1 有站點柵格融合校正結(jié)果精度評價

根據(jù)圖2,基于人工神經(jīng)網(wǎng)絡和協(xié)變量要素,將IMERG- F和GSMaP- G迭代融合成多源融合校正數(shù)據(jù)集IG- F(IMERG- F & GSMaP- G Fusion),1 367次后達到最優(yōu),迭代尋優(yōu)收斂過程如圖3。

圖3(a)表示神經(jīng)網(wǎng)絡在迭代過程中的梯度下降過程,在1 000次迭代后下降趨勢顯著減?。粓D3(b)表示神經(jīng)網(wǎng)絡收斂參數(shù)在迭代過程中的結(jié)果,在第1 367次迭代時收斂參數(shù)超出閾值使得迭代停止；圖3(c)表示有效參數(shù)個數(shù)的變化,在400次迭代后趨于穩(wěn)定；圖3(d)表示參數(shù)平方和變化,在800次迭代后趨于穩(wěn)定；圖3(e)表示神經(jīng)網(wǎng)絡的泛化能力檢查過程,可以看出未被包含在訓練過程中的數(shù)據(jù)在該網(wǎng)絡中都可以得到合理的輸出。綜合以上,本研究建立的神經(jīng)網(wǎng)絡模型對此類多柵格多屬性數(shù)據(jù)融合問題的尋優(yōu)搜索是較為有效的。

將多源融合校正數(shù)據(jù)集在有站點柵格上的精度分布結(jié)果繪于圖4。從圖4中可以看出,衛(wèi)星降水產(chǎn)品IMERG- F和GSMaP- G原始數(shù)據(jù)在季節(jié)尺度上與實測值符合性最佳,在年尺度上最差,總體均高估海河流域降水。多源融合校正降水數(shù)據(jù)集IG- F在4個統(tǒng)計指標的表現(xiàn)均有顯著提升,平均r在月、季節(jié)、年尺度上分別為0.94、0.97、0.86,均大于原始產(chǎn)品；平均ERMS在月、季節(jié)、年尺度上分別降低至19.18 mm、33.19 mm、68.63 mm,均小于原產(chǎn)品；平均SBIA在月、季節(jié)、年尺度上均為-0.02%,優(yōu)于IMERG- F的8.6%和GSMaP- G的7.3%；平均EM在月、季節(jié)、年尺度上分別為-0.01 mm、-0.02 mm、-0.08 mm,同樣優(yōu)于原產(chǎn)品。

圖3 模型尋優(yōu)收斂過程Fig.3 Optimization convergence process diagram

圖4 降水數(shù)據(jù)集站點精度分布評價箱線圖Fig.4 Distribution of the precision of precipitation datasets in stations

3.2 無站點柵格融合校正精度結(jié)果評價

對比發(fā)現(xiàn),區(qū)域Ⅰ降水與位置參數(shù)呈現(xiàn)緊密的聯(lián)系,區(qū)域Ⅱ和區(qū)域Ⅲ降水與高程關系最為緊密。因此,式(3) 中降水- 主要因子掩膜在區(qū)域Ⅰ表現(xiàn)為降水- 位置參數(shù)掩膜,區(qū)域Ⅱ和區(qū)域Ⅲ為降水- 高程掩膜。為更直觀地顯示各分區(qū)高程(位置參數(shù))與降水之間的關系,圖5給出了海河流域等高線分布、降水數(shù)據(jù)融合分區(qū)、多年平均年降水等值線和位置參數(shù)等值線分布。在經(jīng)過經(jīng)度值與緯度值不同組合試驗后,各柵格中心經(jīng)度平方值與緯度平方值之和與降水關系最佳,將其開平方形式作為位置參數(shù)，計算過程見式(4) :

(4)

式中:P為柵格位置參數(shù)， °；X為柵格中心經(jīng)度值， °；Y為柵格中心緯度值， °。

圖5 海河流域等高線、降水等值線、位置參數(shù)等值線分布及區(qū)域劃分Fig.5 Contour line,precipitation contour line,location contour line and regional division in Haihe River basin

為了對比衛(wèi)星降水產(chǎn)品在有站點和無站點柵格上的校正效果,本研究把有站點柵格處的降水數(shù)據(jù)集同樣進行了一致性評價,評價結(jié)果繪于圖6。圖6顯示,未經(jīng)融合校正的降水產(chǎn)品在平原區(qū)位置參數(shù)中值區(qū)嚴重高估了降水分布,在山區(qū)迎風坡低估高程高值區(qū)的降水分布,在山區(qū)背風坡高估高程低值區(qū)的降水分布,多源融合數(shù)據(jù)集分布與掩膜更為一致。

對比圖6中有、無地面觀測站點柵格的評價結(jié)果,不難發(fā)現(xiàn),有站點處RC系數(shù)評價結(jié)果均比無站點處的評價結(jié)果要好,表明RC系數(shù)具有客觀性。同時,對比校正前后的評價結(jié)果可知,校正后的降水數(shù)據(jù)集與掩膜分布較為一致,改進后的RC系數(shù)能更好地描述無站點柵格處的一致性。將3個區(qū)域在海河流域的面積占比與其RC值乘積之和作為海河流域的RC值,IMERG- F和GSMaP- G在有站點柵格上的RC分別為55.3%和60.4%,無站點柵格上RC分別為42.7%和48.3%；IG- F在有、無站點柵格上RC分別為80.4%和56.2%。根據(jù)Xia等[26]的研究認為,有站點柵格和無站點柵格的RC值并非完全一致,兩者差異在30%內(nèi)均可接受,因此融合數(shù)據(jù)集精度結(jié)果在理想范圍內(nèi)。對原始衛(wèi)星數(shù)據(jù)進行融合校正后,落入掩膜內(nèi)的數(shù)量均有明顯增加,在有、無站點柵格上的一致性系數(shù)RC值均有提高。以上分析可以說明:改進后的RC系數(shù)通過選取主要因子建立掩膜,能夠反映流域不同主要因子與降水之間的分布關系；生成的多源融合數(shù)據(jù)集在海河流域表現(xiàn)要優(yōu)于原始衛(wèi)星降水數(shù)據(jù)集,研究構建的融合校正模型能夠用于整個研究區(qū)的衛(wèi)星降水校正。

圖6 降水數(shù)據(jù)集一致性評價結(jié)果Fig.6 Evaluation of the consistency of precipitation datasets

3.3 海河流域面降水分析

海河流域氣象站點、原始衛(wèi)星降水和多源融合校正數(shù)據(jù)集2001—2019年平均降水空間分布如圖7所示。分析原始衛(wèi)星產(chǎn)品及多源融合數(shù)據(jù)集在海河流域的空間分布情況不難發(fā)現(xiàn),原始衛(wèi)星降水產(chǎn)品不能很好地識別位于平原中部的低降水量中心,高估平原區(qū)降水,這與圖6(a)、圖6(d)結(jié)果一致；相較于2種原始數(shù)據(jù),IG- F捕捉的降水空間分布與地表站點更為一致,且不存在明顯高估的情況。融合數(shù)據(jù)集顯示海河流域東北、東南、西南、中部偏西降水較多,在西北和中部偏東降水較少,降水空間分布差異較大。

進一步分析圖6可以發(fā)現(xiàn),海河流域各降水分區(qū)的降水存在以下分布規(guī)律:

(1) 平原區(qū)高程變化不明顯,降水空間分布主要受海洋、氣壓帶和風帶影響,與其空間位置參數(shù)存在明顯關系。圖6(a)、圖6(d)顯示位置參數(shù)等值線121°—123°區(qū)間內(nèi)存在1個明顯的降水低值區(qū),與圖4中的降水等值線分布位置一致,說明設定的位置參數(shù)變量能夠較好地描述平原區(qū)降水空間分布。

(2) 山區(qū)迎風坡降水與高程關系密切,高程升高易成云致雨。圖6(b)、圖6(e)顯示海河流域山區(qū)迎風坡高程500 m以下的降水較多,但在高程500～1 000 m間存在1個降水低值區(qū),之后隨高程升高降水再次增加。

(3) 山區(qū)背風坡降水與高程同樣具有緊密聯(lián)系,背風坡因地形阻擋,氣流下沉增溫難以成云致雨,所以降水較少。圖6(c)、圖6(f)顯示背風坡在高程1 200～1 600 m間存在1個降水低值區(qū),主要分布在背風坡山谷內(nèi),顯示出與迎風面不同的降水特征。

圖7 海河流域2001—2019年多年平均降水量空間分布Fig.7 Spatial distribution of annual precipitation in Haihe River basin from 2001 to 2019

3.4 海河流域降水資源評價

計算不同數(shù)據(jù)源下海河流域2001—2019年多年平均面降水量(表1),顯示融合校正數(shù)據(jù)集的多年平均面降水量為515.2 mm,相較于2種原始產(chǎn)品IMERG- F和GSMaP- G分別小4.4%和3.6%,較水資源公報數(shù)據(jù)大1.7%。水資源公報顯示海河流域水資源量占降水資源量的比(產(chǎn)流系數(shù))為17.3%,若將原始產(chǎn)品IMERG- F和GSMaP- G顯示的降水值作為海河流域降水資源量進行計算,產(chǎn)流系數(shù)將分別降低為16.2%和16.4%。綜上,若用衛(wèi)星降水原始數(shù)據(jù)進行流域管理,會高估海河流域降水量,低估產(chǎn)流系數(shù)。

3.5 多源融合變量影響分析

2種衛(wèi)星產(chǎn)品IMERG- F和GSMaP- G在中國大陸地區(qū)降水表達的表現(xiàn)優(yōu)于其他降水產(chǎn)品[33],并且通常認為多衛(wèi)星產(chǎn)品融合可提高結(jié)果的可靠性[34]。本文選用IMERG- F和GSMaP- G生成的融合數(shù)據(jù)集IG- F的表現(xiàn)也證實了這一點。同時,在無協(xié)變量的情況下,衛(wèi)星產(chǎn)品融合存在ERMS誤差增大的現(xiàn)象(增加至48.25 mm),加入?yún)f(xié)變量后,r從0.93上升到0.94,ERMS從48.25 mm下降到19.18 mm，多源融合月尺度試驗方案評價結(jié)果見表2。因此，研究認為加入氣象要素、植被因子等作為協(xié)變量能夠有效地改善融合結(jié)果。

盡管本研究已通過增加協(xié)變量因子提高了數(shù)據(jù)融合的可靠性,但選取的協(xié)變量因子仍然較為有限,忽略了土壤含水量等因素。恰當?shù)娜谌雲(yún)f(xié)變量仍是此類方法今后研究的重點之一。

表1 不同數(shù)據(jù)源及方法計算的海河流域2001—2019年降水資源量

表2 多源融合月尺度試驗方案評價結(jié)果

4 結(jié) 論

為準確評價海河流域降水資源量及其空間分布,本研究提出了海河流域衛(wèi)星降水產(chǎn)品融合校正框架,并利用多源數(shù)據(jù)融合生成了2001—2019年降水數(shù)據(jù)集,評價了海河流域降水資源量。主要結(jié)論如下:

(1) 衛(wèi)星降水產(chǎn)品原始數(shù)據(jù)在海河流域存在一定高估,本研究生成的多源融合降水數(shù)據(jù)集在有無站點柵格上的評價結(jié)果均優(yōu)于原始衛(wèi)星降水產(chǎn)品。

(2) 海河流域降水存在3個較為明顯的空間分區(qū)及其分布規(guī)律。海河流域平原區(qū)降水與空間位置參數(shù)關系明顯,在位置參數(shù)中值區(qū)存在明顯降水低值區(qū)；山區(qū)迎風坡和背風坡降水均與高程關系明顯,在迎風坡高程500～1 000 m、背風坡高程1 200～1 600 m存在降水低值區(qū)。

(3) 多源融合數(shù)據(jù)集IG- F顯示海河流域2001—2019年多年平均面降水量為515.2 mm,合降水資源量1 639.4億m3,較2種原始衛(wèi)星產(chǎn)品分別小4.4%和3.6%,較水資源公報數(shù)據(jù)大1.7%。