融合氣象要素時空特征的深度學(xué)習(xí)水文模型

李 步，田富強(qiáng)，李鈺坤，倪廣恒

(清華大學(xué)水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

20世紀(jì)50年代以來，在氣候變化與人類活動的疊加影響下，洪澇災(zāi)害日益加劇，嚴(yán)重威脅人民群眾的生命財產(chǎn)安全[1-2]。高精度的徑流預(yù)報與模擬是有效應(yīng)對洪澇災(zāi)害和提高水資源管理水平的重要手段[3]。概念性或有物理基礎(chǔ)的水文模型是常用的水文預(yù)報與模擬模型[4]，主要包括新安江和HBV等集總式模型，以及SWAT、MIKE和THREW等分布式模型。但是由于部分水文過程機(jī)制認(rèn)識不足，目前物理水文模型精度仍難以滿足實(shí)際洪水應(yīng)對與高質(zhì)量水資源管理的需求[5-6]。

深度學(xué)習(xí)近年來發(fā)展快速，從最初的多層感知器(即人工神經(jīng)網(wǎng)絡(luò))發(fā)展到結(jié)構(gòu)更加高級的卷積神經(jīng)網(wǎng)絡(luò)[7]和循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)，這些模型廣泛應(yīng)用于水文模擬中[8-9]。然而，由于模型結(jié)構(gòu)限制，這些模型均難以準(zhǔn)確刻畫長時序過程規(guī)律，因此在水文模擬中難以準(zhǔn)確反映諸如融雪、地下水等長時間尺度水文過程[10-11]。長短時記憶神經(jīng)網(wǎng)絡(luò)(LSTM)模型具有獨(dú)特的“門”結(jié)構(gòu)[12]，很好地解決了長時序過程訓(xùn)練梯度爆炸和梯度消失的問題，提高了長時序過程模擬精度。Kratzert等[11]將LSTM模型應(yīng)用在美國241個小流域的徑流模擬中，模擬精度優(yōu)于物理水文模型和RNN。后續(xù)大量研究將LSTM模型成功應(yīng)用在世界不同流域，例如Lees等[13]應(yīng)用于英國669個小流域，Liu等[14]應(yīng)用于中國漢江流域等。然而，目前在基于LSTM模型的徑流模擬研究中，模型輸入數(shù)據(jù)均為流域空間的平均值，損失了流域氣象數(shù)據(jù)的空間信息。以往的研究表明，流域氣象空間異質(zhì)性是影響徑流過程的重要因素[15]，如何將流域氣象要素的空間特征降維輸入至LSTM模型中，以提高LSTM模型徑流模擬精度，值得進(jìn)一步研究。此外，黃河源區(qū)是黃河流域重要的產(chǎn)流區(qū)和水源涵養(yǎng)區(qū)，徑流變化將直接影響黃河中下游的用水安全。黃河源區(qū)屬于青藏高原高寒流域，水分多相態(tài)變化頻繁，高寒下墊面產(chǎn)匯流機(jī)制尚未完全厘清，給黃河源區(qū)高精度徑流模擬預(yù)報與徑流變化研究帶來了挑戰(zhàn)[16]。

本文將利用主成分分析(PCA)和LSTM搭建一種綜合考慮氣象要素時空特性的深度學(xué)習(xí)水文模型(PCA-LSTM)，以黃河源區(qū)為研究區(qū)域，以LSTM模型和THREW模型為比對模型，系統(tǒng)評估PCA-LSTM模型在水文模擬中的適用性和魯棒性。

1 模型原理和搭建

1.1 PCA-LSTM模型

本文提出一種耦合PCA和LSTM的水文模型，如圖1所示，基于PCA提取氣象要素的空間信息，利用LSTM學(xué)習(xí)長時序過程規(guī)律。模型的輸入數(shù)據(jù)為流域降水和氣溫等網(wǎng)格氣象數(shù)據(jù)，輸出為徑流過程。耦合模型內(nèi)部結(jié)構(gòu)與比對模型簡介如下。

1.1.1 主成分分析法

PCA是一種重要的降維方法，通過正交變換將原始變量變換為一組不相關(guān)的變量，實(shí)現(xiàn)利用較少變量來代替原始的較多變量，變換后的較少變量稱之為主成分。主成分變量一般基于累計方差貢獻(xiàn)率確定，當(dāng)方差的累計貢獻(xiàn)率大于或等于一定閾值時(以下簡稱為主成分閾值)，可利用這些成分來刻畫原始變量。PCA具體計算步驟可參考文獻(xiàn)[17]。本文利用PCA提取流域降水和氣溫等空間氣象數(shù)據(jù)的主成分，作為LSTM模型的輸入。主成分閾值設(shè)置為85%、90%和95%，基于徑流模擬性能獲取最優(yōu)主成分閾值。

