基于LSTM的青藏高原凍土區(qū)典型小流域徑流模擬及預(yù)測

黃克威，王根緒，宋春林，俞祁浩

（1.中國科學(xué)院、水利部成都山地災(zāi)害與環(huán)境研究所，四川成都610041；2.中國科學(xué)院大學(xué)，北京100049；3.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川成都610065；4.中國科學(xué)院西北生態(tài)環(huán)境資源研究院凍土工程國家重點(diǎn)實(shí)驗(yàn)室，甘肅蘭州730000）

0 引言

凍土作為冰凍圈的主要組成要素之一，由于其對氣候變化的高度敏感性而受到廣泛關(guān)注［1-3］，其不同于非凍土區(qū)的水循環(huán)和三水轉(zhuǎn)化關(guān)系決定了多年凍土區(qū)具有完全不同的徑流形成過程、機(jī)制與季節(jié)動(dòng)態(tài)［4］。凍土由于其低滲透性作為隔水層而形成了凍結(jié)層上地下水，同時(shí)大部分的水文過程都被限制在凍土活動(dòng)層中，凍融過程中伴隨的水熱變化直接影響土壤水再分配和土壤水儲(chǔ)量［5-6］，且土壤的儲(chǔ)水能力及導(dǎo)水系數(shù)隨凍融過程而變化［7］。在多年凍土區(qū)，春季活動(dòng)層開始融化，初期以蓄滿產(chǎn)流為主，隨著氣溫升高，活動(dòng)層融化深度增加，壤中流開始出現(xiàn)，當(dāng)夏季融化深度較大時(shí)，土壤下滲能力增加，此時(shí)以超滲產(chǎn)流為主導(dǎo)作用，當(dāng)秋季隨著氣溫降低，凍結(jié)層上地下水迅速減小，同時(shí)隨著凍結(jié)鋒面的上升，產(chǎn)流以蓄滿產(chǎn)流為主［4，8-11］。因此，凍土區(qū)多種產(chǎn)流方式并存，且受溫度因素控制［8-9，11］，在凍土水文模擬過程中，應(yīng)重視溫度因素的影響。

然而，傳統(tǒng)的水文模型因未考慮溫度因素引起活動(dòng)層的變化而不適用于凍土區(qū)，因此學(xué)者們相繼提出了適用于凍土地區(qū)的水文模型。一方面，學(xué)者們對傳統(tǒng)水文模型加以改進(jìn)使其適用于凍土地區(qū)，如Kang等［12］、Lindstrom等［13］分別采用概念水庫、增加土壤凍結(jié)模型的方式改進(jìn)20世紀(jì)70年代瑞典氣象水文局開發(fā)的集總式概念性模型HBV，并將其應(yīng)用于我國黑河山區(qū)流域、瑞典北部小流域；關(guān)志成等［14］考慮凍土形成是累積負(fù)氣溫的函數(shù)，改進(jìn)了新安江模型，建成了具有物理基礎(chǔ)的概念性寒區(qū)流域水文模型。近年來，周劍等［15］利用BP神經(jīng)網(wǎng)絡(luò)識(shí)別水文單元凍土面積改進(jìn)了USGS開發(fā)的半分布式模型PRMS；Qi等［16-18］開發(fā)一個(gè)具有物理基礎(chǔ)的溫度模塊用以改進(jìn)USDA開發(fā)的分布式水文模型SWAT；李明亮等［19］引入凍結(jié)土壤導(dǎo)水率隨氣溫呈指數(shù)變化的概化模型改進(jìn)了楊大文開發(fā)的大尺度分布式水文模型GBHM。另一方面，學(xué)者們也發(fā)展了新的包含凍土模塊的水文模型，如VIC模型［20］、GEOtop模 型［21］、CRHM模型［22］、GBEHM模 型［23］、CBHM模型［24-25］和WEB-DHM-SF模型［26］，充分考慮了因溫度變化引起的凍土活動(dòng)層凍融循環(huán)，從而能較好地模擬凍土水文過程。然而，由于考慮了凍土因素，上述模型大多需要更多的、更精細(xì)的輸入條件［27-30］，但多年凍土區(qū)往往分布在高緯度、高海拔地區(qū)，氣候條件惡劣、人煙稀少，觀測站點(diǎn)分布有限、極不均勻且觀測要素有限，因此上述凍土水文模型在缺少相應(yīng)實(shí)測數(shù)據(jù)的小流域區(qū)域應(yīng)用受到了很大的限制。

