李福威，包愛美，疏杏勝，丁 偉

（1.國電電力和禹水電開發(fā)公司，遼寧本溪 117201；2.大連理工大學(xué)水利工程學(xué)院，遼寧大連 116024）

0 引言

長期徑流預(yù)報對掌握未來徑流信息，實現(xiàn)水資源的高效利用、水電站優(yōu)化運行、防汛抗旱等具有重要意義。然而，由于長期預(yù)報的預(yù)見期較長，其受大氣環(huán)流、下墊面情況、人類活動等多重因素影響［1］，徑流過程的不確定性大，預(yù)報難度較大，預(yù)報精度難以保證。目前，對于降低較長預(yù)見期徑流預(yù)報的不確定性、提高預(yù)報精度的研究，主要考慮兩個方面，一是預(yù)報因子篩選方面，探究徑流形成的物理基礎(chǔ)和影響徑流變化的各因素及其間的相互作用，從眾多因子中篩選出與預(yù)報對象密切相關(guān)的因子；二是預(yù)報模型方面，研究預(yù)報模型的方法原理及其適用性，進行多模型的比較分析和綜合利用。

隨著中長期徑流預(yù)報研究的快速發(fā)展，可利用預(yù)報方法較多，主要分為兩類，一類是根據(jù)徑流序列自身規(guī)律，僅考慮徑流單要素作為預(yù)報因子的歷史演變法、時間序列法等傳統(tǒng)統(tǒng)計方法［2-4］；一類是挖掘水文和氣象數(shù)據(jù)的潛在規(guī)律，考慮水文、氣象等多要素作為預(yù)報因子的回歸分析法和機器學(xué)習(xí)法。后者可選用一切與預(yù)報對象存在潛在關(guān)系的影響因素［5］，包括考慮入海平面氣壓、大氣環(huán)流因子等氣象因子用作預(yù)報因子來實現(xiàn)對流域徑流的預(yù)報，具有較強的物理基礎(chǔ)，獲得較長的預(yù)見期和相對較高的預(yù)報精度［6-8］。朱春苗等［9］利用相關(guān)系數(shù)法、主成分分析法和互信息法優(yōu)選松花江流域上游3個水文代表站的預(yù)報因子，采用SVR（Support Vector Regression）模型進行流域月徑流預(yù)報，得出適用于該流域月徑流預(yù)報的最優(yōu)預(yù)報因子組合模型。雷莉和王超［10］從130項氣候系統(tǒng)指數(shù)中篩選出預(yù)報因子后，基于篩選的因子構(gòu)建了石羊河流域BP-ANN（Back Propaga‐tion-Artificial Neural Network）、Elman 和PSO-SVR（Particle Swarm Optimization-SVR）三種年徑流預(yù)報模型并比較各模型預(yù)報結(jié)果，確定Elman 和PSO-SVR 中長期徑流預(yù)報模型可為該流域中長期徑流預(yù)報提供支撐。李伶杰等［11］以龍江水庫入庫徑流預(yù)報為研究對象，從環(huán)流指數(shù)、海溫、氣壓和前期月徑流中選取關(guān)鍵預(yù)報因子，建立隨機森林與支持向量機模型，發(fā)現(xiàn)太平洋中北部與西部氣候因子對徑流預(yù)報的影響較大。

考慮到長期徑流預(yù)報受影響因素多，預(yù)報不確定性較大的問題，本研究以渾江桓仁水庫年入庫徑流和汛期月入庫徑流預(yù)報為研究對象，考慮將前期降雨、徑流作為預(yù)報因子的同時引入大氣環(huán)流因子，建立包括統(tǒng)計分析法和機器學(xué)習(xí)法的多個長期徑流預(yù)報模型，比較各模型在桓仁流域的適用性，并分析各模型對年徑流預(yù)報及汛期各月徑流預(yù)報的預(yù)報水平，給出最優(yōu)預(yù)報方案，提高桓仁流域長期徑流預(yù)報的精度。

1 材料與方法

1.1 研究區(qū)概況及數(shù)據(jù)資料

