沈嘉聚，楊漢波，劉志武，楊大文

(1.中國長江三峽集團(tuán)有限公司科學(xué)技術(shù)研究院，北京 100038； 2.清華大學(xué)水利水電工程系，北京 100084；3.清華大學(xué)水沙科學(xué)與水利水電工程國家重點實驗室，北京 100084)

由于氣候變化、社會經(jīng)濟(jì)發(fā)展以及長江上游水資源開發(fā)利用的影響，近年來長江上游徑流發(fā)生了較大變化[1-5]。開展長江上游徑流對氣候變化的響應(yīng)研究可以為長江上游水電開發(fā)、梯級水庫優(yōu)化調(diào)度運行提供依據(jù)，對流域水資源規(guī)劃、開發(fā)利用與保護(hù)具有重要意義。目前，徑流變化的分析方法可歸納為3類：水文模型情景模擬法[6-10]、統(tǒng)計模型相關(guān)分析方法[11-14]和基于Budyko框架的彈性系數(shù)法[15-17]。水文模型情景模擬法采用的模型通常具有良好的物理基礎(chǔ)，但模型結(jié)構(gòu)與模型參數(shù)存在一定的不確定性，且對輸入數(shù)據(jù)要求較高[18-20]；統(tǒng)計模型運算簡單，但自變量之間的相關(guān)性往往會影響歸因結(jié)果；基于Budyko框架的彈性系數(shù)法，因其物理意義清楚、所需資料易于獲得而被廣泛應(yīng)用，但長時間序列中的噪聲通常會對歸因結(jié)果產(chǎn)生影響[21-22],因此在歸因分析時通常采用多種分析方法進(jìn)行對比驗證[7]。根據(jù)水文模型情景模擬法中采用模型的建模原理可以分為過程驅(qū)動模型和數(shù)據(jù)驅(qū)動模型，過程驅(qū)動模型在不同區(qū)域參數(shù)需要分別率定，可能產(chǎn)生較大不確定性，且應(yīng)用較為復(fù)雜；隨著數(shù)據(jù)科學(xué)與高性能計算的快速發(fā)展，基于機(jī)器學(xué)習(xí)方法的數(shù)據(jù)驅(qū)動模型能快速從大量數(shù)據(jù)中訓(xùn)練學(xué)習(xí)得到有效信息，并用于未來變化情景的預(yù)測，因而受到廣泛關(guān)注[23]。BP神經(jīng)網(wǎng)絡(luò)是機(jī)器學(xué)習(xí)中較為成熟的一種算法，通過反向傳播算法訓(xùn)練由輸入層、輸出層及若干隱藏層的節(jié)點相互連接而成的多層前饋網(wǎng)絡(luò)，具有很強(qiáng)的非線性映射能力，但其初始權(quán)值與閾值為隨機(jī)選取，容易在網(wǎng)絡(luò)訓(xùn)練過程中陷入局部最小值[24]；而遺傳算法則是模仿生物演替過程中的選擇、交叉、變異和優(yōu)勝劣汰的篩選過程，具有很好的全局搜索性，能夠為神經(jīng)網(wǎng)絡(luò)中不確定的初始權(quán)值與閾值尋找全局最優(yōu)解[25]，從而提高模型準(zhǔn)確性。

在全球變暖的大背景下，氣候變化對長江上游徑流產(chǎn)生了重要影響[26]。以往對于長江上游氣象要素對徑流影響的研究主要通過統(tǒng)計方法[13-14，26-27]進(jìn)行定性分析，以及利用傳統(tǒng)水文模型情景模擬法[7-8，10]進(jìn)行定量歸因，而運用機(jī)器學(xué)習(xí)算法通過情景模擬進(jìn)行敏感性分析的研究較少。因此，本文運用遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)(genetic algorithm optimize BP neural network，GA-BP)算法替代傳統(tǒng)水文模型，分析長江上游徑流對降水與氣溫變化的敏感性，并與多元線性回歸法及基于Budyko框架法的結(jié)果進(jìn)行對比，探究近年來長江上游年徑流量變化趨勢及其氣象驅(qū)動因子。

1 研究區(qū)概況與研究數(shù)據(jù)

1.1 研究區(qū)概況

