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

0 引 言

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

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

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

1 研究方法

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

人工神經(jīng)網(wǎng)絡(luò)是具有很強(qiáng)的處理非線性問題能力的數(shù)據(jù)處理系統(tǒng)，適合于解決因果關(guān)系復(fù)雜的非確定性推理等問題，且自身具有自適應(yīng)性能力、容錯性等優(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個輸入到輸出的權(quán)重；βj為連接第j個隱含層到輸出的權(quán)重；γj0為第j個隱含層閥值；γji為連接第i個輸入到第j個隱含層的權(quán)重；Ψ(t)為轉(zhuǎn)移函數(shù)，取用雙曲正切函數(shù)，即：Ψ(t)=tanh(z)。

通常，人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)是固定的，即各神經(jīng)元之間連接的個數(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，此鏈接無效。

使∧為包含全部指示函數(shù)的矢量，使α=(α0,α1,…,αp)，β=(β1,…,βM)，γj=(γj0,γj1,…,γjp),γ=(γ1,…,γM) ，θ=(α∧,β∧,γ∧,σ∧)，其中α∧,β∧,γ∧分別表示α,β,γ中非零子集，因此可以用(θ,∧)來表示貝葉斯神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)。為了對模型進(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)練中，選出一個最能反映觀測數(shù)據(jù)的最佳網(wǎng)絡(luò)結(jié)構(gòu)；貝葉斯方法主要原理為在已知觀測數(shù)據(jù)條件下，指導(dǎo)(θ,∧)的后驗概率分布。在貝葉斯訓(xùn)練結(jié)構(gòu)中，應(yīng)用觀測數(shù)據(jù)和參數(shù)及模型結(jié)構(gòu)的先驗知識來導(dǎo)出模型結(jié)構(gòu)和參數(shù)的后驗分布。已知數(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時，δ(Z)=1，否則δ(Z)=0。

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

(4)

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

1.3 EMC

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

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

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

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

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

本文選用了葉爾羌河卡群站月平均徑流(1954-2008年)，其中1954-2000年資料用來率定，2001-2008年用來驗證統(tǒng)計降尺度模型。本文采用NCEP全球再分析日資料作為觀測的大尺度氣候資料，未來氣候數(shù)據(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ù)測。

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

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

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

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

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

2.3 評價指標(biāo)

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

3 結(jié)果與分析

3.1 模擬能力對比分析

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

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

3.2 未來徑流情景預(yù)測

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

圖3 率定期及驗證期內(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ù)測平均值及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

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

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

4 結(jié) 語

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

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

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

□

[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ù)報湖北汛期降水量的應(yīng)用研究[J]．氣象學(xué)報，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年來徑流演變規(guī)律研究[J]．干旱區(qū)資源與環(huán)境，2012，26(1)：93-97．

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