人工神經(jīng)網(wǎng)絡(luò)作為一種黑箱模型，不需要流域相關(guān)物理參數(shù)即可進(jìn)行模擬，大大減少了模型的參數(shù)，較常規(guī)凍土水文模型更加適用于無資料、缺資料地區(qū)。相較于其他經(jīng)典線性黑箱模型如自回歸（AR）、滑動(dòng)平均（MA）、自回歸滑動(dòng)平均（ARMA）、線性回歸（LR）、多元線性回歸（MLR）等模型，人工神經(jīng)網(wǎng)絡(luò)（ANN）由于具有非線性、自適應(yīng)學(xué)習(xí)和容錯(cuò)性等特點(diǎn)而被廣泛地應(yīng)用于水文領(lǐng)域［28，31-32］，包括凍土水文［33-35］。然而，普通的前饋人工神經(jīng)網(wǎng)絡(luò)，如BP、RBF模型在處理時(shí)間序列時(shí)會(huì)丟失時(shí)序信息，不能記憶之前的輸入信息［28，36］。循環(huán)神經(jīng)網(wǎng)絡(luò)（RNN）是一種特殊的神經(jīng)網(wǎng)絡(luò)，其通過內(nèi)部自循環(huán)神經(jīng)元存儲(chǔ)和提取時(shí)間序列中的動(dòng)態(tài)信息，既實(shí)現(xiàn)了序列信息的記憶，又可將之前的信息用于之后時(shí)刻的計(jì)算中，非常適用于處理水文中的時(shí)間序列數(shù)據(jù)［36］。但是，經(jīng)典的RNN模型存在梯度消失、爆炸的問題［37］，難以真正有效地利用長距離的時(shí)序信息。長短期記憶人工神經(jīng)網(wǎng)絡(luò)（long short-term memory，LSTM）的出現(xiàn)解決了這一問題，其與標(biāo)準(zhǔn)的RNN模型結(jié)構(gòu)基本相同，但擁有更加細(xì)化的內(nèi)部處理單元，能真正有效地利用長距離的時(shí)序信息［28，36，38］。LSTM包含特殊的細(xì)胞狀態(tài)和門結(jié)構(gòu)，可動(dòng)態(tài)地控制時(shí)間序列信息的流動(dòng)和保存，能夠捕獲水文變量間長時(shí)間的依賴性和水流路徑連通性的變化［39］，從而提高了徑流的模擬精度。黨池恒等［36］將LSTM應(yīng)用于受季節(jié)性降雪影響的岷江源頭區(qū)，并與RNN、BP對比，結(jié)果表明LSTM實(shí)現(xiàn)了流域狀態(tài)的長期存儲(chǔ)和更新，徑流模擬結(jié)果明顯優(yōu)于RNN、BP模型。Kratzert等［28］將LSTM應(yīng)用于美國CAMELS數(shù)據(jù)集中的241個(gè)流域，結(jié)果表明在受降雪影響的流域中，LSTM相較于RNN能較好地模擬融雪徑流，且模擬結(jié)果與具有物理基礎(chǔ)的SACSMA+Snow-17模型模擬結(jié)果相當(dāng)。與受降雪影響的流域相比，凍土區(qū)徑流受降水和氣溫的雙重影響，同時(shí)活動(dòng)層的凍融循環(huán)增加了凍土水文模擬的復(fù)雜性，LSTM特殊的結(jié)構(gòu)能夠有助于提高模擬精度，但目前關(guān)于這方面的研究還很少。

青藏高原，素有“第三極”之稱，海拔普遍在4 000 m以上，面積約250萬km2，是中低緯度地區(qū)多年凍土最大分布區(qū)，多年凍土面積約為106萬km2，約占42.4%［2，40-41］。青藏高原對氣候變化敏感，而氣溫升高將導(dǎo)致多年凍土退化［42-44］，從根本上改變凍土區(qū)的水文地質(zhì)條件，導(dǎo)致地下水動(dòng)態(tài)產(chǎn)生顯著的變化，從而改變凍土流域的徑流過程［40，45］。目前，關(guān)于氣候變化下的青藏高原上的大型流域的徑流變化研究較多，如三江源、黑河流域等［34，46-49］，但青藏高原凍土小流域由于觀測條件的限制，具有物理基礎(chǔ)的凍土水文模型難以應(yīng)用，而一般的黑箱模型也難以精確模擬凍土流域徑流，因此相關(guān)研究還較少。然而凍土退化對于高凍土覆蓋率區(qū)域（＞60%）的徑流影響較大，且隨著凍土覆蓋率升高而增加［40，50-51］。長江源凍土覆蓋率達(dá)76%，氣溫升高對徑流過程影響顯著［34］，因此位于長江源的凍土覆蓋率100%的小流域，氣溫升高對徑流的影響不可忽視。本文選取位于青藏高原腹地的長江源區(qū)風(fēng)火山流域?yàn)檠芯繉ο螅荚谝郧嗖馗咴陷^易獲取的降水、氣溫作為模型輸入，基于LSTM及凍土流域產(chǎn)流機(jī)制，建立一個(gè)適用于凍土小流域的水文模型，并探究研究氣候變化下，風(fēng)火山流域的徑流變化。同時(shí)，為了驗(yàn)證模型的可靠性，將模型應(yīng)用于沱沱河流域。