長江上游是指宜昌水文站控制的流域，包括金沙江、岷江、烏江、嘉陵江等水系的流域區(qū)域(圖1)，流域面積約100萬km2。長江上游區(qū)域地形復(fù)雜，包括青藏高原、四川盆地，涉及青海、云南、甘肅、四川、貴州等多個省份，氣候受青藏高原、東南季風(fēng)以及西南季風(fēng)的綜合影響，區(qū)域內(nèi)空間差異較大。長江上游年均氣溫10～15℃，呈自西向東遞增的分布規(guī)律；年平均降水量800～1 200 mm，呈自西北向東南遞增的分布規(guī)律。除源區(qū)、高山區(qū)域多降雪，長江上游徑流以降水補(bǔ)給為主，徑流主要集中在夏季，雨熱同期，洪水具有洪峰高、洪量大、持續(xù)時間長等特點。

圖1 研究區(qū)域及水文站點分布Fig.1 Study area and distribution of hydrological stations

1.2 研究數(shù)據(jù)

收集了長江上游8個水文站(直門達(dá)、石鼓、屏山、高場、北碚、武隆、寸灘和宜昌)1979—2015年逐日徑流數(shù)據(jù)。同時期的氣象數(shù)據(jù)來源于國家青藏高原科學(xué)數(shù)據(jù)中心(National Tibetan Plateau Third Pole Environment Data Center, TPDC)的中國區(qū)域地面氣象要素驅(qū)動數(shù)據(jù)集(China Meteorological Forcing Dataset, CMFD)[28-29]，包括近地面氣溫、地面向下短波輻射與降水強(qiáng)度，時間分辨率為1 d，水平空間分辨率為0.1°。潛在蒸散發(fā)量ET0采用國際糧農(nóng)組織(FAO)推薦的潛在蒸散發(fā)計算方法[30]簡化參數(shù)推導(dǎo)得到的Irmak-Allen公式[31](式(1))進(jìn)行估算，該公式在Penman-Monteith公式基礎(chǔ)上根據(jù)美國濕潤地區(qū)的數(shù)據(jù)進(jìn)行線性簡化得到[31]，在長江流域有較好的適用性[32-33]。

ET0=-0.611+0.149RS+0.079T

(1)

式中：RS為太陽輻射量，MJ/m2；T為日平均氣溫，℃。

2 研究方法

2.1 GA-BP算法

GA-BP算法的具體流程如圖2所示。神經(jīng)網(wǎng)絡(luò)輸入變量為當(dāng)日降水量、前6 d降水量和日平均氣溫，輸出變量為日徑流量，神經(jīng)網(wǎng)絡(luò)由輸入層、隱藏層和輸出層構(gòu)成，各層神經(jīng)元數(shù)量分別為8、16、1。與水文模型情景模擬法類似，前80%的數(shù)據(jù)為訓(xùn)練集，剩余20%數(shù)據(jù)為測試集。采用納什效率系數(shù)(NSE)、相關(guān)系數(shù)(R2)和均方根誤差(RMSE)評估GA-BP算法的模擬效果。計算徑流對氣象要素變化敏感性時，對測試集輸入變量施加一微小變化量，考慮實際降水與氣溫的變化，取該氣象要素的3%進(jìn)行擾動，得到新的輸出。敏感性系數(shù)為輸出變量的變化量與輸入變量的變化量的比值。

圖2 GA-BP算法流程Fig.2 GA-BP algorithm flowchart

2.2 多元線性回歸法

趨勢分析采用多元線性回歸法，顯著性水平p取0.01，該方法通過多個自變量進(jìn)行線性組合共同預(yù)估因變量。假設(shè)年徑流R與年降水P、年氣溫Ty之間存在線性相關(guān)關(guān)系：

R=β0+βPP+βTTy+ε

(2)

式中：β0為回歸常數(shù)；βP、βT分別為徑流對降水和氣溫變化的敏感性系數(shù)；ε為回歸殘差。

2.3 Budyko框架法

Budyko框架法假設(shè)流域?qū)嶋H蒸散發(fā)量可以表示為流域降水量和潛在蒸散發(fā)量的函數(shù)[34]，據(jù)此Choudhury等[35-36]提出了Choudhury-Yang公式：

(3)

式中：ET為實際蒸散發(fā)量，mm/a；n為下墊面參數(shù)，綜合反映了地形、土壤、植被等因素影響。結(jié)合多年流域水量平衡方程：

R=P-ET

(4)

n可根據(jù)R、P、ET0的多年平均值反算求得，并假定在研究期間不變。將年徑流量R的變化趨勢表示為全微分形式：

(5)

3 結(jié)果與分析

3.1 徑流、降水及氣溫的變化趨勢

