資料匱乏地區(qū)徑流降尺度模型構(gòu)建及預(yù)測(cè)

夏依拉·居爾艾提 ，楊 濤(河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 南京 210098)

0 引 言

水是干旱區(qū)的命脈，中亞干旱區(qū)四面環(huán)山，徑流補(bǔ)給主要來(lái)自于冰川和積雪，受氣候影響較大，隨著氣候變化的進(jìn)一步加劇，冰川退縮、雪蓋減少，對(duì)水資源產(chǎn)生了重大的影響。

目前研究氣候變化對(duì)水資源影響的常規(guī)方法是: 由全球氣候模型GCM結(jié)合降尺度方法用來(lái)驅(qū)動(dòng)水文模型而產(chǎn)生徑流的計(jì)算結(jié)果。但由于一般流域水文模型需要大量的水文氣象觀測(cè)資料(降水、氣溫、徑流)，而我國(guó)大部分西部地區(qū)，水文氣象觀測(cè)站網(wǎng)非常低，很多中高山帶(2 500 m以上)基本上屬于無(wú)資料地區(qū)，水文氣象資料嚴(yán)重匱乏，故常規(guī)的流域水文模型根本不能使用，迫切需要研究一種簡(jiǎn)便快捷同時(shí)具有一定精度的徑流預(yù)測(cè)模型。徑流降尺度方法簡(jiǎn)單易行，所需資料少，同時(shí)精度不低于常規(guī)方法，適合于資料匱乏地區(qū)的氣候變化影響研究，近來(lái)日益受到國(guó)內(nèi)外學(xué)者的重視。

研究表明，很多方法都可用于徑流降尺度模型構(gòu)建，例如：人工神經(jīng)網(wǎng)絡(luò)(ANN)，廣義線性模型(GLM)，關(guān)聯(lián)向量機(jī)(RVM)，典型相關(guān)分析(CCA)等。Cannon[1]在加拿大不列顛哥倫比亞省21個(gè)流域利用人工神經(jīng)網(wǎng)絡(luò)構(gòu)建徑流降尺度模型，并與線性回歸方法進(jìn)行了對(duì)比，得出人工神經(jīng)網(wǎng)絡(luò)可以較好地捕捉到輸入和輸出之間復(fù)雜的非線性關(guān)系；Ghosh[2]應(yīng)用關(guān)聯(lián)向量機(jī)構(gòu)建統(tǒng)計(jì)降尺度模型，模擬了印度的雨季徑流，并結(jié)合GCM預(yù)測(cè)了未來(lái)氣候情景下的徑流趨勢(shì)。其中人工神經(jīng)網(wǎng)絡(luò)(ANNs)可以很好地捕捉到輸入和輸出之間的非線性關(guān)系并且不需要對(duì)其中基本物理過(guò)程有深入的了解[3]，常被用來(lái)解決很多水文問(wèn)題，但因其無(wú)法給出預(yù)報(bào)量模擬結(jié)果的不確定性，限制了ANNs在解決水文問(wèn)題中的可用性。貝葉斯神經(jīng)網(wǎng)絡(luò)BNNs作為近來(lái)的一種較為新的方法，不僅能算出預(yù)報(bào)量的平均的模擬值，還可以估計(jì)預(yù)測(cè)的置信區(qū)間。最近，貝葉斯神經(jīng)網(wǎng)絡(luò)廣泛地應(yīng)用于模型的訓(xùn)練、選擇和預(yù)測(cè)等。Zhang[4]等用4種類型的貝葉斯神經(jīng)網(wǎng)絡(luò)BNNs來(lái)估計(jì)美國(guó)兩個(gè)流域的徑流模擬的不確定性。

1 研究方法

1.1 神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)