桓仁水庫作為渾江梯級電站的龍頭水庫，總庫容為36.4 億m3，電站裝機容量為222.5 MW，兼有防洪、灌溉等綜合作用?；溉仕畮煲陨狭饔蛴兄行∷畮旃?7 座，總庫容合計2.43 億m3，其中中型水庫6 座，總庫容1.96 億m3，占比80.7%，小型水庫41座，庫容0.48 億m3，占比19.3%。流域多年平均年降水量860 mm，降水年內(nèi)分配不均，70%的雨量集中在6-9 月間，且在7、8月達到最大，大洪水主要發(fā)生在7 月下旬至8 月中旬；冬季一般從11 月份開始到翌年3 月末或4 月初結(jié)束，期間主要以降雪為主，受積雪影響，冬季徑流一般較枯。因此，桓仁水庫來水從5月份流域發(fā)生春汛開始一直持續(xù)到10月份，汛期作為桓仁水庫年內(nèi)主要來水階段，充分把握汛期各月來水情況，對于電站做出全局較優(yōu)的發(fā)電調(diào)度過程，提升發(fā)電效益具有重要意義；而年徑流預(yù)報可為水庫運行管理提供更長遠的參考信息，幫助管理者更早地制定調(diào)度計劃，從而更高效地利用水力資源。

徑流資料來源于桓仁水庫1967-2017 年月實測入庫流量，1967-2017 年74 項大氣環(huán)流因子的逐月數(shù)據(jù)源自國家氣候中心網(wǎng)站（http：∕∕www.ncc-cma.net∕cn∕），桓仁水庫流域1967-2017年逐月實測降雨量源自中國氣象數(shù)據(jù)網(wǎng)（http：∕∕data.cma.cn∕）。

1.2 研究方法

1.2.1 預(yù)報因子篩選

（1）因子初選。采用相關(guān)系數(shù)法對預(yù)報因子進行初選，按照相關(guān)性大小分別從諸多水文、氣象因子選取預(yù)報因子。由于預(yù)報因子對徑流的影響存在滯后效應(yīng)，預(yù)報年徑流時，水文因子取桓仁水庫預(yù)報年份前3 年的實測年降雨量與年徑流量，氣象因子取預(yù)報時刻前一年內(nèi)國家氣候中心各月份發(fā)布的74 項大氣環(huán)流因子。預(yù)報月徑流時，前12 個月的水文、氣象因子作為輸入因子。相關(guān)性系數(shù)計算公式為：

式中：r為相關(guān)系數(shù)；n為資料樣本數(shù)；Xi與分別為輸入因子及其序列均值；Yi與分別為預(yù)報對象及其序列均值。r的取值范圍為［-1，1］，其正負號表示兩變量之間是正相關(guān)還是負相關(guān)，絕對值越接近1表明兩變量的線性相關(guān)程度越高。

（2）因子精選。在使用相關(guān)系數(shù)法初選出因子后，進一步結(jié)合預(yù)報模型分析不同因子組合對模擬和檢驗結(jié)果的影響，以確定最終的預(yù)報因子。本文采用向前搜索包裹法，以率定期殘差平方和最小為目標(biāo)，通過評價預(yù)報模型在不同因子組合情況下的預(yù)報結(jié)果，確定最優(yōu)預(yù)報因子組合，具體步驟如下：

①將初選因子相關(guān)性系數(shù)排序第一的因子x0作為固定因子，分別與其他因子組合輸入到預(yù)報模型，選擇率定期的殘差平方和f(x0，xi)最小時加入的因子xi作為選定因子。

②以x0、xi為固定因子，分別與其他因子組合輸入到預(yù)報模型，確定使f(x0，xi，xj)最小時的因子xj，若f(x0，xi，xj)小 于f(x0，xi)，則將xj因子加入模型，進行步驟③；否則，排除該因子。

③重復(fù)步驟②，不斷增加輸入因子個數(shù)，直至f不再減小或者f的遞減速度小于一定程度時，從而確定出最終預(yù)報因子集。