1 研究區(qū)概況

本文選取位于青藏高原腹地、長江源多年凍土區(qū)的風(fēng)火山流域（92°50′～93°03′E、34°40′～34°47′N）作為研究區(qū)域，風(fēng)火山流域是北麓河的二級支流、通天河的三級支流，干流河長17.07 km，河道平均比降18.92‰，流域面積117 km2，如圖1所示。流域內(nèi)東南高、西北低，植被主要為高寒草甸和高寒沼澤草甸［52］，主要分布在海拔5 000 m以下，土壤層厚度約為30～80 cm，主要為壤土和砂質(zhì)壤土［53］。流域內(nèi)多年凍土覆蓋率100%，屬于長江河源高平原連續(xù)多年凍土區(qū)丘陵亞區(qū)，多年凍土厚度60～120 m，凍土活動(dòng)層厚度1.3～2.5 m［54-58］。風(fēng)火山流域寒冷、干燥，屬于典型的高原內(nèi)陸氣候，近十年年平均氣溫為-5.2℃，年平均降水量為328.9 mm，降水主要集中在6—9月，占總降水量的80%以上；11月至次年4月降水小于5%［9，53-54］。流域內(nèi)徑流受積雪、凍土的影響，年內(nèi)可劃分為春汛期（5月初至6月下旬）、夏季退水期（6月下旬至7月底）、夏汛期（8月初至9月初）、秋季退水期（9月初至10月中旬）、冬季凍結(jié)期（11月至次年4月初）共5個(gè)時(shí)段［59］。模型使用的風(fēng)火山流域日降水、氣溫及徑流數(shù)據(jù)來源于流域內(nèi)氣象站、水文站（如圖1所示），時(shí)間范圍為2017—2019年。

圖1 研究區(qū)域位置Fig.1 Location of the study area

風(fēng)火山流域鄰近的水文站為沱沱河水文站（34°13′N、92°27′E），是國家重要水文站，屬沱沱河流域（33°22′～35°12′N、89°48′～92°54′E），流域內(nèi)有與水文站相鄰的沱沱河氣象站（34°13′N、92°26′E）。沱沱河流域?yàn)殚L江正源，流域面積15 924 km2，位于青藏高原腹地，氣候寒冷、干燥，多年平均氣溫為-4.2℃，多年平均降水量為283 mm，多年平均流量為26.2 m3·s-1［60-61］，流域內(nèi)多年凍土覆蓋率極高（如圖1所示）。本研究采用沱沱河流域1990—2019年逐日降水、氣溫及徑流數(shù)據(jù)用以驗(yàn)證基于LSTM的凍土水文模型，沱沱河氣象站降水、氣溫?cái)?shù)據(jù)來源于國家氣象科學(xué)數(shù)據(jù)中心，沱沱河水文站徑流數(shù)據(jù)來源于青海省水文水資源勘測局。

2 研究方法

2.1 基于LSTM的凍土水文模型

不同于一般的RNN，LSTM增加了細(xì)胞狀態(tài)（Ct）這個(gè)關(guān)鍵變量來存儲(chǔ)長期記憶信息，并由遺忘門（Ft）、輸入門（INt）和輸出門（Ot）這三個(gè)門結(jié)構(gòu)來調(diào)整細(xì)胞狀態(tài)［37，62］。其中，遺忘門（Ft）決定了t時(shí)刻細(xì)胞狀態(tài)需要移除的t-1時(shí)刻細(xì)胞狀態(tài)的信息，輸入門（INt）決定了t時(shí)刻細(xì)胞狀態(tài)需要存儲(chǔ)的新信息，輸出門（Ot）決定了t時(shí)刻細(xì)胞狀態(tài)需要輸出的信息，而t時(shí)刻的細(xì)胞狀態(tài)則記錄了t時(shí)刻的輸入、門結(jié)構(gòu)信息及t-1時(shí)刻隱藏層狀態(tài)、t-1時(shí)刻細(xì)胞狀態(tài)。本文基于LSTM建立了適用于多年凍土區(qū)徑流模擬的水文模型，模型結(jié)構(gòu)如圖2所示，圖中LSTM神經(jīng)單元中的Ⅰ、Ⅱ、Ⅲ分別表示遺忘門、輸入門和輸出門。模型具體計(jì)算過程如下：