人工神經(jīng)網(wǎng)絡(luò)是具有很強(qiáng)的處理非線性問(wèn)題能力的數(shù)據(jù)處理系統(tǒng)，適合于解決因果關(guān)系復(fù)雜的非確定性推理等問(wèn)題，且自身具有自適應(yīng)性能力、容錯(cuò)性等優(yōu)良性能[5]。神經(jīng)網(wǎng)普遍應(yīng)用于輸入輸出之間復(fù)雜的非線性關(guān)系，輸入矢量xt到目標(biāo)函數(shù)的yi的映射，即yt=f(xt)+ε，其中f(xt)為神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)：

(1)

式中：p為輸入層節(jié)點(diǎn)數(shù)；M為隱含層節(jié)點(diǎn)數(shù)；αb為輸出層閥值；αi為連接第i個(gè)輸入到輸出的權(quán)重；βj為連接第j個(gè)隱含層到輸出的權(quán)重；γj0為第j個(gè)隱含層閥值；γji為連接第i個(gè)輸入到第j個(gè)隱含層的權(quán)重；Ψ(t)為轉(zhuǎn)移函數(shù)，取用雙曲正切函數(shù)，即：Ψ(t)=tanh(z)。

通常，人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)是固定的，即各神經(jīng)元之間連接的個(gè)數(shù)固定的。

1.2 貝葉斯神經(jīng)網(wǎng)絡(luò)

一組指示函數(shù)可以表示特定的連接的有效性，神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)就轉(zhuǎn)變成為：

(2)

令I(lǐng)(ζ)為連接ζ的指示函數(shù)，若I(ζ)=1，此鏈接有效；若I(ζ)=0，此鏈接無(wú)效。

使∧為包含全部指示函數(shù)的矢量，使α=(α0,α1,…,αp)，β=(β1,…,βM)，γj=(γj0,γj1,…,γjp),γ=(γ1,…,γM) ，θ=(α∧,β∧,γ∧,σ∧)，其中α∧,β∧,γ∧分別表示α,β,γ中非零子集，因此可以用(θ,∧)來(lái)表示貝葉斯神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)。為了對(duì)模型進(jìn)行貝葉斯分析，需要指定以下分布：αi～N(0,σ2a)(i=0,…,p)，βj～N(0,σ2β)(j=1,…,M)，γji～N(0,σ2γ)(j=1,…,M)(i=0,…,p)，σ2～I(xiàn)G(v1,v2)，其中σ2α,σ2β,σ2γ,λ是需要人為確定的超參數(shù)，本文中取σ2α=σ2β=σ2γ=5,λ=10。

圖1 貝葉斯神經(jīng)網(wǎng)絡(luò)模型示意圖Fig.1 Schematic of bayesian neural networks

在傳統(tǒng)神經(jīng)網(wǎng)絡(luò)訓(xùn)練中，選出一個(gè)最能反映觀測(cè)數(shù)據(jù)的最佳網(wǎng)絡(luò)結(jié)構(gòu)；貝葉斯方法主要原理為在已知觀測(cè)數(shù)據(jù)條件下，指導(dǎo)(θ,∧)的后驗(yàn)概率分布。在貝葉斯訓(xùn)練結(jié)構(gòu)中，應(yīng)用觀測(cè)數(shù)據(jù)和參數(shù)及模型結(jié)構(gòu)的先驗(yàn)知識(shí)來(lái)導(dǎo)出模型結(jié)構(gòu)和參數(shù)的后驗(yàn)分布。已知數(shù)據(jù)集D={(x1,y1),(x2,y2),…,(xn,yn)},則：

(3)

認(rèn)為人工神經(jīng)網(wǎng)絡(luò)模擬誤差服從均值為0，將標(biāo)準(zhǔn)差為σ的高斯白噪聲分布，則得到的似然函數(shù)即為下式所示的表達(dá)式，可根據(jù)以下表達(dá)式進(jìn)行抽樣：

其中，m=∑pi=0Iαi+∑Mj=1Iβjδ(∑pi=0Iγj1)+∑Mj=1∑pi=0IβjIγji，當(dāng)Z>0時(shí)，δ(Z)=1，否則δ(Z)=0。