1.2.2 預(yù)報模型構(gòu)建

選用多元線性回歸（Multiple Linear Regression，MLR）［12］、多元門限回歸（TR）［13］、逐步回歸（Stepwise Regression，SR）［14］、人工神經(jīng)網(wǎng)絡(luò)（Back Propagation-Artificial Neural Network，BPANN）［15］、支持向量機（Support Vector Machine，SVM）［16，17］和基于主成分分析的人工神經(jīng)網(wǎng)絡(luò)模型（Artificial Neural Network Model based on Principal Component Analysis，PCA-BP-ANN）六個模型，以最終篩選的水文、氣象因子為模型輸入，以預(yù)報徑流的合格率（QR）為評定指標(biāo)，分別構(gòu)建桓仁水庫流域年徑流預(yù)報模型和汛期月徑流預(yù)報模型。據(jù)《水文情報預(yù)報規(guī)范GB∕T 22482-2008》，許可誤差限為多年同期變幅的20%，桓仁水庫年徑流多年變幅的20%為14.9 億m3，汛期各月的允許誤差（已轉(zhuǎn)化為流量）見表1。

表1 汛期各月允許預(yù)報誤差 m3∕sTab.1 Allowable forecast error of each month in the flood season

多元線性回歸（MLR）、多元門限回歸（TR）、逐步回歸（SR）等傳統(tǒng)統(tǒng)計模型在中長期徑流預(yù)報中應(yīng)用較早，其結(jié)構(gòu)簡單、易于操作，只需確定各預(yù)報模型數(shù)學(xué)方程，即可得到較好的預(yù)報效果。人工神經(jīng)網(wǎng)絡(luò)（BP-ANN）和支持向量機（SVM）作為在中長期徑流預(yù)報中應(yīng)用最廣泛且較成熟的機器學(xué)習(xí)模型，具有較強的非線性映射能力，相對傳統(tǒng)的統(tǒng)計方法，能夠更好地刻畫徑流過程非線性及非穩(wěn)態(tài)性的特點。選用具有3層網(wǎng)絡(luò)結(jié)構(gòu)（輸入層、隱含層和輸出層）的人工神經(jīng)網(wǎng)絡(luò)（BP-ANN）模型，以年徑流預(yù)報合格率為指標(biāo)，采用試錯法確定網(wǎng)絡(luò)隱含層，并經(jīng)反復(fù)訓(xùn)練確定每層結(jié)構(gòu)所含神經(jīng)元的個數(shù)；SVM 模型能夠快速處理小樣本的非線性問題，具有較強的泛化能力，模型需確定的參數(shù)有懲罰函數(shù)C、不敏感損失系數(shù)ε、核函數(shù)及其所含參數(shù)g，本文采用試錯法確定核函數(shù)，并采用網(wǎng)格搜索法優(yōu)化確定懲罰系數(shù)C和核函數(shù)參數(shù)g。

主成分分析（Principal Component Analysis，PCA）法是通過解析原來具有一定相關(guān)性的多變量，將原變量進行線性組合導(dǎo)出一組新的綜合變量，使這些新的綜合變量能夠盡可能多地反映原始變量的信息，以達到簡化數(shù)據(jù)和降維的目的［18］。本文構(gòu)建基于主成分分析的人工神經(jīng)網(wǎng)絡(luò)模型（PCA-BP-ANN）時，采用PCA 法將篩選后的預(yù)報因子重新組合為新的變量，并根據(jù)計算的特征值、方差貢獻率和累積貢獻率確定主成分個數(shù)，以新確定的主成分作為輸入因子輸入到BP-ANN 模型進行徑流預(yù)測。

2 結(jié)果與分析

2.1 預(yù)報因子篩選