1.1.2 長短時記憶神經(jīng)網(wǎng)絡(luò)

LSTM模型作為一種RNN的變式[12]，引入了特殊的控制單元，其內(nèi)部結(jié)構(gòu)如圖1(b)所示。其中，細(xì)胞狀態(tài)(ct)和隱藏狀態(tài)(ht)分別捕捉緩慢和快速的演變過程，輸入門(i)、遺忘門(f)和輸出門(o)分別用于控制信息的存儲、丟失和傳遞。這些特殊單元解決了RNN模型在處理長序列數(shù)據(jù)時存在的梯度爆炸或者消失的問題，使得LSTM模型擅于處理長序列過程。LSTM模型的詳細(xì)運(yùn)算過程與在水文模擬中的應(yīng)用可參考相關(guān)文獻(xiàn)[11]。

1.1.3 PCA-LSTM模型

基于以上算法，本文構(gòu)建了耦合PCA和LSTM的水文模型，模型流程主要包括(圖1(a))：

圖1 PCA-LSTM模型結(jié)構(gòu)與LSTM模型內(nèi)部結(jié)構(gòu)Fig.1 PCA-LSTM model architectures and the internals of LSTM cells

(1) 將降水、氣溫等氣象空間數(shù)據(jù)分別重構(gòu)為一維數(shù)據(jù)，并進(jìn)行“最大值-最小值”歸一化等預(yù)處理；

(2) 基于PCA和設(shè)定主成分閾值，分別篩選出降水和氣溫等氣象數(shù)據(jù)的主成分變量，代表氣象要素的空間特征；

(3) 多時刻(1～L)氣象數(shù)據(jù)的主成分變量輸入至LSTM模型，LSTM模型的輸出變量引入全連接層(即線性回歸模型)，獲得待預(yù)測時刻(L)的徑流變量。

1.2 比對模型

為了評估PCA-LSTM模型的適用性，本文搭建了2個比對模型：物理水文模型THREW和深度學(xué)習(xí)模型LSTM。其中LSTM模型的輸入數(shù)據(jù)是流域空間平均的氣象數(shù)據(jù)，其余設(shè)置與PCA-LSTM模型一致。

物理水文模型THREW采用代表性單元流域方法(REW)對流域進(jìn)行空間離散，對每個REW劃分為地上和地下2層，并基于下墊面和土壤特征劃分為8種水文子流域類型，THREW模型具體描述可參考文獻(xiàn)[18-20]。模型已成功應(yīng)用于國內(nèi)外多個流域，包括阿爾卑斯山以及青藏高原多個典型山區(qū)大流域[21-23]。THREW模型參數(shù)率定采用pySOT優(yōu)化算法，該算法應(yīng)用異步并行優(yōu)化框架，充分利用計算資源，采用代理模型的方法尋找最優(yōu)解，減少優(yōu)化模型中的模型運(yùn)行次數(shù)，提升計算效率。率定過程會在運(yùn)行參數(shù)達(dá)到某一限定次數(shù)后停止，本文設(shè)置為3 000。

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

2.1 研究區(qū)域

如圖2所示，研究區(qū)域黃河源區(qū)位于32°N—36°N、95°E—103°E之間，青藏高原東部，流域面積為12.3萬km2，高程范圍為2 656～6 253 m，下墊面以草地、林地為主。黃河源區(qū)屬于典型的高原大陸性氣候，冷熱兩季交替、干濕兩季分明。黃河源區(qū)處于印度洋南亞季風(fēng)和太平洋東亞季風(fēng)控制區(qū)的邊緣交錯帶上，降水同時受暖濕西南季風(fēng)和東亞季風(fēng)的支配，多年平均降水量約為510 mm，呈現(xiàn)由東南向西北遞減的分布特征。流域出口為唐乃亥水文站。黃河源區(qū)是黃河流域的重要產(chǎn)流區(qū)和水源涵養(yǎng)區(qū)，高精度的徑流模擬與預(yù)報可為黃河流域水旱災(zāi)害防御和水資源管理提供參考。

圖2 黃河源區(qū)流域及其水文站位置Fig.2 Study area and the location of the hydrological station

2.2 數(shù)據(jù)

本文構(gòu)建3個模型的輸入數(shù)據(jù)主要包括：