圖2 基于LSTM的凍土水文模型結(jié)構(gòu)Fig.2 Structure of permafrost hydrology model based on LSTM

①將t時(shí)刻的降水Pt、氣溫Tt標(biāo)準(zhǔn)化后作為模型輸入。

②通過遺忘門（Ft）移除t-1時(shí)刻細(xì)胞狀態(tài)的信息。

③通過輸入門（INt）確定用以更新細(xì)胞狀態(tài)的信息Ct。

④通過輸出門（Ot）計(jì)算細(xì)胞狀態(tài)的輸出及生成的隱藏層狀態(tài)變量Ht。

⑤將模型輸出去標(biāo)準(zhǔn)化后，得到模擬的凍土區(qū)徑流RMt。

式中：Pt為t時(shí)刻的降水（mm）；Tt為t時(shí)刻的氣溫（℃）；It為模型在t時(shí)刻的輸入（包含t時(shí)刻的降水和氣溫）；Rt為模型模擬的t時(shí)刻的徑流（m3·s-1）；WF、bF分別為遺忘門（Ft）的權(quán)重矩陣和偏置；Ht-1、Ht分別為t-1和t時(shí)刻的隱藏層狀態(tài)；WIN、bIN分別為輸入門（INt）的權(quán)重矩陣和偏置；C~t為細(xì)胞狀態(tài)（Ct）的更新向量；WC、bC分別為細(xì)胞狀態(tài)更新向量C~t的權(quán)重矩陣和偏置；WO、bO分別為輸出門（Ot）的權(quán)重矩陣和偏置；σ為sigmoid激活函數(shù)tanh為sigmoid函數(shù)的變形函數(shù)tanh（x）激活函數(shù)，為 標(biāo) 準(zhǔn) 化 函 數(shù)，fn(x)=為去標(biāo)準(zhǔn)化函數(shù)，zi=fd(yi)=yizsd+zm，其中xi、zi為x、z數(shù)組的任意值，xm、zm為x、z數(shù)組的均值，xsd、zsd為x、z數(shù)組的標(biāo)準(zhǔn)差，yi為標(biāo)準(zhǔn)化后的zi。

2.2 模型評價(jià)指標(biāo)

采用決定系數(shù)（R2）、納什效率系數(shù)（NSE）及均方根誤差（RMSE）指標(biāo)來評價(jià)模型的模擬效果，具體計(jì)算公式為

式中：ROt為t時(shí)刻實(shí)測徑流（m3·s-1）；RMt為t時(shí)刻模擬徑流（m3·s-1）為實(shí)測徑流的均值為模擬徑流的均值（m3·s-1）；n為自模擬開始的第n時(shí)刻。

3 結(jié)果與分析

3.1 模型訓(xùn)練與驗(yàn)證

本研究以風(fēng)火山流域2017—2018年作為模型訓(xùn)練期、2019年作為模型驗(yàn)證期。模型訓(xùn)練期，除了需要率定模型細(xì)胞狀態(tài)（Ct）、遺忘門（Ft）、輸入門（INt）和輸出門（Ot）相關(guān)參數(shù)，還需要率定LSTM隱藏層層數(shù)（numLayer）、單個(gè)隱藏層所含神經(jīng)元的數(shù)量（numHiddenUnit）、模型進(jìn)行完整訓(xùn)練的最大次數(shù)（MaxEpochs）、學(xué)習(xí)率下降周期（LearnRateDrop-Period）和學(xué)習(xí)率下降因子（LearnRateDropFactor）等超參數(shù)，模型主要超參數(shù)及其取值如表1所示。如圖3（a）所示，模型訓(xùn)練期內(nèi)模擬結(jié)果較好，R2、NSE達(dá)0.93，RMSE為0.63 m3·s-1，且模型在春汛期、夏季退水期、夏汛期、秋季退水期、冬季凍結(jié)期都能較好地模擬徑流［圖3（b）］。將訓(xùn)練好的模型應(yīng)用于2019年徑流模擬，用以驗(yàn)證模型的可靠性。如圖4所示，雖然2019年降水年內(nèi)分配不同于2017年、2018年，最大降水最主要集中在7月初和8月底，但模型仍能較好的模擬年內(nèi)各時(shí)段的徑流，R2、NSE分別為0.81、0.77，RMSE為0.69 m3·s-1。模型訓(xùn)練期、驗(yàn)證期在風(fēng)火山流域都能較好地模擬徑流，因此模型能夠用于凍土水文模擬中。

圖3 模型訓(xùn)練期模擬結(jié)果Fig.3 Simulation results during model training：comparison between simulated runoff and measured runoff during model training（a），and simulation results of runoff process during model training（b）