（1）因子初選?？紤]水文和氣象兩類因子對徑流變化的影響，采用相關(guān)性分析法分別對流域水文、氣象因子進行初選，取相關(guān)性系數(shù)較大的因子作為備選因子。對桓仁水庫流域月徑流量與前期各月份實測降雨、徑流和大氣環(huán)流因子進行相關(guān)性分析，得到該流域長期徑流預(yù)報的備選因子，結(jié)果見表2 和表3。由表可知，對于月徑流，除了上月月降雨量相關(guān)性系數(shù)達0.57，其他水文因子與預(yù)報對象的相關(guān)性都低于備選的大氣環(huán)流因子；對于預(yù)見期更長的年徑流，雖然大氣環(huán)流因子與年徑流的相關(guān)性系數(shù)相對月徑流有所減小，但其與預(yù)報對象的相關(guān)性依舊遠高于水文要素。

表2 相關(guān)性較高的水文因子Tab.2 Hydrological factors with high correlation

（2）因子精選。結(jié)合所選預(yù)報模型，采用向前搜索包裹法，對表2 和表3 確定的備選因子進一步篩選，分別確定各預(yù)報模型的最優(yōu)組合因子作為各模型的最終輸入因子，表4 列出了結(jié)合多元線性回歸預(yù)報模型確定的最終預(yù)報因子，以此作為多元線性回歸預(yù)報模型的輸入，同理可確定其他預(yù)報模型的最終預(yù)報因子。

表3 相關(guān)性較高的大氣環(huán)流因子Tab.3 Atmospheric circulation factors with high correlation

采用主成分分析（PCA）法將表4 中選定的預(yù)報因子進行線性組合，使其成為一組新的線性無關(guān)的綜合變量，并根據(jù)特征值的累積貢獻率達90%確定主成分個數(shù)，各成分的特征值及其方差貢獻率和累積貢獻率結(jié)果見表5。由表5 可知，對于汛期月徑流預(yù)報，前3 個主成分的累積貢獻率為90.5%，而后3 個主成分的貢獻率較低，對汛期月徑流預(yù)報的影響不大，因此選定前3 個主成分為模型輸入，式（2）～（4）分別給出了3 個主成分與標(biāo)準(zhǔn)化原始變量的關(guān)系；對于年徑流預(yù)報，前5個主成分的累積貢獻率才能達90%以上，即需要5 個主成分方可有效地代表原始變量90%以上的信息。此外，年徑流預(yù)報確定的主成分個數(shù)比月徑流預(yù)報多，說明預(yù)報預(yù)見期更長的年徑流所受影響因素要比預(yù)報月徑流更為復(fù)雜，其預(yù)報不確定性更大。

表4 多元線性回歸預(yù)報模型最終預(yù)報因子Tab.4 Determined predictors

表5 應(yīng)用PCA法計算的各成分特征值、方差貢獻率和累積貢獻率%Tab.5 The eigenvalues，variance contribution rate and cumulative contribution rate of each component calculated by PCA method

式中：Fi（i=1，2，3）表示主成分；Vj（j=1，…，6）表示標(biāo)準(zhǔn)化原始變量。

2.2 模型構(gòu)建及預(yù)報效果分析

將桓仁水庫流域?qū)崪y徑流資料以及選中的預(yù)報因子劃分為率定期和驗證期兩部分，率定期為1967-2000 年，驗證期為2001-2017 年，以預(yù)報徑流的合格率（QR）為評定指標(biāo)，分別構(gòu)建多個桓仁水庫流域長期徑流預(yù)報模型并驗證。式（5）和式（6）分別給出了流域汛期月徑流和年徑流的多元線性回歸模型（MLR）預(yù)報方程，而多元門限回歸（TR）和逐步回歸（SR）方法類似；表6 展示了人工神經(jīng)網(wǎng)絡(luò)（BP-ANN）和支持向量機（SVM）兩機器學(xué)習(xí)模型的參數(shù)結(jié)果；采用主成分分析（PCA）結(jié)果，將月預(yù)報確定的3 個主成分和年預(yù)報確定的5 個主成分分別輸入到BP-ANN模型中，構(gòu)建月徑流預(yù)報PCA-BP-ANN 模型和年徑流預(yù)報PCA-BP-ANN模型。

表6 BP-ANN和SVM模型參數(shù)Tab.6 The parameters of BP-ANN and SVM models

