長江上游徑流對氣象要素變化的敏感性分析

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

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

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

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

1 研究區概況與研究數據

1.1 研究區概況

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

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

1.2 研究數據

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

ET0=-0.611+0.149RS+0.079T

(1)

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

2 研究方法

2.1 GA-BP算法

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

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

2.2 多元線性回歸法

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

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

(2)

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

2.3 Budyko框架法

Budyko框架法假設流域實際蒸散發量可以表示為流域降水量和潛在蒸散發量的函數[34]，據此Choudhury等[35-36]提出了Choudhury-Yang公式：

(3)

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

R=P-ET

(4)

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

(5)

3 結果與分析

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

表1為長江上游各水文站1979—2015年徑流、降水量及氣溫的變化趨勢?？梢?，徑流增加趨勢最大的為直門達站，達0.69 mm/a；減少趨勢最大的為高場站，為2.81 mm/a。沿長江干流從上向下，徑流變化由直門達站、石鼓站、屏山站的增加趨勢，逐漸轉變為高場站、北碚站、寸灘站、武隆站和宜昌站的下降趨勢；值得注意的是，在p=0.01的顯著性下，除高場站外其余站點變化趨勢均不顯著。降水量方面，直門達站和石鼓站控制流域顯著增加(p<0.01)，分別為4.28 mm/a和3.13 mm/a；屏山、寸灘和宜昌站點控制流域呈不顯著增加，其余站點控制流域呈不顯著減少。氣溫方面，各站點控制流域內氣溫均呈現顯著上升趨勢(p<0.01)，變化范圍從武隆站的0.021℃/a到直門達站的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算法徑流模擬結果。8個水文站訓練集的NSE均值為0.71，R2均值為0.71，RMSE均值為0.59 mm/d；測試集NSE均值為0.61，R2均值為0.64，RMSE均值為0.55 mm/d，整體模擬效果良好。直門達站、石鼓站以及宜昌站模擬效果最好，訓練集NSE值均高于0.7，測試集NSE值均高于0.6；北碚站和武隆站的模擬效果相對較差，主要是因為日徑流數據波動較大，GA-BP算法對洪峰模擬效果較差，從而拉低整體模擬效果。

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

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

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

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

3.4 氣象要素對徑流變化的貢獻量

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

表4 長江上游各水文站降水和氣溫對徑流變化的貢獻量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算法進行徑流對氣象要素變化的敏感性分析時，輸入變量的微小擾動都會影響敏感性系數的計算結果。為探究GA-BP算法計算徑流對氣象要素變化的敏感性結果的穩定性，降水擾動分別取測試集降水量均值的1%、2%、3%、5%和10%，氣溫擾動分別取測試集氣溫均值的1%、2%、3%、5%和10%，計算不同擾動得到的徑流對氣象要素變化的敏感性，結果見表5?？梢钥闯?，隨著擾動程度的增加，GA-BP算法得到的徑流對降水與氣溫的敏感性系數逐漸變小，且變化幅度同樣變小并趨于穩定?？紤]到采用較小擾動時計算誤差對于敏感性結果影響較大，并且氣象要素在1979—2015年的實際變化程度平均約3%，而且擾動在3%時其機器學習敏感性結果的波動幅度較小且趨于穩定，因此在利用GA-BP算法計算敏感性時，選擇3%的擾動進行分析。

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

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

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

4.2 長江上游徑流敏感性與干旱指數的關系

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

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

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

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

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

4.3 不確定性分析

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

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

5 結 論

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

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

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

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