圖4 模型驗(yàn)證期模擬結(jié)果Fig.4 Simulation results during model validation：comparison between simulated runoff and measured runoff during model validation（a），and simulation results of runoff process during model validation（b）

表1 模型主要超參數(shù)Table 1 Main hyperparameters of the model

為了進(jìn)一步驗(yàn)證LSTM模型在凍土流域的可靠性，將模型應(yīng)用于同屬于長江源的沱沱河流域，以1990—2009年作為模型訓(xùn)練期、2010—2019年作為驗(yàn)證期。模型模擬結(jié)果如圖5所示，雖然相較于風(fēng)火山流域，沱沱河流域模擬過程中，汛期峰值的模擬結(jié)果稍差，但是考慮沱沱河流域面積較大，汛期降水空間分布不均勻，本次模擬中僅采用了沱沱河雨量站的降水、氣溫資料可能會(huì)導(dǎo)致一定程度的誤差，因此結(jié)果是較為合理的。總體上，模型訓(xùn)練 期R2、NSE均 為0.73，驗(yàn) 證 期R2、NSE分 別 為0.66、0.64，與寒區(qū)水文模型CRHM、WEB-DHMSF在長江源區(qū)模擬結(jié)果相當(dāng)［26，63］，且模型結(jié)果在豐水年、枯水年、暖年、冷年中均較好，因此模型可靠。

圖5 沱沱河流域模擬結(jié)果Fig.5 Simulation results of Tuotuohe watershed

3.2 不同氣候變化情景下徑流的模擬

近10年來，青藏高原腹地氣溫呈升高趨勢，變化率約為0.5℃·（10a）-1，同時(shí)降水也成呈上升趨勢［34，40］。為了分析不同氣候變化條件下風(fēng)火山流域徑流的變化規(guī)律，依據(jù)實(shí)際的氣候變化情況，本研究以2019年為基準(zhǔn)，設(shè)置了10種可能的氣候變化情景：①降水不變，氣溫分別增加0.5℃、1.0℃、1.5℃、2.0℃；②降水增加10%，同時(shí)氣溫分別增加0℃、1.0℃、2.0℃；③降水增加20%，同時(shí)氣溫分別增加0℃、1.0℃、2.0℃。將10種氣候變化情況下徑流模擬結(jié)果按月統(tǒng)計(jì)，并與基準(zhǔn)年相比較，得到不同氣候變化情景下風(fēng)火山流域相對于基準(zhǔn)年的各月平均徑流增加幅度。由表2可知，降水增加將導(dǎo)致全年徑流增加，但年內(nèi)各月增加幅度不同，其中8—9月增幅最大，6—7月次之，4—5月及10月增幅較小，11—12月徑流未變化，總體上降水每增加10%，年徑流約增加12%；氣溫增加將導(dǎo)致全年除8月外各月徑流的不同幅度的增加，其中11—12月最大，4月、7月次之，5—6月、9—10月較小，而8月徑流隨氣溫的增加而減少，總體上氣溫每升高0.5℃，年徑流約增加1%左右。因此，未來降水增加、氣溫升高的情景下，總體上徑流是呈增加趨勢的。

表2 不同氣候變化情景下風(fēng)火山流域徑流變化幅度Table 2 Range of runoff change in Fenghuoshan watershed under different climate change scenarios

4 討論

4.1 LSTM模型的凍土水文學(xué)意義