式（5）和式（6）中xi分別為該模型最終確定的月徑流預(yù)報輸入因子和年徑流預(yù)報輸入因子（見表4）。

圖1給出了各月徑流預(yù)報模型在率定期和驗證期的模擬預(yù)報精度，由圖可知，無論是在率定期還是在驗證期，不同模型對不同月份的模擬預(yù)報精度均有一定差異。在率定期，TR模型在8、9 月份的模擬預(yù)報效果最好，合格率達90%，而在5、6、7 月份的表現(xiàn)卻不如SVM模型；雖然SVM和BP-ANN機器學(xué)習(xí)模型的合格率在5、6 月份的合格率明顯優(yōu)于傳統(tǒng)回歸模型，但在8、9、10 月份卻不如回歸模型；結(jié)合主成分分析法的PCA-BP-ANN模型在6月份表現(xiàn)最差，但在7月和9月份的合格率要比沒有采用PCA 方法的BP-ANN 模型高10%左右。在驗證期，TR 模型在5月和10月份預(yù)報效果最差，預(yù)報合格率不足60%，而在7、8月份卻能達85%以上，其中7 月份合格率達94%；PCA-BPANN模型雖在6、7月份預(yù)報合格率不足50%，但在5月份和8月份表現(xiàn)最佳，合格率均在85%以上，且該模型相對BP-ANN 模型精度提升了10%左右，說明結(jié)合主成分分析法（PCA）進行該流域月徑流預(yù)報可有效提高汛期某些月份的預(yù)報精度。對比驗證期和率定期，各模型在驗證期的預(yù)報精度相對率定期均有所降低，而SVM 和BP-ANN 等機器學(xué)習(xí)模型在汛期各月表現(xiàn)相對穩(wěn)定。

圖1 汛期月預(yù)報各模型在率定期和驗證期的模擬精度Fig.1 The simulation accuracy of monthly forecast models in flood season for both calibration and verification periods

由此可見，各模型均有“擅長”預(yù)報的月份，為充分發(fā)揮模型各自價值，實現(xiàn)多模型間的相互協(xié)調(diào)和性能互補，綜合各模型的預(yù)報效果，表7 給出了汛期各月份推薦使用模型的最優(yōu)預(yù)報方案。在預(yù)報各月來水時，采用該月份對應(yīng)的推薦模型，除6月和10 月份，其他月份合格率均能在70%以上，并且在汛期來水最重要的7、8月份預(yù)報合格率能達90%以上。

表7 汛期各月份推薦使用模型Tab.7 Recommended model for each month in the flood season

從各模型在不同月份預(yù)報的整體效果來看，6 月預(yù)報效果較差，主要是受流域初始狀態(tài)的影響，6月初流域內(nèi)中小水庫基本呈庫空狀態(tài)，同時該時段為流域灌溉期，因而該月徑流受流域蓄水影響較大。7、8月份預(yù)報精度較高，是因為該流域在7、8月份降水量較大，徑流量也大，從而受其他影響因素（包括流域初始狀態(tài)、水利工程等因素）較小，允許誤差大（見表1），因此預(yù)報合格率較高，預(yù)報效果較好。而9 月和10 月份預(yù)報精度較8月份降低，是因該流域于9月份開始降水逐漸減少，而上游眾多水庫會在9 月和10 月份蓄水截留導(dǎo)致預(yù)報水平降低。為進一步分析各月份預(yù)報效果的差異，計算汛期各月徑流的變差系數(shù)（Cv）分析徑流的不確定性，結(jié)果見圖2，Cv值越大表示該月來水不確定性越大。由圖2 知，5 月份Cv值最低，考慮是5 月徑流大多由流域退水產(chǎn)生，來水相對穩(wěn)定；而在6 月進入汛期后，徑流受天然降水影響較大，不確定性較大，因此5月份預(yù)報效果優(yōu)于6 月份；而10 月份Cv值最大，反應(yīng)該月受降水和水利工程蓄水等影響導(dǎo)致來水不確定性最大，因而預(yù)報效果較差。

