沈嘉聚,楊漢波,劉志武,楊大文
(1.中國(guó)長(zhǎng)江三峽集團(tuán)有限公司科學(xué)技術(shù)研究院,北京 100038; 2.清華大學(xué)水利水電工程系,北京 100084;3.清華大學(xué)水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100084)
由于氣候變化、社會(huì)經(jīng)濟(jì)發(fā)展以及長(zhǎng)江上游水資源開發(fā)利用的影響,近年來長(zhǎng)江上游徑流發(fā)生了較大變化[1-5]。開展長(zhǎng)江上游徑流對(duì)氣候變化的響應(yīng)研究可以為長(zhǎng)江上游水電開發(fā)、梯級(jí)水庫優(yōu)化調(diào)度運(yùn)行提供依據(jù),對(duì)流域水資源規(guī)劃、開發(fā)利用與保護(hù)具有重要意義。目前,徑流變化的分析方法可歸納為3類:水文模型情景模擬法[6-10]、統(tǒng)計(jì)模型相關(guān)分析方法[11-14]和基于Budyko框架的彈性系數(shù)法[15-17]。水文模型情景模擬法采用的模型通常具有良好的物理基礎(chǔ),但模型結(jié)構(gòu)與模型參數(shù)存在一定的不確定性,且對(duì)輸入數(shù)據(jù)要求較高[18-20];統(tǒng)計(jì)模型運(yùn)算簡(jiǎn)單,但自變量之間的相關(guān)性往往會(huì)影響歸因結(jié)果;基于Budyko框架的彈性系數(shù)法,因其物理意義清楚、所需資料易于獲得而被廣泛應(yīng)用,但長(zhǎng)時(shí)間序列中的噪聲通常會(huì)對(duì)歸因結(jié)果產(chǎn)生影響[21-22],因此在歸因分析時(shí)通常采用多種分析方法進(jìn)行對(duì)比驗(yàn)證[7]。根據(jù)水文模型情景模擬法中采用模型的建模原理可以分為過程驅(qū)動(dòng)模型和數(shù)據(jù)驅(qū)動(dòng)模型,過程驅(qū)動(dòng)模型在不同區(qū)域參數(shù)需要分別率定,可能產(chǎn)生較大不確定性,且應(yīng)用較為復(fù)雜;隨著數(shù)據(jù)科學(xué)與高性能計(jì)算的快速發(fā)展,基于機(jī)器學(xué)習(xí)方法的數(shù)據(jù)驅(qū)動(dòng)模型能快速從大量數(shù)據(jù)中訓(xùn)練學(xué)習(xí)得到有效信息,并用于未來變化情景的預(yù)測(cè),因而受到廣泛關(guān)注[23]。BP神經(jīng)網(wǎng)絡(luò)是機(jī)器學(xué)習(xí)中較為成熟的一種算法,通過反向傳播算法訓(xùn)練由輸入層、輸出層及若干隱藏層的節(jié)點(diǎn)相互連接而成的多層前饋網(wǎng)絡(luò),具有很強(qiáng)的非線性映射能力,但其初始權(quán)值與閾值為隨機(jī)選取,容易在網(wǎng)絡(luò)訓(xùn)練過程中陷入局部最小值[24];而遺傳算法則是模仿生物演替過程中的選擇、交叉、變異和優(yōu)勝劣汰的篩選過程,具有很好的全局搜索性,能夠?yàn)樯窠?jīng)網(wǎng)絡(luò)中不確定的初始權(quán)值與閾值尋找全局最優(yōu)解[25],從而提高模型準(zhǔn)確性。
在全球變暖的大背景下,氣候變化對(duì)長(zhǎng)江上游徑流產(chǎn)生了重要影響[26]。以往對(duì)于長(zhǎng)江上游氣象要素對(duì)徑流影響的研究主要通過統(tǒng)計(jì)方法[13-14,26-27]進(jìn)行定性分析,以及利用傳統(tǒng)水文模型情景模擬法[7-8,10]進(jìn)行定量歸因,而運(yùn)用機(jī)器學(xué)習(xí)算法通過情景模擬進(jìn)行敏感性分析的研究較少。因此,本文運(yùn)用遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)(genetic algorithm optimize BP neural network,GA-BP)算法替代傳統(tǒng)水文模型,分析長(zhǎng)江上游徑流對(duì)降水與氣溫變化的敏感性,并與多元線性回歸法及基于Budyko框架法的結(jié)果進(jìn)行對(duì)比,探究近年來長(zhǎng)江上游年徑流量變化趨勢(shì)及其氣象驅(qū)動(dòng)因子。
長(zhǎng)江上游是指宜昌水文站控制的流域,包括金沙江、岷江、烏江、嘉陵江等水系的流域區(qū)域(圖1),流域面積約100萬km2。長(zhǎng)江上游區(qū)域地形復(fù)雜,包括青藏高原、四川盆地,涉及青海、云南、甘肅、四川、貴州等多個(gè)省份,氣候受青藏高原、東南季風(fēng)以及西南季風(fēng)的綜合影響,區(qū)域內(nèi)空間差異較大。長(zhǎng)江上游年均氣溫10~15℃,呈自西向東遞增的分布規(guī)律;年平均降水量800~1 200 mm,呈自西北向東南遞增的分布規(guī)律。除源區(qū)、高山區(qū)域多降雪,長(zhǎng)江上游徑流以降水補(bǔ)給為主,徑流主要集中在夏季,雨熱同期,洪水具有洪峰高、洪量大、持續(xù)時(shí)間長(zhǎng)等特點(diǎn)。
圖1 研究區(qū)域及水文站點(diǎn)分布Fig.1 Study area and distribution of hydrological stations
收集了長(zhǎng)江上游8個(gè)水文站(直門達(dá)、石鼓、屏山、高場(chǎng)、北碚、武隆、寸灘和宜昌)1979—2015年逐日徑流數(shù)據(jù)。同時(shí)期的氣象數(shù)據(jù)來源于國(guó)家青藏高原科學(xué)數(shù)據(jù)中心(National Tibetan Plateau Third Pole Environment Data Center, TPDC)的中國(guó)區(qū)域地面氣象要素驅(qū)動(dòng)數(shù)據(jù)集(China Meteorological Forcing Dataset, CMFD)[28-29],包括近地面氣溫、地面向下短波輻射與降水強(qiáng)度,時(shí)間分辨率為1 d,水平空間分辨率為0.1°。潛在蒸散發(fā)量ET0采用國(guó)際糧農(nóng)組織(FAO)推薦的潛在蒸散發(fā)計(jì)算方法[30]簡(jiǎn)化參數(shù)推導(dǎo)得到的Irmak-Allen公式[31](式(1))進(jìn)行估算,該公式在Penman-Monteith公式基礎(chǔ)上根據(jù)美國(guó)濕潤(rùn)地區(qū)的數(shù)據(jù)進(jìn)行線性簡(jiǎn)化得到[31],在長(zhǎng)江流域有較好的適用性[32-33]。
ET0=-0.611+0.149RS+0.079T
(1)
式中:RS為太陽輻射量,MJ/m2;T為日平均氣溫,℃。
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ù)為測(cè)試集。采用納什效率系數(shù)(NSE)、相關(guān)系數(shù)(R2)和均方根誤差(RMSE)評(píng)估GA-BP算法的模擬效果。計(jì)算徑流對(duì)氣象要素變化敏感性時(shí),對(duì)測(cè)試集輸入變量施加一微小變化量,考慮實(shí)際降水與氣溫的變化,取該氣象要素的3%進(jìn)行擾動(dòng),得到新的輸出。敏感性系數(shù)為輸出變量的變化量與輸入變量的變化量的比值。
圖2 GA-BP算法流程Fig.2 GA-BP algorithm flowchart
趨勢(shì)分析采用多元線性回歸法,顯著性水平p取0.01,該方法通過多個(gè)自變量進(jìn)行線性組合共同預(yù)估因變量。假設(shè)年徑流R與年降水P、年氣溫Ty之間存在線性相關(guān)關(guān)系:
R=β0+βPP+βTTy+ε
(2)
式中:β0為回歸常數(shù);βP、βT分別為徑流對(duì)降水和氣溫變化的敏感性系數(shù);ε為回歸殘差。
Budyko框架法假設(shè)流域?qū)嶋H蒸散發(fā)量可以表示為流域降水量和潛在蒸散發(fā)量的函數(shù)[34],據(jù)此Choudhury等[35-36]提出了Choudhury-Yang公式:
(3)
式中:ET為實(shí)際蒸散發(fā)量,mm/a;n為下墊面參數(shù),綜合反映了地形、土壤、植被等因素影響。結(jié)合多年流域水量平衡方程:
R=P-ET
(4)
n可根據(jù)R、P、ET0的多年平均值反算求得,并假定在研究期間不變。將年徑流量R的變化趨勢(shì)表示為全微分形式:
(5)
表1為長(zhǎng)江上游各水文站1979—2015年徑流、降水量及氣溫的變化趨勢(shì)??梢姡瑥搅髟黾于厔?shì)最大的為直門達(dá)站,達(dá)0.69 mm/a;減少趨勢(shì)最大的為高場(chǎng)站,為2.81 mm/a。沿長(zhǎng)江干流從上向下,徑流變化由直門達(dá)站、石鼓站、屏山站的增加趨勢(shì),逐漸轉(zhuǎn)變?yōu)楦邎?chǎng)站、北碚站、寸灘站、武隆站和宜昌站的下降趨勢(shì);值得注意的是,在p=0.01的顯著性下,除高場(chǎng)站外其余站點(diǎn)變化趨勢(shì)均不顯著。降水量方面,直門達(dá)站和石鼓站控制流域顯著增加(p<0.01),分別為4.28 mm/a和3.13 mm/a;屏山、寸灘和宜昌站點(diǎn)控制流域呈不顯著增加,其余站點(diǎn)控制流域呈不顯著減少。氣溫方面,各站點(diǎn)控制流域內(nèi)氣溫均呈現(xiàn)顯著上升趨勢(shì)(p<0.01),變化范圍從武隆站的0.021℃/a到直門達(dá)站的0.062℃/a。
表1 長(zhǎng)江上游各水文站控制流域1979—2015年徑流、降水及氣溫的變化趨勢(shì)Table 1 Trends in runoff, precipitation and air temperature at hydrological stations in upper reaches of the Yangtze River from 1979 to 2015
表2為GA-BP算法徑流模擬結(jié)果。8個(gè)水文站訓(xùn)練集的NSE均值為0.71,R2均值為0.71,RMSE均值為0.59 mm/d;測(cè)試集NSE均值為0.61,R2均值為0.64,RMSE均值為0.55 mm/d,整體模擬效果良好。直門達(dá)站、石鼓站以及宜昌站模擬效果最好,訓(xùn)練集NSE值均高于0.7,測(cè)試集NSE值均高于0.6;北碚站和武隆站的模擬效果相對(duì)較差,主要是因?yàn)槿諒搅鲾?shù)據(jù)波動(dòng)較大,GA-BP算法對(duì)洪峰模擬效果較差,從而拉低整體模擬效果。
表2 GA-BP算法徑流模擬結(jié)果Table 2 Results of GA-BP algorithm runoff simulation
表3為3種方法得到的長(zhǎng)江上游各水文站徑流對(duì)降水、氣溫變化的敏感性系數(shù)。對(duì)于徑流對(duì)降水的敏感性,多元線性回歸法得到的各站點(diǎn)敏感性系數(shù)均值為0.65 mm/mm,GA-BP算法和Budyko框架法得到的結(jié)果分別為0.52 mm/mm和0.72 mm/mm;對(duì)于徑流對(duì)氣溫的敏感性,多元線性回歸法、GA-BP算法以及Budyko框架法的各站點(diǎn)敏感性系數(shù)均值分別為-20.94 mm/℃、-17.99 mm/℃和-18.96 mm/℃。3種方法計(jì)算徑流對(duì)降水的敏感性均值最小的是直門達(dá)站,為0.37 mm/mm,最大的是高場(chǎng)站,為0.76 mm/mm;對(duì)氣溫的敏感性均值最小的時(shí)直門達(dá)站,為-4.77 mm/℃,最大的是高場(chǎng)站,為-33.53 mm/℃。沿長(zhǎng)江干流從上向下,各水文站徑流對(duì)降水的敏感性總體呈現(xiàn)增加趨勢(shì),即下游區(qū)域徑流對(duì)降水的變化更為敏感,徑流對(duì)氣溫的敏感性同樣呈現(xiàn)增加趨勢(shì)。
表3 長(zhǎng)江上游各水文站徑流對(duì)降水和氣溫的敏感性系數(shù)Table 3 Sensitivity coefficient of runoff of hydrological stations in upper reaches of the Yangtze River to precipitation and temperature
表4為長(zhǎng)江上游各水文站降水和氣溫對(duì)徑流變化的貢獻(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é)果偏大,在徑流減少的高場(chǎng)站、北碚站、寸灘站、武隆站和宜昌站結(jié)果偏小。3種方法得到的降水對(duì)徑流變化的貢獻(xiàn)量均值的絕對(duì)值最小的是宜昌站,為0.15 mm/a,最大的是直門達(dá)站,為1.58 mm/a;氣溫對(duì)徑流變化的貢獻(xiàn)量均值的絕對(duì)值最小的直門達(dá)站,為-0.30 mm/a,最大的時(shí)北碚站,為-1.17 mm/a。沿長(zhǎng)江干流從上到下各水文站降水對(duì)徑流變化的貢獻(xiàn)量整體上呈現(xiàn)減少趨勢(shì),氣溫對(duì)徑流變化的貢獻(xiàn)量呈現(xiàn)增加趨勢(shì)。對(duì)于降水和氣溫對(duì)徑流變化的貢獻(xiàn)量,定義絕對(duì)值大者為主導(dǎo)因素。可以發(fā)現(xiàn),直門達(dá)站、石鼓站、屏山站降水對(duì)徑流變化的貢獻(xiàn)量大于氣溫對(duì)徑流變化的貢獻(xiàn)量,而高場(chǎng)站、北碚站、寸灘站、武隆站以及宜昌站氣溫對(duì)徑流變化的貢獻(xiàn)量逐漸超過降水的貢獻(xiàn)量,成為徑流變化趨勢(shì)的主導(dǎo)因素。
表4 長(zhǎng)江上游各水文站降水和氣溫對(duì)徑流變化的貢獻(xiàn)量Table 4 Contribution of precipitation and temperature at hydrologic stations in upper reaches of the Yangtze River to runoff change
采用GA-BP算法進(jìn)行徑流對(duì)氣象要素變化的敏感性分析時(shí),輸入變量的微小擾動(dòng)都會(huì)影響敏感性系數(shù)的計(jì)算結(jié)果。為探究GA-BP算法計(jì)算徑流對(duì)氣象要素變化的敏感性結(jié)果的穩(wěn)定性,降水?dāng)_動(dòng)分別取測(cè)試集降水量均值的1%、2%、3%、5%和10%,氣溫?cái)_動(dòng)分別取測(cè)試集氣溫均值的1%、2%、3%、5%和10%,計(jì)算不同擾動(dòng)得到的徑流對(duì)氣象要素變化的敏感性,結(jié)果見表5??梢钥闯?,隨著擾動(dòng)程度的增加,GA-BP算法得到的徑流對(duì)降水與氣溫的敏感性系數(shù)逐漸變小,且變化幅度同樣變小并趨于穩(wěn)定。考慮到采用較小擾動(dòng)時(shí)計(jì)算誤差對(duì)于敏感性結(jié)果影響較大,并且氣象要素在1979—2015年的實(shí)際變化程度平均約3%,而且擾動(dòng)在3%時(shí)其機(jī)器學(xué)習(xí)敏感性結(jié)果的波動(dòng)幅度較小且趨于穩(wěn)定,因此在利用GA-BP算法計(jì)算敏感性時(shí),選擇3%的擾動(dòng)進(jìn)行分析。
表5 不同程度擾動(dòng)GA-BP算法敏感性計(jì)算結(jié)果Table 5 Sensitivity calculation results of GA-BP algorithm with different disturbances
為研究GA-BP算法利用不同時(shí)間尺度數(shù)據(jù)進(jìn)行敏感性分析的適用性,分別采用日尺度和月尺度數(shù)據(jù)進(jìn)行計(jì)算,結(jié)果見表6。可見,日尺度與月尺度計(jì)算結(jié)果中個(gè)別站點(diǎn)有所差異,但總體相差不大,計(jì)算結(jié)果穩(wěn)定,所有站點(diǎn)的敏感性系數(shù)均值近乎相等。這意味著GA-BP算法不僅可以利用日尺度數(shù)據(jù)計(jì)算徑流對(duì)降水與氣溫變化的敏感性,也可利用月尺度數(shù)據(jù)進(jìn)行計(jì)算,從而可以在缺少日數(shù)據(jù)的流域進(jìn)行應(yīng)用。
表6 不同時(shí)間尺度數(shù)據(jù)GA-BP算法敏感性計(jì)算結(jié)果Table 6 Sensitivity calculation results of GA-BP algorithm using data of different time scales
沿長(zhǎng)江干流從上向下,各水文站徑流對(duì)降水與氣溫的敏感性總體呈增加趨勢(shì),高場(chǎng)站徑流對(duì)氣溫的敏感性偏高,以下進(jìn)一步探究這種分布形成的原因。直門達(dá)站、石鼓站、屏山站、高場(chǎng)站、北碚站、寸灘站、武隆站、宜昌站的干旱指數(shù)(ET0/P)分別為1.75、1.60、1.28、0.82、0.92、1.08、0.73和1.00,可見,高場(chǎng)站控制的流域干旱指數(shù)明顯低于除武隆外的其他流域。圖3為長(zhǎng)江上游各水文站徑流敏感性與干旱指數(shù)的關(guān)系。如圖3(a)所示,徑流對(duì)降水的敏感性與干旱指數(shù)間存在很強(qiáng)的線性關(guān)系,敏感性系數(shù)隨著干旱指數(shù)的增大而減小;如圖3(b)所示,徑流對(duì)氣溫的敏感性隨著干旱指數(shù)的增大而減小。對(duì)于氣象要素的貢獻(xiàn)量,如圖3(c)所示,降水對(duì)徑流變化的貢獻(xiàn)量與干旱指數(shù)之間也存在顯著的線性關(guān)系,敏感性系數(shù)隨著干旱指數(shù)的增大而增大;如圖3(d)所示,氣溫對(duì)徑流變化的貢獻(xiàn)量隨著干旱指數(shù)的增大而減小。即降水對(duì)徑流變化的影響在越干旱的區(qū)域影響越大,而氣溫對(duì)徑流變化的影響在越濕潤(rùn)的區(qū)域越大,這與之前研究結(jié)果一致[37]。
(a) 徑流對(duì)降水的敏感性
(b) 徑流對(duì)氣溫的敏感性
(c) 降水對(duì)徑流變化貢獻(xiàn)量
(d) 氣溫對(duì)徑流變化貢獻(xiàn)量圖3 長(zhǎng)江上游各水文站徑流敏感性與干旱指數(shù)的關(guān)系Fig.3 Relationship between runoff sensitivity of hydrologic stations and drought index in upper reaches of the Yangtze River
GA-BP算法對(duì)徑流過程模擬效果較好,但對(duì)洪峰的模擬效果較差,NSE指標(biāo)表現(xiàn)一般。以往研究運(yùn)用機(jī)器學(xué)習(xí)法對(duì)長(zhǎng)江上游進(jìn)行徑流模擬時(shí),NSE值一般為0.6~0.8,如熊一橙等[38]利用LSTM算法模擬北碚站、高場(chǎng)站徑流的NSE值分別為0.77、0.62;Zhu等[39]利用ANN算法模擬北碚站、高場(chǎng)站日徑流NSE值為0.78、0.74;黃鈺瀚[40]利用VIC模型對(duì)寸灘站日徑流進(jìn)行模擬,NSE值為0.80。對(duì)比以往研究,本文GA-BP算法對(duì)長(zhǎng)江上游日徑流模擬的結(jié)果相差不大。除此之外,在計(jì)算年徑流量對(duì)年降水量與年均氣溫的敏感性系數(shù)時(shí),日徑流或洪峰的模擬效果不佳對(duì)于年徑流量的計(jì)算結(jié)果影響不大,因此可以認(rèn)為該方法盡管存在一定的不確定性,但得出的長(zhǎng)江上游流域徑流變化的敏感性總體是可信的。
本文計(jì)算得到的降水與氣溫對(duì)徑流變化的貢獻(xiàn)量,與實(shí)測(cè)值之間存在一定差異。導(dǎo)致這一差異的可能原因是只考慮了降水與氣溫對(duì)徑流變化的影響,對(duì)其他氣象要素、下墊面植被的變化以及人類活動(dòng)的影響并未加以考慮。從結(jié)果來看,徑流增加的直門達(dá)站和石鼓站,實(shí)測(cè)的增加值小于估計(jì)值;徑流減少的高場(chǎng)站、北碚站、寸灘站、武隆站和宜昌站,實(shí)測(cè)的減少值大于估計(jì)值,可能原因是人類活動(dòng)及植被變化導(dǎo)致耗水的增加、徑流減少。
a.基于GA-BP算法的敏感性分析,利用日徑流數(shù)據(jù)得到的結(jié)果整體較好,與多元線性回歸法以及Budyko框架法的結(jié)果總體一致;對(duì)于缺少日徑流數(shù)據(jù)的流域,該方法也可利用月數(shù)據(jù)進(jìn)行計(jì)算。
b.年徑流量方面,高場(chǎng)站顯著減少(p<0.01),趨勢(shì)為-2.81 mm/a,其余站點(diǎn)無顯著變化;年降水量方面,直門達(dá)站和石鼓站顯著增加(p<0.01),趨勢(shì)分別為4.28 mm/a和3.13 mm/a,其余站點(diǎn)無顯著變化;年氣溫方面,各站點(diǎn)均顯著上升,幅度為0.021~0.062℃/a。
c.沿長(zhǎng)江干流自上而下,徑流變化對(duì)降水的敏感性系數(shù)由直門達(dá)站的0.37 mm/mm逐漸增大至宜昌站的0.74 mm/mm,最高為高場(chǎng)站的0.76 mm/mm;而徑流變化對(duì)氣溫的敏感性也呈現(xiàn)增大趨勢(shì),范圍從-4.77~-33.53 mm/℃;敏感性與干旱指數(shù)相關(guān),隨著干旱指數(shù)的增加,徑流變化對(duì)降水和氣溫的敏感性減小。
d.從對(duì)徑流變化的貢獻(xiàn)量來看,降水的貢獻(xiàn)量沿干流自上而下逐漸變小,氣溫的貢獻(xiàn)量無明顯規(guī)律。其中直門達(dá)站、石鼓站、屏山站徑流的變化主要由降水主導(dǎo);向下至高場(chǎng)站、北碚站、寸灘站、武隆站、宜昌站,氣溫的貢獻(xiàn)量超過了降水,居于主導(dǎo)地位。呈現(xiàn)出隨著干旱指數(shù)增加降水對(duì)徑流變化的影響增大、氣溫對(duì)徑流變化的影響減小的規(guī)律。