貝葉斯神經(jīng)網(wǎng)絡(luò)中，取選出的樣本的平均值為模擬值，即：

(4)

式中：K為所選的模型(θ,∧)的個(gè)數(shù)。

1.3 EMC

本文所選擇的抽樣方法為EMC方法，是結(jié)合了MCMC重要性抽樣、遺傳算法的學(xué)習(xí)能力、平行算法的空間多鏈搜索能力的一種貝葉斯抽樣方法。EMC的一次迭代抽樣中包含三個(gè)操作，分別為變異、交叉、交換，如圖2所示，操作是否有效的判斷服從各自的Metropolis原則[6]。

圖2 一次EMC迭代過(guò)程示意圖Fig.2 Schematic of one iteration of EMC

2 統(tǒng)計(jì)降尺度方法的建立

2.1 研究區(qū)及數(shù)據(jù)

葉爾羌河是塔里木河四源頭之一，發(fā)源于克什米爾北部喀喇昆侖山脈的喀喇昆侖山口，流域范圍東經(jīng)74°28′～80°54′，北緯35°40′～40°31′之間。葉爾羌河流域由于南、西、北三面分別受喀喇昆侖山、帕米爾高原和南天山的阻擋，外來(lái)水汽較少，為典型的干旱大陸性氣候，高山區(qū)的冰雪融水成為河川徑流的最為主要的補(bǔ)給來(lái)源，流域多年平均冰川消融量約占出山口卡群站多年平均徑流量的64.0%，雨雪混合補(bǔ)給占13.4%，地下水補(bǔ)給占22.6%[7]。

本文選用了葉爾羌河卡群站月平均徑流(1954-2008年)，其中1954-2000年資料用來(lái)率定，2001-2008年用來(lái)驗(yàn)證統(tǒng)計(jì)降尺度模型。本文采用NCEP全球再分析日資料作為觀測(cè)的大尺度氣候資料，未來(lái)氣候數(shù)據(jù)來(lái)自世界氣候研究計(jì)劃耦合模式比較計(jì)劃第四代模式(CMIP5)中的幾種模式，根據(jù)相關(guān)研究結(jié)果[8]，本文中選GFDL_ESM2G，GFDL_ESM2M，MIROC5等三種模式的RCP 4.5, RCP 6.0, RCP 8.5等3種情景的2011-2100年月平均數(shù)據(jù)進(jìn)行預(yù)測(cè)。

2.2 氣象預(yù)報(bào)因子的選擇

逐步回歸是多元回歸中用以選擇自變量的一種常用方法。在變量很多時(shí)，其中有的因素可能對(duì)因變量的影響不是很大，而且變量之間不完全相互獨(dú)立的，可能有種種互作關(guān)系，在這種情況下，可以用逐步回歸分析，進(jìn)行因子的篩選。

基于流域范圍內(nèi)6個(gè)NCEP網(wǎng)格的所有待選預(yù)報(bào)因子，應(yīng)用逐步回歸方法(其中設(shè)置置信度α=0.15)進(jìn)行篩選,最終優(yōu)選出徑流降尺度所需要的預(yù)報(bào)因子，如表1所示。

表1 應(yīng)用逐步回歸篩選出來(lái)的預(yù)報(bào)因子Tab.1 The input variables using stepwise selection method

注：air為近地面平均氣溫；tmax為近地面日最高氣溫；shum為近地面比濕；uwnd為地表經(jīng)向風(fēng)；vwnd為地表緯向風(fēng)。

2.3 評(píng)價(jià)指標(biāo)