(1) 降水。融合了GLDAS、TRMM 等多種數(shù)據(jù)集及國家氣象站點(diǎn)數(shù)據(jù)的CMFD(China Meteorological Forcing Dataset)數(shù)據(jù)[24-25]，空間分辨率為0.1°。

(2) 氣溫。來源于全球再分析數(shù)據(jù)集ERA5的2m AGL溫度數(shù)據(jù)[26]，空間分辨率為0.1°。

(3) 潛在蒸散發(fā)。來源于全球再分析數(shù)據(jù)集ERA5的潛在蒸散發(fā)數(shù)據(jù)，空間分辨率為0.1°。

(4) 高程。來源于SRTM的90 m分辨率的數(shù)字高程模型數(shù)據(jù)(http:∥srtm.csi.cgiar.org/)。

(5) 土壤。來源于聯(lián)合國糧食和農(nóng)業(yè)組織(FAO)和國際應(yīng)用系統(tǒng)分析研究所構(gòu)建的世界土壤數(shù)據(jù)庫(HWSD)。

(6) 歸一化植被指數(shù)(NDVI)和葉面積指數(shù)(LAI)分別為MODIS的MOD13A2和MOD15A2H數(shù)據(jù)。

其中，LSTM和PCA-LSTM模型輸入均為降水和氣溫數(shù)據(jù)，LSTM模型的輸入為流域空間平均數(shù)據(jù)；PCA-LSTM模型的輸入為流域二維網(wǎng)格數(shù)據(jù)。THREW模型輸入為以上所有數(shù)據(jù)。此外，用于深度學(xué)習(xí)模型訓(xùn)練與物理水文模型率定的目標(biāo)為唐乃亥水文站徑流實(shí)測數(shù)據(jù)(1982—2014年逐日數(shù)據(jù)，來源于西南源區(qū)重大研究計劃集成項目)。氣象與徑流數(shù)據(jù)的時間長度為1982—2014年，時間分辨率為日。

2.3 模型參數(shù)設(shè)置

(1) PCA-LSTM模型。基于相關(guān)研究和前期實(shí)驗(yàn)，設(shè)置降水和氣溫為輸入數(shù)據(jù)，實(shí)測徑流為訓(xùn)練目標(biāo)。其中，1983—2004年為訓(xùn)練期，2006—2009年為驗(yàn)證期，2011—2014年為測試期。為了評估不同主成分閾值對徑流模擬的影響，主成分閾值分別設(shè)置為85%、90%和95%。PCA-LSTM模型中LSTM的主要參數(shù)設(shè)置為：LSTM層數(shù)為2，隱藏層數(shù)目為32，時間序列長度為120 d。模型損失函數(shù)為均方根誤差，優(yōu)化器采用Adaptive Moment Estimation (Adam)。模型批尺寸為32，迭代輪次為200。LSTM輸入數(shù)據(jù)采用“最大值-最小值”方法進(jìn)行歸一化。

(2) LSTM模型。模型輸入為流域空間平均的降水和氣溫數(shù)據(jù)，其余設(shè)置與PCA-LSTM模型的LSTM部分保持一致。

(3) THREW模型。黃河源區(qū)劃分83個REWs，率定參數(shù)包括上、下層土壤飽和導(dǎo)水率平均值等17個參數(shù)[19]，率定期設(shè)置為1983—2004年，測試期設(shè)置為2011—2014年，預(yù)熱期為2 a。

(4) LSTM和PCA-LSTM模型對輸入數(shù)據(jù)不確定性的魯棒性評估。為模型中LSTM的輸入數(shù)據(jù)添加不同高斯噪音，評估模型性能隨著噪音增加的變化趨勢，若模型性能隨著噪音平滑增大平滑下降，表明模型未出現(xiàn)數(shù)據(jù)過擬合的現(xiàn)象，模型具有很好的魯棒性。參考Kratzert等[27]的相關(guān)研究，設(shè)置10組符合正態(tài)分布的噪音，均值為0,標(biāo)準(zhǔn)差分別為0.02、0.04、…、0.18、0.20。對于每組噪音，生成50條隨機(jī)噪音向量，以此降低噪音自身隨機(jī)性對模型性能的影響。

以上所有深度學(xué)習(xí)程序均由Python的Pytorch庫編譯實(shí)現(xiàn)。

2.4 評價指標(biāo)

本文選取納什效率系數(shù)[28](ENS)及其3個分解變量[29](相關(guān)系數(shù)(r)，標(biāo)準(zhǔn)差偏差(α)和均值偏差(β))用于系統(tǒng)評估模型性能。r用于評估模擬值與實(shí)測值的趨勢一致性，α與β分別評估模擬值標(biāo)準(zhǔn)差與平均值的偏差。