不同于一般的神經(jīng)網(wǎng)絡(luò)模型只能作為黑箱模型用于水文過程模擬，LSTM由于其特殊的細(xì)胞狀態(tài)和門結(jié)構(gòu)，使其具有一定的水文學(xué)意義［28］。LSTM與一般水文模型相似，對于模型降水和氣溫輸入是逐時(shí)間步長處理的，如本研究中以天為時(shí)間步長，每天的降水和氣溫輸入進(jìn)LSTM中后都被用來更新當(dāng)前步長內(nèi)的細(xì)胞狀態(tài)。細(xì)胞狀態(tài)是LSTM中用來存儲(chǔ)長期記憶信息的關(guān)鍵變量，類比于一般水文模型，細(xì)胞狀態(tài)可以理解為積雪深度、土壤含水量、地下水儲(chǔ)量等水文過程中的狀態(tài)變量；而遺忘門、輸入門和輸出門則可以類比理解為積雪深度、土壤含水量、地下水儲(chǔ)量的狀態(tài)變量的消耗、增長和出流［28］。特別地，在凍土地區(qū)，由于包含氣溫作為輸入，而活動(dòng)層的凍融過程是氣溫的函數(shù)，因此細(xì)胞狀態(tài)也可以類比于活動(dòng)層的凍融狀態(tài)，而門結(jié)構(gòu)則控制著活動(dòng)層內(nèi)的能量變化：當(dāng)氣溫大于0℃時(shí)，隨著融化指數(shù)的增加活動(dòng)層融化深度逐漸加深；當(dāng)氣溫小于0℃時(shí)，隨著凍結(jié)指數(shù)的增加活動(dòng)層凍結(jié)深度逐漸加深。從凍土水文學(xué)的角度分析，LSTM模型以降水和氣溫作為輸入，細(xì)胞狀態(tài)同時(shí)體現(xiàn)了流域蓄水量和凍土活動(dòng)層凍融狀態(tài)兩個(gè)狀態(tài)變量的變化情況，即受活動(dòng)層凍融過程影響的流域蓄水量。在模型訓(xùn)練期內(nèi)確定的LSTM參數(shù)和超參數(shù)則類比于一般水文模型的參數(shù)率定過程；而在模型驗(yàn)證期，LSTM與一般水文模型驗(yàn)證一樣，采用訓(xùn)練期內(nèi)已確定的參數(shù)，僅依靠當(dāng)前步長的輸入和當(dāng)前細(xì)胞狀態(tài)來更新細(xì)胞狀態(tài)。然而，LSTM相對于一般水文模型來說，其沒有具有物理意義的數(shù)學(xué)公式來描述凍土水文過程，只能通過數(shù)據(jù)在模型訓(xùn)練期學(xué)習(xí)凍土水文過程特征。

如圖6所示，LSTM中兩個(gè)神經(jīng)元細(xì)胞狀態(tài)值在驗(yàn)證期內(nèi)隨時(shí)間的變化過程類似活動(dòng)層凍融過程［圖6（b）］、土壤含水量［圖6（c）］的年內(nèi)變化過程。圖6（b）中，4月初至5月初，氣溫雖然仍低于0℃，但總體氣溫呈快速上升趨勢，且此時(shí)流域內(nèi)地表已經(jīng)開始逐漸融化，第Ⅰ部分中細(xì)胞狀態(tài)值的下降與表層土壤向下融化深度逐漸加深的過程相一致；5月初，氣溫開始高于0℃，隨著氣溫的升高，活動(dòng)層融化深度逐漸加深，直至9月中旬氣溫降至0℃附近波動(dòng)，活動(dòng)層融化深度達(dá)最大值，而在多年凍土區(qū)活動(dòng)層存在雙向凍結(jié)的過程，活動(dòng)層底部土壤開始逐漸由下向上凍結(jié)，第Ⅱ部分中細(xì)胞狀態(tài)值的變化與此過程相一致；10月初氣溫開始低于0℃且逐漸降低，表層土壤開始凍結(jié)，且凍結(jié)深度隨著氣溫的降低逐漸加深，第Ⅲ部分中細(xì)胞狀態(tài)值的變化與此過程相一致。圖6（c）中，可以明顯地看到細(xì)胞狀態(tài)值在融化、凍結(jié)過程中（0℃附近）呈S型曲線變化，與土壤水分在活動(dòng)層融化過程、凍結(jié)過程中的變化一致［64］。類似的，Kratzert等［28］發(fā)現(xiàn)將LSTM模型應(yīng)用于積雪影響的流域時(shí)，細(xì)胞狀態(tài)值的變化能夠體現(xiàn)積雪、融雪過程，當(dāng)氣溫低于0℃時(shí)細(xì)胞狀態(tài)值開始逐漸增大，直到氣溫升至0℃細(xì)胞狀態(tài)值迅速減小。因此，盡管LSTM僅僅利用降水、氣溫和徑流來訓(xùn)練模型參數(shù)用以模擬凍土區(qū)的徑流過程，但模型學(xué)習(xí)到了活動(dòng)層凍土凍融變化過程及其土壤水分變化過程特征，從而具有了一定的凍土水文學(xué)意義。因此，可以根據(jù)LSTM、凍土水文的特點(diǎn)，進(jìn)一步改進(jìn)LSTM，建立更具凍土水文學(xué)意義的模型，如Chen等［29］依據(jù)短期徑流在水文預(yù)報(bào)中重要作用，引入自注意力機(jī)制改進(jìn)了LSTM，建立了更適用于水文預(yù)報(bào)的SA-LSTM模型。

圖6 模型驗(yàn)證期LSTM神經(jīng)元細(xì)胞狀態(tài)值變化過程Fig.6 Cell state change processes of LSTM neurons during model validation：measured precipitation，air temperature and runoff（a），cell state change processes of Neuron 1（b），and cell state change processes of Neuron 2（c）