表1為長江上游各水文站1979—2015年徑流、降水量及氣溫的變化趨勢?？梢?，徑流增加趨勢最大的為直門達(dá)站，達(dá)0.69 mm/a；減少趨勢最大的為高場站，為2.81 mm/a。沿長江干流從上向下，徑流變化由直門達(dá)站、石鼓站、屏山站的增加趨勢，逐漸轉(zhuǎn)變?yōu)楦邎稣?、北碚站、寸灘站、武隆站和宜昌站的下降趨勢；值得注意的是，在p=0.01的顯著性下，除高場站外其余站點變化趨勢均不顯著。降水量方面，直門達(dá)站和石鼓站控制流域顯著增加(p<0.01)，分別為4.28 mm/a和3.13 mm/a；屏山、寸灘和宜昌站點控制流域呈不顯著增加，其余站點控制流域呈不顯著減少。氣溫方面，各站點控制流域內(nèi)氣溫均呈現(xiàn)顯著上升趨勢(p<0.01)，變化范圍從武隆站的0.021℃/a到直門達(dá)站的0.062℃/a。

表1 長江上游各水文站控制流域1979—2015年徑流、降水及氣溫的變化趨勢Table 1 Trends in runoff, precipitation and air temperature at hydrological stations in upper reaches of the Yangtze River from 1979 to 2015

3.2 GA-BP算法徑流模擬

表2為GA-BP算法徑流模擬結(jié)果。8個水文站訓(xùn)練集的NSE均值為0.71，R2均值為0.71，RMSE均值為0.59 mm/d；測試集NSE均值為0.61，R2均值為0.64，RMSE均值為0.55 mm/d，整體模擬效果良好。直門達(dá)站、石鼓站以及宜昌站模擬效果最好，訓(xùn)練集NSE值均高于0.7，測試集NSE值均高于0.6；北碚站和武隆站的模擬效果相對較差，主要是因為日徑流數(shù)據(jù)波動較大，GA-BP算法對洪峰模擬效果較差，從而拉低整體模擬效果。

表2 GA-BP算法徑流模擬結(jié)果Table 2 Results of GA-BP algorithm runoff simulation

3.3 徑流對氣象要素變化的敏感性

表3為3種方法得到的長江上游各水文站徑流對降水、氣溫變化的敏感性系數(shù)。對于徑流對降水的敏感性，多元線性回歸法得到的各站點敏感性系數(shù)均值為0.65 mm/mm，GA-BP算法和Budyko框架法得到的結(jié)果分別為0.52 mm/mm和0.72 mm/mm；對于徑流對氣溫的敏感性，多元線性回歸法、GA-BP算法以及Budyko框架法的各站點敏感性系數(shù)均值分別為-20.94 mm/℃、-17.99 mm/℃和-18.96 mm/℃。3種方法計算徑流對降水的敏感性均值最小的是直門達(dá)站，為0.37 mm/mm，最大的是高場站，為0.76 mm/mm；對氣溫的敏感性均值最小的時直門達(dá)站，為-4.77 mm/℃，最大的是高場站，為-33.53 mm/℃。沿長江干流從上向下，各水文站徑流對降水的敏感性總體呈現(xiàn)增加趨勢，即下游區(qū)域徑流對降水的變化更為敏感，徑流對氣溫的敏感性同樣呈現(xiàn)增加趨勢。

表3 長江上游各水文站徑流對降水和氣溫的敏感性系數(shù)Table 3 Sensitivity coefficient of runoff of hydrological stations in upper reaches of the Yangtze River to precipitation and temperature

3.4 氣象要素對徑流變化的貢獻(xiàn)量

表4為長江上游各水文站降水和氣溫對徑流變化的貢獻(xiàn)量，可見，多元線性回歸法、GA-BP算法和Budyko框架法得到的降水導(dǎo)致徑流變化的范圍分別為-0.89～1.19 mm/a、-0.55～1.54 mm/a和-0.85～2.10 mm/a；氣溫導(dǎo)致徑流變化的范圍分別為-1.51～0.02 mm/a、-1.26～-0.01 mm/a和-1.11～-0.65 mm/a。與多元線性回歸法相比，GA-BP算法得到的降水和氣溫變化導(dǎo)致徑流的變化，在徑流增加的石門達(dá)站和石鼓站結(jié)果偏大，在徑流減少的高場站、北碚站、寸灘站、武隆站和宜昌站結(jié)果偏小。3種方法得到的降水對徑流變化的貢獻(xiàn)量均值的絕對值最小的是宜昌站，為0.15 mm/a，最大的是直門達(dá)站，為1.58 mm/a；氣溫對徑流變化的貢獻(xiàn)量均值的絕對值最小的直門達(dá)站，為-0.30 mm/a，最大的時北碚站，為-1.17 mm/a。沿長江干流從上到下各水文站降水對徑流變化的貢獻(xiàn)量整體上呈現(xiàn)減少趨勢，氣溫對徑流變化的貢獻(xiàn)量呈現(xiàn)增加趨勢。對于降水和氣溫對徑流變化的貢獻(xiàn)量，定義絕對值大者為主導(dǎo)因素?？梢园l(fā)現(xiàn)，直門達(dá)站、石鼓站、屏山站降水對徑流變化的貢獻(xiàn)量大于氣溫對徑流變化的貢獻(xiàn)量，而高場站、北碚站、寸灘站、武隆站以及宜昌站氣溫對徑流變化的貢獻(xiàn)量逐漸超過降水的貢獻(xiàn)量，成為徑流變化趨勢的主導(dǎo)因素。