(1)

(2)

(3)

(4)

式中：Qsim，l(Qobs，l)分別為時刻l的模擬(觀測)徑流;L為模擬期長度;σsim(σobs)和Qsim(Qobs)分別為模擬(觀測)的徑流的標(biāo)準(zhǔn)值和均值。此外，基于觀測值的超越概率設(shè)置洪峰偏差(Bp)、中值偏差(Bm)和低值偏差(Bl)，評估模型在不同流量條件下的性能表現(xiàn)[30]。

(5)

式中：i分別代表p(洪峰，超越概率為0～0.02)、m(中值，超越概率為0.3～0.7)和l(低值，超越概率為0.7～1);Li為模擬期中以上不同超越概率對應(yīng)的數(shù)據(jù)長度。

3 結(jié)果與討論

3.1 PCA-LSTM模型適用性評估

由于氣溫數(shù)據(jù)的空間相關(guān)性較高，第一主成分變量的方差累計貢獻(xiàn)率達(dá)96%，本文設(shè)置的3個主成分閾值對應(yīng)的氣溫主成分?jǐn)?shù)量均為1。在主成分閾值設(shè)置為85%、90%和95%時，降水主成分?jǐn)?shù)量分別為9、14和29。圖3展示了降水和氣溫的第一主成分變量與空間平均值的相關(guān)關(guān)系，二者在降水與氣溫的相關(guān)系數(shù)r2>0.99，表明降水與氣溫的第一主成分變量包含信息與其流域空間平均值基本一致。

圖3 黃河源區(qū)降水與氣溫的第一主成分變量與流域空間平均值比較Fig.3 Comparisons between precipitation and temperature data of first principal component and basin spatially average

圖4展示了PCA-LSTM模型在主成分閾值分別設(shè)置為85%、90%和95%時唐乃亥水文站徑流的模擬與實(shí)測徑流結(jié)果，評價指標(biāo)結(jié)果如表1所示。3種不同主成分閾值的PCA-LSTM模型模擬ENS分別為0.92、0.92和0.93，r分別為0.96、0.96與0.97，模擬徑流與實(shí)測徑流具有很好的一致性；此外，均值偏差結(jié)果說明3個主成分閾值對應(yīng)的模擬徑流均存在略微低估，低估小于5%；其中，洪峰、中值與低值分別出現(xiàn)低估、低估與高估，偏差絕對值均低于20%。說明模型偏差較小，能夠準(zhǔn)確模擬徑流量，且不同主成分閾值的模擬結(jié)果之間差異較小。綜合以上指標(biāo)表明本文提出的PCA-LSTM模型能夠準(zhǔn)確捕捉水文氣象要素的時空特征，在黃河源區(qū)徑流模擬表現(xiàn)優(yōu)異。此外，主成分閾值分別為85%、90%和95%時，ENS均超過0.90，模擬表現(xiàn)無明顯差異，說明主成分閾值為85%時PCA方法能夠提取滿足黃河源區(qū)降水—徑流過程模擬的氣象空間特征，進(jìn)一步增加主成分閾值對徑流模擬影響較小。

圖4 不同主成分閾值下PCA-LSTM模型模擬和實(shí)測徑流過程線Fig.4 Performance of streamflow simulations of PCA-LSTM models with various thresholds in the test period

表1 不同主成分閾值下PCA-LSTM模型評價指標(biāo)結(jié)果Table 1Evaluation metrics results of PCA-LSTM models with different PCA thresholds

3.2 不同模型結(jié)果的對比

圖5和表2展示了PCA-LSTM模型(主成分閾值為85%)與比對模型THREW和LSTM的模擬徑流過程和評價指標(biāo)結(jié)果。THREW、LSTM與PCA-LSTM模型的ENS分別為0.75、0.88和0.92，r分別為0.91、0.95和0.96，表明3個模型均能較好模擬徑流趨勢，PCA-LSTM模型的表現(xiàn)最為優(yōu)異，LSTM模型次之，THREW模型最差。此外，PCA-LSTM模型徑流量偏差(-2%)遠(yuǎn)低于LSTM模型(-11%)與THREW模型(-23%)；比較不同階段的徑流偏差可以發(fā)現(xiàn)，PCA-LSTM與LSTM模型的低值和中值偏差遠(yuǎn)低于THREW模型，PCA-LSTM模型峰值與低值偏差在3個模型中最低，但3個模型在徑流峰值模擬中均存在高于10%的低估。