為了系統(tǒng)評(píng)價(jià)降尺度效果，使用納什效率系數(shù)Ens，相關(guān)系數(shù)R，均方根誤差RMSE來(lái)衡量模型對(duì)徑流量及其年際變化的捕捉能力。應(yīng)用貝葉斯神經(jīng)網(wǎng)絡(luò)，可以估計(jì)預(yù)測(cè)結(jié)果的不確定性區(qū)間。為了定量地評(píng)估不同的隱含層節(jié)點(diǎn)條件下的不確定性區(qū)間，應(yīng)用以下評(píng)價(jià)指標(biāo)：90%不確定性區(qū)間所包含觀測(cè)值的百分比(POC), 平均相對(duì)寬度(ARIL), 預(yù)測(cè)界限關(guān)于相應(yīng)的觀測(cè)值的平均對(duì)稱度(AAD) (若AAD值小于0.5表明，觀測(cè)徑流基本在不確定性區(qū)間之內(nèi)，AAD值越大，不確定性邊界關(guān)于觀測(cè)值越不對(duì)稱)和平均偏離度(ADA)。

3 結(jié)果與分析

3.1 模擬能力對(duì)比分析

以篩選的預(yù)報(bào)因子為神經(jīng)網(wǎng)絡(luò)的輸入，用卡群站月平均徑流作為網(wǎng)絡(luò)的期望輸出,分別用貝葉斯神經(jīng)網(wǎng)絡(luò)和BP神經(jīng)網(wǎng)絡(luò)構(gòu)建徑流降尺度模型，評(píng)價(jià)神經(jīng)網(wǎng)絡(luò)統(tǒng)計(jì)降尺度法對(duì)卡群站月徑流的模擬精度(表2)?？梢钥闯?，BP神經(jīng)網(wǎng)絡(luò)雖然在率定期模擬效果較優(yōu)，但驗(yàn)證期較差，因?yàn)锽P神經(jīng)網(wǎng)絡(luò)訓(xùn)練中，選出一個(gè)最能反映觀測(cè)數(shù)據(jù)的最佳網(wǎng)絡(luò)結(jié)構(gòu)，不一定能很好地反映其他的數(shù)據(jù)序列，而貝葉斯神經(jīng)網(wǎng)絡(luò)取符合條件的多個(gè)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的平均值作為模擬值，彌補(bǔ)了BP神經(jīng)網(wǎng)絡(luò)的不足，即率定期及驗(yàn)證期相關(guān)系數(shù)都達(dá)到0.9以上，效率系數(shù)也達(dá)到0.8。不同的隱含層節(jié)數(shù)對(duì)模擬效果影響不是很大，相對(duì)而言當(dāng)隱含層節(jié)點(diǎn)M=4時(shí)，模擬效果較優(yōu),故本文中取隱含層節(jié)點(diǎn)數(shù)為4。貝葉斯神經(jīng)網(wǎng)絡(luò)最大的優(yōu)點(diǎn)是，不僅能給出平均的模擬值，還可以定量估計(jì)模擬的不確定性區(qū)間(表2)，不同隱含層條件下的不確定區(qū)間，難以同時(shí)滿足較高的覆蓋度、較小的相對(duì)寬度、較大的對(duì)稱性和較小的偏離度，當(dāng)隱含層節(jié)點(diǎn)數(shù)M=4時(shí)，率定期覆蓋度達(dá)到95%以上，驗(yàn)證期達(dá)到91%以上，相對(duì)M=5、6來(lái)說(shuō)比較高，觀測(cè)值基本在90%不確定性區(qū)間內(nèi)；率定期及驗(yàn)證期AAD值都等于0.25，對(duì)稱性也比較好；率定期相對(duì)寬度為2.56及偏離度為73.33，相對(duì)M=5、6來(lái)說(shuō)都比較低；圖3為當(dāng)M=4時(shí)，率定期及驗(yàn)證期卡群站月平均徑流觀測(cè)值，基于貝葉斯神經(jīng)網(wǎng)絡(luò)的平均模擬值及90%不確定性區(qū)間。