4.2 氣溫升高對凍土區(qū)徑流過程的影響

氣溫升高，將導(dǎo)致多年凍土退化［42-44］，使活動(dòng)層加深、融化期延長、凍結(jié)期縮短及地下冰融化，進(jìn)而通過影響地下水補(bǔ)給、徑流路徑和排泄過程及地下水與地表水的交換等方式改變徑流過程［50，65-69］。地下冰融化對徑流的補(bǔ)給有限，Yang等［66］發(fā)現(xiàn)其對徑流的貢獻(xiàn)占13.2%～16.7%，因此凍土退化雖然增加了冬季基流，但氣溫升高導(dǎo)致的地下冰的融化的增量對全年的徑流增加不大（圖7），更多的是通過活動(dòng)層的變化改變產(chǎn)流過程。由表2、圖7可知，降水增加將導(dǎo)致全年徑流增加，而氣溫升高雖然總體上使風(fēng)火山流域年徑流增加，但由于活動(dòng)層的存在，對年內(nèi)不同月份的徑流的影響是不同的。

圖7 不同氣候變化情景下的模擬結(jié)果Fig.7 Simulation results under different climate change scenarios

氣溫升高，一方面加速了積雪的融化從而直接增加了春季的徑流，另一方面改變了活動(dòng)層的凍融過程間接改變了徑流過程。春季融化期，風(fēng)火山流域4月初開始，地表出現(xiàn)緩慢的融化，同時(shí)積雪開始融化，此時(shí)表層融化較淺易于形成飽和狀態(tài)，地表以蓄滿產(chǎn)流為主，氣溫升高，一方面使積雪、地表融化時(shí)間提前，同時(shí)也加速了二者的融化過程，從而使蓄滿產(chǎn)流出現(xiàn)的時(shí)間提前且增加，此時(shí)氣溫主導(dǎo)了流域的產(chǎn)流，增溫對徑流的增幅顯著；5—6月積雪已完全融化，而活動(dòng)層的融化主要發(fā)生在土壤淺層，此時(shí)壤中流開始出現(xiàn)，流域以壤中流、蓄滿產(chǎn)流并存方式產(chǎn)流，此時(shí)氣溫升高將加速淺層土壤的融化從而增加壤中流、流域的最大蓄水容量，降水和冰融化的需補(bǔ)充更多的土壤水分才能使流域蓄滿產(chǎn)流，另一方面由于此時(shí)氣溫已經(jīng)大于0℃，蒸發(fā)不可忽視，氣溫的升高也將增大淺層土壤的蒸發(fā)，因此增溫對徑流的增幅較小，但由于此時(shí)降水較4月增加較多，故氣溫升高導(dǎo)致的徑流增量大于4月；7月，由于土壤融化深度加深、超過60 cm，此時(shí)地表產(chǎn)流以超滲產(chǎn)流為主，氣溫升高對流域超滲產(chǎn)流影響甚微，而地下產(chǎn)流以壤中流和地下水（即凍結(jié)層上水）為主，氣溫升高將加速土壤融化，增加壤中流和地下水，同時(shí)由于土壤融化深度較深，對深層土壤蒸發(fā)影響有限，因此氣溫升高對徑流的增幅較為明顯；8月，活動(dòng)層融化深度繼續(xù)加深，接近活動(dòng)層最大厚度，此時(shí)，產(chǎn)流方式與7月相同，但由于此時(shí)氣溫達(dá)年內(nèi)最大值，蒸發(fā)量大，氣溫的增加將加劇蒸發(fā)，同時(shí)土壤持水能力也大大增加，下滲的降水、冰融化水優(yōu)先補(bǔ)給土壤，氣溫升高將加大這兩部分的耗水，從而將導(dǎo)致徑流的減少，如表2所示。