表2 不同模型評價指標(biāo)結(jié)果Table 2Evaluation metrics results of various models

圖5 不同模型模擬和實(shí)測徑流過程線Fig.5 Performance of streamflow simulations of various models in the test period

物理水文模型由于對下墊面特征、含水層結(jié)構(gòu)與部分水文過程刻畫不足，徑流過程模擬總體存在一定誤差。徑流低水過程階段較小的絕對誤差會產(chǎn)生較大的相對誤差，使得THREW模型對于低水過程模擬的相對偏差顯著高于其他階段。而深度學(xué)習(xí)模型是直接基于歷史氣象與徑流觀測數(shù)據(jù)訓(xùn)練擬合降水—徑流過程[6]，由于強(qiáng)大的數(shù)據(jù)挖掘和擬合能力，LSTM和PCA-LSTM模型徑流模擬精度高于THREW模型。深度學(xué)習(xí)模型的優(yōu)異性能證明了基于現(xiàn)有水文氣象數(shù)據(jù)可以挖掘更為準(zhǔn)確的降水—徑流關(guān)系，解析深度學(xué)習(xí)水文模型內(nèi)部結(jié)構(gòu)，揭示隱藏的降水—徑流過程機(jī)理規(guī)律，將有助于推動對水文過程的認(rèn)識與提高物理水文模型的模擬精度。另一方面，對于深度學(xué)習(xí)模型，研究區(qū)枯水期的氣象與徑流年際差異較小，降水—徑流過程關(guān)系較為單一，而在豐水期降水與徑流的年際差異顯著，降水—徑流過程關(guān)系更為復(fù)雜，使得豐水期降水—徑流過程訓(xùn)練難度更大。此外，研究區(qū)下墊面空間異質(zhì)性明顯，產(chǎn)匯流機(jī)制差異大，因此降水的空間分布對徑流高水過程影響顯著。LSTM模型只利用黃河源區(qū)降水與氣溫的流域平均值作為輸入，忽略了降水的空間特征，導(dǎo)致徑流峰值模擬存在較大偏差。本文提出的PCA-LSTM模型考慮了降水與氣溫的空間特征，與LSTM模型相比，ENS、徑流偏差與峰值偏差分別提高了5%、9%和7%，表明氣象要素空間特征可以提升深度學(xué)習(xí)水文模擬表現(xiàn)。

3.3 模型魯棒性評估

圖6展示了LSTM和PCA-LSTM模型在添加不同標(biāo)準(zhǔn)差的高斯噪音后的模擬結(jié)果。隨著高斯噪音標(biāo)準(zhǔn)差的增加，ENS降低，變化范圍增大，表明模擬表現(xiàn)逐漸降低，噪音自身的隨機(jī)性對模擬結(jié)果的影響增加。但是模擬表現(xiàn)隨著噪音光滑增加呈現(xiàn)光滑下降而不是驟降趨勢，表明PCA-LSTM模型未對作為LSTM模型輸入的主成分變量過擬合，具有良好的魯棒性。

圖6 基于高斯噪音的LSTM與PCA-LSTM模型魯棒性評估結(jié)果Fig.6 Performance of LSTM and PCA-LSTM models with increasing Gaussian noises added to LSTM inputs

4 結(jié) 論

(1) 本文提出了融合氣象要素時空特性的深度學(xué)習(xí)水文模型PCA-LSTM，利用主成分分析提取了氣象要素的空間特征，基于長短時記憶神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)了長時序過程規(guī)律，提高了深度學(xué)習(xí)水文模型對氣象要素空間異質(zhì)性的刻畫能力。

(2) 利用比對模型THREW和LSTM在黃河源區(qū)對PCA-LSTM模型進(jìn)行系統(tǒng)評估，結(jié)果表明PCA-LSTM模型徑流模擬納什效率系數(shù)達(dá)到0.92，較比對模型THREW和LSTM分別提高了23%和5%，證明PCA-LSTM模型具有良好的適用性與魯棒性，在深度學(xué)習(xí)模型中引入氣象要素空間特征能夠提升徑流模擬精度。

本文主旨是在LSTM水文模型中考慮氣象要素空間信息，但氣象要素空間特征提取方法較多，后續(xù)將利用圖像識別領(lǐng)域神經(jīng)網(wǎng)絡(luò)開展進(jìn)一步研究。此外，物理水文模型由于具有明確的物理機(jī)制，仍然是實(shí)際水文預(yù)報中最常用的工具，如何深入耦合物理水文模型與深度學(xué)習(xí)模型，同時發(fā)揮二者優(yōu)勢，值得未來深入研究。