表2 兩種模型對(duì)月徑流模擬能力及不確定性區(qū)間的評(píng)估Tab.2 Assessment of the performance in monthly streamflow simulating and the uncertainty intervals for BNNs and ANNs

3.2 未來(lái)徑流情景預(yù)測(cè)

應(yīng)用已建立的貝葉斯神經(jīng)網(wǎng)絡(luò)，選用GFDL_ESM2G,GFDL_ESM2M, MIROC5等三種模式RCP4.5,RCP6.0,RCP8.5三種情景下的輸出結(jié)果為氣候強(qiáng)迫數(shù)據(jù)，對(duì)葉爾羌河未來(lái)3個(gè)時(shí)段[2011-2040年(2020年代),2041-2070年(2050年代)2071-2100年(2080年代)]卡群站的月平均徑流進(jìn)行預(yù)測(cè)。3種模式，3種情景下的徑流預(yù)測(cè)值及90%不確定性區(qū)間如圖4所示。由圖4可見(jiàn)，GFDL_ESM2G,GFDL_ESM2M在相同情景下的月徑流及其不確定性區(qū)間較為相似；MIROC5未來(lái)3個(gè)時(shí)段枯水期(10-4月)月徑流普遍比前兩個(gè)模式大，且其不確定區(qū)間寬度也較大。

圖3 率定期及驗(yàn)證期內(nèi)徑流降尺度模型模擬結(jié)果及90%不確定性區(qū)間Fig.3 The simulations of the streamflow downscaling in both calibration and validation period with 90% uncertainty intervals

圖4 3種模式3種情景下的徑流預(yù)測(cè)平均值及90%不確定性區(qū)間Fig.4 Projected scenarios of streamflow with a 90% uncertainty interval in three future time horizons under three RCPs based on three selected GCMs

在未來(lái)氣候變化情景下，卡群站年徑流隨時(shí)間推移均呈增加的趨勢(shì)，與基準(zhǔn)期(1961-1990年)相比，3種模式3種情景下，未來(lái)3個(gè)時(shí)段年徑流變化范圍分別為75%～92%，83%～110%，88%～127%。不同情景下的徑流預(yù)測(cè)結(jié)果有差異，其中RCP8.5情景下的徑流增加幅度比其他情景較明顯，尤其是在2050年代和2080年代；主要原因是葉爾羌河徑流以冰川補(bǔ)給為主，溫度的逐年代升高，加速了山區(qū)冰雪資源的消融，對(duì)徑流的影響也越來(lái)越顯著。對(duì)選取的3個(gè)模型3個(gè)情景2010-2100年的氣溫?cái)?shù)據(jù)進(jìn)行分析，可知葉爾羌河流域平均氣溫明顯升高，其中RCP8.5情景下氣溫升高的幅度較大，尤其是在2050年代及2080年代；以GFDL_ESM2G為例，RCP4.5情景下，與基準(zhǔn)期(1961-1990年)平均氣溫相比，未來(lái)3個(gè)時(shí)段平均氣溫增加幅度分別為3.8,4.7,4.9 ℃；RCP6.0情景下，未來(lái)三個(gè)時(shí)段氣溫增加幅度為3.1,4.5,5.7 ℃；RCP8.5情景下，未來(lái)三個(gè)時(shí)段氣溫增加幅度為4.1,5.5,7.6 ℃；故RCP8.5情景下，年徑流增加趨勢(shì)較明顯，例如：MIROC5 RCP8.5情景，在2080年代時(shí)段的年徑流增加速率為2.2 m3/(s·a)。

不同月份的徑流存在不同程度的增加趨勢(shì)，其中9-4月份月徑流變化趨勢(shì)相對(duì)不明顯，例如：GFDL_ESM2G RCP4.5情景下，1月份的徑流增加速率為0.27 m3/(s·a)，MIROC5 RCP6.0情景下，10月份的徑流增加速率為0.53 m3/(s·a)；5-8月份變化趨勢(shì)較明顯，例如：GFDL_ESM2G RCP6.0情景下，7月份的徑流增加速率為1.77 m3/(s·a)，GFDL_ESM2M RCP8.5情景下，6月份的徑流增加速率為2.95 m3/(s·a)。