9月，流域徑流過程進(jìn)入秋季退水階段，活動(dòng)層融化深度達(dá)最大值，同時(shí)由于氣溫開始降低，流域內(nèi)部分區(qū)域地表開始凍結(jié)，隨著凍結(jié)面積的增加，使地表開始出現(xiàn)不透水層，地表出現(xiàn)零厚度包氣帶性質(zhì)的蓄滿產(chǎn)流，而壤中流因地表逐漸凍結(jié)而缺少入滲補(bǔ)給迅速減少并逐漸消失，地下水因?yàn)槎嗄陜鐾羺^(qū)雙向凍結(jié)的存在而隨氣溫的降低而減少直至因活動(dòng)層完全凍結(jié)而消失。一方面，氣溫升高延緩了地表開始凍結(jié)的時(shí)間，使地表產(chǎn)流仍以超滲產(chǎn)流為主，同時(shí)也使壤中流消失的時(shí)間延遲；另一方面，氣溫升高延遲了活動(dòng)層的雙向凍結(jié)過程，延緩了活動(dòng)層導(dǎo)水系數(shù)因凍結(jié)而減小的過程［7］，增加了同時(shí)期的地下水出流。如圖7所示，10月下旬至12月，由于氣溫已經(jīng)低于0℃，上述的氣溫升高導(dǎo)致地下水增幅明顯，為年內(nèi)因氣溫升高導(dǎo)致的徑流增幅之最；而9月初至中旬，由于氣溫仍大于0℃，氣溫升高對上述過程影響較小，且由于氣溫升高也將導(dǎo)致蒸發(fā)增加，故徑流反而略有減小。王根緒等［4］基于風(fēng)火山流域徑流的研究提出了溫度變源產(chǎn)流的概念，即在多年凍土區(qū)產(chǎn)流過程并非由土壤水分條件唯一決定而更多是由溫度條件控制的［9］；本研究中，氣溫升高對徑流不同月份的影響反映了這一規(guī)律。因此，在多年凍土區(qū)溫度變源產(chǎn)流條件下，氣溫升高，一方面促進(jìn)活動(dòng)層的融化過程、減緩凍結(jié)過程而改變產(chǎn)流過程，使春、冬季節(jié)徑流增加，另一方面夏季氣溫的升高也會(huì)促進(jìn)蒸發(fā)，使徑流出現(xiàn)一定的減少。在未來氣候變暖、凍土退化的情況下，最終將使凍土區(qū)流域年內(nèi)徑流過程趨于平緩。

5 結(jié)論

本文依據(jù)多年凍土區(qū)產(chǎn)流機(jī)制，以降水、氣溫作為輸入，建立了基于LSTM的凍土水文模型，探究了其凍土水文學(xué)意義，并分析了氣溫升高對凍土區(qū)徑流過程的影響，得到以下結(jié)論：

（1）LSTM模型的細(xì)胞狀態(tài)能夠存儲(chǔ)凍土水文過程中的長期記憶信息，如活動(dòng)層凍融狀態(tài)、土壤含水量等，而遺忘門、輸入門和輸出門則可以類比理解為活動(dòng)層中的能量、土壤含水量等狀態(tài)變量的消耗、增長和出流。因此，通過降水、氣溫和徑流來訓(xùn)練模型時(shí)，其能夠?qū)W習(xí)活動(dòng)層凍融過程、土壤含水量變化規(guī)律，并且在驗(yàn)證時(shí)能通過細(xì)胞狀態(tài)值反映出這些凍土水文關(guān)鍵變量隨氣溫、降水的變化，從而使模型具有了一定的凍土水文學(xué)意義，并在風(fēng)火山流域表現(xiàn)出了良好的適用性。

（2）將模型應(yīng)用于預(yù)測未來氣候變化情景下風(fēng)火山流域徑流的變化規(guī)律，總體上風(fēng)火山流域的徑流是呈現(xiàn)增加趨勢的，降水每增加10%，年徑流約增加12%；而氣溫升高對年徑流增加不大，每升高0.5℃，年徑流僅增加1%。但是氣溫升高將改變活動(dòng)層的凍融過程，對年內(nèi)不同時(shí)期的徑流過程產(chǎn)生不同的影響，在春季融化期、秋季凍結(jié)期時(shí)徑流增幅明顯，而由于蒸發(fā)加劇、活動(dòng)層加深，徑流在8月出現(xiàn)了減少；同時(shí)，氣溫升高將延長融化期、縮短凍結(jié)期，并且改變徑流組分。例如，在秋季凍結(jié)過程中，本來因凍結(jié)地表出現(xiàn)不透水層，從而形成了零厚度包氣帶性質(zhì)的蓄滿產(chǎn)流，同時(shí)壤中流因地表逐漸凍結(jié)而缺少入滲補(bǔ)給迅速減少并逐漸消失；而因氣溫升高，延緩了凍結(jié)時(shí)間，使地表產(chǎn)流仍以超滲為主，同時(shí)也使壤中流消失的時(shí)間延遲。這些現(xiàn)象反映了溫度變源產(chǎn)流規(guī)律，即在多年凍土區(qū)產(chǎn)流過程并非由土壤水分條件唯一決定而更多是由溫度條件控制的。

綜上所述，基于LSTM的凍土水文模型具有一定的凍土水文學(xué)意義，在缺少土壤溫度、水分觀測資料的條件下，能夠依靠有限的降水、氣溫、徑流資料模擬青藏高原上的多年凍土區(qū)受活動(dòng)層凍融過程影響的徑流過程，為凍土水文模擬研究提供了一種簡單有效、具有一定物理意義的方法。