圖2 汛期各月實測徑流不確定性分析Fig.2 Uncertainty analysis of monthly measured runoff in flood season

圖3為年徑流預(yù)報各模型精度，對于年徑流預(yù)報，各模型在率定期的預(yù)報合格率除了TR 模型外均達到90%以上，取得了較高的預(yù)報精度。驗證期的預(yù)報精度有所下降，除了SR 模型外，其他模型預(yù)報合格率都在80%以上，但BP-ANN 和SVM 等機器學(xué)習(xí)模型預(yù)報精度明顯高于傳統(tǒng)的統(tǒng)計模型，其中SVM 模型的合格率為88%，BP-ANN 和PCA-BP-ANN 模型的合格率在90%以上，比統(tǒng)計模型高10%左右。此外，由于年徑流預(yù)報的預(yù)見期較長，不確定性大，為了防止個別模型的預(yù)報誤差較大而誤導(dǎo)決策，推薦選用BP-ANN 和PCA-BP-ANN 模型，同時綜合考慮其他模型的預(yù)報結(jié)果，最終確定合適的年徑流預(yù)報結(jié)果。

圖3 年徑流預(yù)報各模型在率定期和驗證期的模擬精度Fig.3 The simulation accuracy of annual runoff forecast models for both calibration and verification periods

3 結(jié)論

耦合相關(guān)性分析方法和向前搜索包裹法確定桓仁流域汛期各月徑流與年徑流預(yù)報模型輸入，并基于主成分分析重組因子以簡化模型輸入，綜合考慮統(tǒng)計方法和機器學(xué)習(xí)方法，分別建立年徑流預(yù)報模型和汛期月徑流預(yù)報模型，對比分析各模型的預(yù)報水平，給出最佳預(yù)報方案，以提高桓仁流域長期徑流預(yù)報的精度。主要結(jié)論如下：

（1）在月徑流預(yù)報中，篩選出6 個關(guān)鍵影響因子，且上個月月降雨量起主要作用，其余6 個大氣環(huán)流因子與預(yù)報對象的相關(guān)性和上月月降雨量水平相當(dāng)，相關(guān)性系數(shù)均在0.5 以上；年徑流預(yù)報中，篩選出7個關(guān)鍵影響因素，均為大氣環(huán)流因子。

（2）基于多元線性回歸（MLR）、多元門限回歸（TR）、逐步回歸（SR）、人工神經(jīng)網(wǎng)絡(luò)（BP-ANN）、支持向量機（SVM）和結(jié)合主成分分析的人工神經(jīng)網(wǎng)絡(luò)（PCA-BP-ANN）六種方法，構(gòu)建了桓仁流域年徑流預(yù)報模型和汛期月徑流預(yù)報模型。對于年徑流預(yù)報，BP-ANN、SVM 和PCA-BP-ANN 三個機器學(xué)習(xí)模型的預(yù)報效果優(yōu)于傳統(tǒng)統(tǒng)計模型，其合格率要比統(tǒng)計模型高10%左右；對于汛期月徑流預(yù)報，各模型均有自己“擅長”預(yù)報的月份，利用最佳預(yù)報方案預(yù)報時，在汛期來水最重要的7、8 月份預(yù)報合格率均能達90%以上。

在進行流域長期徑流預(yù)報時，考慮徑流形成的物理基礎(chǔ)和影響徑流變化的各因素選擇合適的預(yù)報因子，同時綜合多種預(yù)報方法擇優(yōu)選擇最佳預(yù)報方案，可有效降低流域長期徑流預(yù)報的不確定性。此外，研究采用預(yù)見期為一個月，可滿足多模型綜合預(yù)報方法在該研究區(qū)有效性的論證，而該方法在不同預(yù)見期的有效性可能會有所區(qū)別，還需結(jié)合實際調(diào)度需求進一步檢驗，探索不同預(yù)見期下的最優(yōu)預(yù)報方案，以更加合理地指導(dǎo)水庫調(diào)度。