表4 長江上游各水文站降水和氣溫對徑流變化的貢獻(xiàn)量Table 4 Contribution of precipitation and temperature at hydrologic stations in upper reaches of the Yangtze River to runoff change

4 討 論

4.1 GA-BP算法的敏感性

采用GA-BP算法進(jìn)行徑流對氣象要素變化的敏感性分析時，輸入變量的微小擾動都會影響敏感性系數(shù)的計算結(jié)果。為探究GA-BP算法計算徑流對氣象要素變化的敏感性結(jié)果的穩(wěn)定性，降水?dāng)_動分別取測試集降水量均值的1%、2%、3%、5%和10%，氣溫擾動分別取測試集氣溫均值的1%、2%、3%、5%和10%，計算不同擾動得到的徑流對氣象要素變化的敏感性，結(jié)果見表5?？梢钥闯觯S著擾動程度的增加，GA-BP算法得到的徑流對降水與氣溫的敏感性系數(shù)逐漸變小，且變化幅度同樣變小并趨于穩(wěn)定。考慮到采用較小擾動時計算誤差對于敏感性結(jié)果影響較大，并且氣象要素在1979—2015年的實際變化程度平均約3%，而且擾動在3%時其機(jī)器學(xué)習(xí)敏感性結(jié)果的波動幅度較小且趨于穩(wěn)定，因此在利用GA-BP算法計算敏感性時，選擇3%的擾動進(jìn)行分析。

表5 不同程度擾動GA-BP算法敏感性計算結(jié)果Table 5 Sensitivity calculation results of GA-BP algorithm with different disturbances

為研究GA-BP算法利用不同時間尺度數(shù)據(jù)進(jìn)行敏感性分析的適用性，分別采用日尺度和月尺度數(shù)據(jù)進(jìn)行計算，結(jié)果見表6?？梢姡粘叨扰c月尺度計算結(jié)果中個別站點有所差異，但總體相差不大，計算結(jié)果穩(wěn)定，所有站點的敏感性系數(shù)均值近乎相等。這意味著GA-BP算法不僅可以利用日尺度數(shù)據(jù)計算徑流對降水與氣溫變化的敏感性，也可利用月尺度數(shù)據(jù)進(jìn)行計算，從而可以在缺少日數(shù)據(jù)的流域進(jìn)行應(yīng)用。

表6 不同時間尺度數(shù)據(jù)GA-BP算法敏感性計算結(jié)果Table 6 Sensitivity calculation results of GA-BP algorithm using data of different time scales

4.2 長江上游徑流敏感性與干旱指數(shù)的關(guān)系

沿長江干流從上向下，各水文站徑流對降水與氣溫的敏感性總體呈增加趨勢，高場站徑流對氣溫的敏感性偏高，以下進(jìn)一步探究這種分布形成的原因。直門達(dá)站、石鼓站、屏山站、高場站、北碚站、寸灘站、武隆站、宜昌站的干旱指數(shù)(ET0/P)分別為1.75、1.60、1.28、0.82、0.92、1.08、0.73和1.00，可見，高場站控制的流域干旱指數(shù)明顯低于除武隆外的其他流域。圖3為長江上游各水文站徑流敏感性與干旱指數(shù)的關(guān)系。如圖3(a)所示，徑流對降水的敏感性與干旱指數(shù)間存在很強(qiáng)的線性關(guān)系，敏感性系數(shù)隨著干旱指數(shù)的增大而減??；如圖3(b)所示，徑流對氣溫的敏感性隨著干旱指數(shù)的增大而減小。對于氣象要素的貢獻(xiàn)量，如圖3(c)所示，降水對徑流變化的貢獻(xiàn)量與干旱指數(shù)之間也存在顯著的線性關(guān)系，敏感性系數(shù)隨著干旱指數(shù)的增大而增大；如圖3(d)所示，氣溫對徑流變化的貢獻(xiàn)量隨著干旱指數(shù)的增大而減小。即降水對徑流變化的影響在越干旱的區(qū)域影響越大，而氣溫對徑流變化的影響在越濕潤的區(qū)域越大，這與之前研究結(jié)果一致[37]。