4 結(jié) 語(yǔ)

本文應(yīng)用貝葉斯神經(jīng)網(wǎng)絡(luò)構(gòu)建氣候變化條件下的徑流降尺度模型，預(yù)測(cè)未來(lái)氣候情景下的徑流，定量估計(jì)預(yù)測(cè)結(jié)果的不確定性，為解決資料缺乏地區(qū)在氣候變化條件下的徑流預(yù)測(cè)問(wèn)題提供科學(xué)依據(jù)，基本結(jié)論為:

(1)與一般的人工神經(jīng)網(wǎng)絡(luò)降尺度模型相比，貝葉斯神經(jīng)網(wǎng)絡(luò)降尺度模型對(duì)葉爾羌河卡群站月平均徑流的降尺度效果是較優(yōu)的，其相關(guān)系數(shù)達(dá)到0.9以上，效率系數(shù)達(dá)到0.8；其90%不確定性區(qū)間對(duì)觀測(cè)徑流的覆蓋度達(dá)到90%以上，區(qū)間上下值關(guān)于觀測(cè)值對(duì)稱性較優(yōu)，故可以用貝葉斯神經(jīng)網(wǎng)絡(luò)對(duì)卡群站在未來(lái)氣候情景下的徑流進(jìn)行預(yù)測(cè)。

(2) GFDL_ESM2G，GFDL_ESM2M及MIROC5三種模式,RCP4.5, RCP 6.0, RCP 8.5三種情景下,葉爾羌河卡群站未來(lái)3個(gè)時(shí)段( 2020年代,2050年代,2080年代)的年徑流受氣溫升高影響隨時(shí)間推移呈現(xiàn)增加的趨勢(shì)，增加幅度為75%～92%，83%～110%，88%～127%，其中RCP8.5情景下的徑流增加趨勢(shì)比其他情景較為明顯；不同月份的徑流變化存在不同程度的增加趨勢(shì),其中9-4月份月徑流變化趨勢(shì)相對(duì)不明顯，而5-8月份變化趨勢(shì)較明顯。

□

[1] Cannon A J，Whitfield P H．Downscaling recent streamflow conditions in British Columbia Canada using ensemble neural network models [J]．Journal of Hydrology，2002，259 (1-4)：136-151．

[2] Ghosh S，Mujumdar P P．Statistical downscaling of GCM simulations to streamflow using relevance vector machine [J]．Advances in Water Resources，2008，31(1)：132-146．

[3] Kingston G B et al．Bayesian training of artificial neural networks used for water resources modeling [J]．Water Resources Research，2005，41(12)．

[4] Xuesong Zhang, Faming Liang et al．Estimating uncertainty of streamflow simulation using Bayesian neural networks [J]． Water Resources Research，2009，45．

[5] 胡江林，張禮平，宇如聰．神經(jīng)網(wǎng)絡(luò)模型預(yù)報(bào)湖北汛期降水量的應(yīng)用研究[J]．氣象學(xué)報(bào)，2001，59 (6)：776-783．

[6] Liang F．Bayesian neural networks for non-linear time series forecasting[J]．Statistics and Computing，2005，15：13- 29．

[7] 張賢芳，舒 強(qiáng)，李 偲．葉爾羌河近48年來(lái)徑流演變規(guī)律研究[J]．干旱區(qū)資源與環(huán)境，2012，26(1)：93-97．

[8] 吳 晶．氣候模式對(duì)中國(guó)西北干旱區(qū)氣溫和降水模擬結(jié)果評(píng)價(jià)[D]．烏魯木齊：中國(guó)科學(xué)院研究新疆地理與生態(tài)研究所，2013．