(a) 徑流對降水的敏感性

(b) 徑流對氣溫的敏感性

(c) 降水對徑流變化貢獻(xiàn)量

(d) 氣溫對徑流變化貢獻(xiàn)量圖3 長江上游各水文站徑流敏感性與干旱指數(shù)的關(guān)系Fig.3 Relationship between runoff sensitivity of hydrologic stations and drought index in upper reaches of the Yangtze River

4.3 不確定性分析

GA-BP算法對徑流過程模擬效果較好，但對洪峰的模擬效果較差，NSE指標(biāo)表現(xiàn)一般。以往研究運用機(jī)器學(xué)習(xí)法對長江上游進(jìn)行徑流模擬時，NSE值一般為0.6～0.8，如熊一橙等[38]利用LSTM算法模擬北碚站、高場站徑流的NSE值分別為0.77、0.62；Zhu等[39]利用ANN算法模擬北碚站、高場站日徑流NSE值為0.78、0.74；黃鈺瀚[40]利用VIC模型對寸灘站日徑流進(jìn)行模擬，NSE值為0.80。對比以往研究，本文GA-BP算法對長江上游日徑流模擬的結(jié)果相差不大。除此之外，在計算年徑流量對年降水量與年均氣溫的敏感性系數(shù)時，日徑流或洪峰的模擬效果不佳對于年徑流量的計算結(jié)果影響不大，因此可以認(rèn)為該方法盡管存在一定的不確定性，但得出的長江上游流域徑流變化的敏感性總體是可信的。

本文計算得到的降水與氣溫對徑流變化的貢獻(xiàn)量，與實測值之間存在一定差異。導(dǎo)致這一差異的可能原因是只考慮了降水與氣溫對徑流變化的影響，對其他氣象要素、下墊面植被的變化以及人類活動的影響并未加以考慮。從結(jié)果來看，徑流增加的直門達(dá)站和石鼓站，實測的增加值小于估計值；徑流減少的高場站、北碚站、寸灘站、武隆站和宜昌站，實測的減少值大于估計值，可能原因是人類活動及植被變化導(dǎo)致耗水的增加、徑流減少。

5 結(jié) 論

a.基于GA-BP算法的敏感性分析，利用日徑流數(shù)據(jù)得到的結(jié)果整體較好，與多元線性回歸法以及Budyko框架法的結(jié)果總體一致；對于缺少日徑流數(shù)據(jù)的流域，該方法也可利用月數(shù)據(jù)進(jìn)行計算。

b.年徑流量方面，高場站顯著減少(p<0.01)，趨勢為-2.81 mm/a，其余站點無顯著變化；年降水量方面，直門達(dá)站和石鼓站顯著增加(p<0.01)，趨勢分別為4.28 mm/a和3.13 mm/a，其余站點無顯著變化；年氣溫方面，各站點均顯著上升，幅度為0.021～0.062℃/a。

c.沿長江干流自上而下，徑流變化對降水的敏感性系數(shù)由直門達(dá)站的0.37 mm/mm逐漸增大至宜昌站的0.74 mm/mm，最高為高場站的0.76 mm/mm；而徑流變化對氣溫的敏感性也呈現(xiàn)增大趨勢，范圍從-4.77～-33.53 mm/℃；敏感性與干旱指數(shù)相關(guān)，隨著干旱指數(shù)的增加，徑流變化對降水和氣溫的敏感性減小。

d.從對徑流變化的貢獻(xiàn)量來看，降水的貢獻(xiàn)量沿干流自上而下逐漸變小，氣溫的貢獻(xiàn)量無明顯規(guī)律。其中直門達(dá)站、石鼓站、屏山站徑流的變化主要由降水主導(dǎo)；向下至高場站、北碚站、寸灘站、武隆站、宜昌站，氣溫的貢獻(xiàn)量超過了降水，居于主導(dǎo)地位。呈現(xiàn)出隨著干旱指數(shù)增加降水對徑流變化的影響增大、氣溫對徑流變化的影響減小的規(guī)律。