劉建華，徐文馨，石昕顏，陳 杰，胡召根，陳 華

（1. 天生橋一級水電開發(fā)有限責任公司水力發(fā)電廠，廣東 廣州 510600； 2. 水資源與水電工程科學國家重點實驗室 武漢大學，湖北 武漢 430072； 3. 海綿城市建設水系統(tǒng)科學湖北省重點實驗室（武漢大學），湖北 武漢 430072）

0 引 言

徑流預報中通常將預見期在一個月以上、一年以內(nèi)的預報稱為長期預報。有效的長期月徑流預報對于水庫的綜合利用有著十分重要的意義，特別是對于以發(fā)電為主的水庫來說，準確預報入庫流量是編制年度以及各季度發(fā)電計劃的重要基礎［1，2］。

大部分水庫都建在山區(qū)，氣象水文條件復雜，雨量站和水文站建設條件較差，缺乏高質量水文氣象實測資料，給長期徑流預報模型的構建帶來了較大的難度［3］，其中，氣象資料不足是限制以水文模型為基礎的過程驅動長期徑流預報模型使用的重要原因。因此，數(shù)理統(tǒng)計方法，以多元線性回歸為代表，因模型簡單，實現(xiàn)方便，對數(shù)據(jù)變化的適應性高，可以在實際工作中有效地指導生產(chǎn)實踐，應用最為廣泛［4-7］?，F(xiàn)有研究多在分析徑流與前期多因子之間的統(tǒng)計相關關系后，通過構建變因子多元線性回歸模型實現(xiàn)對不同預見期的徑流預報。

一般而言，長期的水文動態(tài)需要考慮水文大循環(huán)，即水在陸地、海洋和大氣3 種介質中的相互轉化。河川徑流的形成和季節(jié)性變化受到海溫和大氣環(huán)流因子的共同作用，故從水文循環(huán)的機理出發(fā)，綜合分析影響區(qū)域水文情勢的物理因素，是長期徑流預報應遵循的基本原則［1］。近代的大氣環(huán)流實測資料無論在種類還是時空尺度上都十分豐富，且獲取方便，已被廣泛用作長期徑流預報因子［8，9］。海洋對大氣變化具有獨特的“記憶功能”和“低通濾波”作用，因此海洋表面溫度場一直被認為是影響大氣環(huán)流和引起氣候異常的重要因素，其也已被引入預報模型指示徑流變化規(guī)律［10-14］。然而現(xiàn)有研究多采取固定海域的海洋表面溫度場［15-19］，包括厄爾尼諾活動區(qū)域和對我國水文情勢影響明顯的北太平洋地區(qū)格點海溫數(shù)據(jù)，忽略了海域溫度場分布對區(qū)域徑流的影響。鑒于此，本研究將考慮海溫場的“偶極型”空間分布特征，探究不同海域溫度場的關聯(lián)性對長期徑流預報的指示價值。

研究以南盤江流域天生橋一級（天一）水庫上游為研究對象，綜合考慮水庫前期徑流、大氣環(huán)流因子和海溫偶極因子作為多元線性回歸模型的關鍵因子，在年內(nèi)各月分別構建預見期為1～12個月的徑流預報模型，并探討因子組合方式對預報效果的影響，以期為其他流域的長期月徑流預報提供參考。

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

1.1 研究區(qū)概況

天生橋一級（天一）水庫位于珠江流域西江水系上游的南盤江流域，其壩址以上集水面積為50 139 km2。南盤江流域屬于亞熱帶季風氣候區(qū)，干濕季節(jié)變化明顯。5 月至10 月（夏半年）受西南季風影響，溫和多雨，且降水集中；11 月至次年4 月（冬半年）受大陸氣團影響，干暖少雨，常出現(xiàn)干季。天一水庫為不完全多年調(diào)節(jié)水庫，多年平均年徑流量為193 億m3，總庫容為102.6 億m3，以發(fā)電為主要任務，是國家“西電東送”的龍頭水庫，在“西部大開發(fā)”中起著重要作用。然而，流域缺乏可靠的長期歷史氣象數(shù)據(jù)，降雨資料缺、漏測時間長且質量不高，為天一水庫構建預報精度滿足電廠需求的長期徑流預報模型面臨著很大的挑戰(zhàn)。天一水庫上游流域如圖1所示。

圖1 天一水庫流域圖Fig.1 Location and characteristics of the Tianyi reservoir

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

（1）月徑流數(shù)據(jù)由天一電廠提供，選用時段為1953 年1 月至2020年12月。

（2）大氣環(huán)流因子下載自中國氣象局國家氣候中心網(wǎng)站（https：//www.ncc-cma.net/），每月初更新上月實測數(shù)據(jù)。網(wǎng)站共提供88項大氣環(huán)流指數(shù)，若指數(shù)缺測月數(shù)大于選用數(shù)據(jù)序列長度10%，則該指數(shù)棄用；其余有缺測項的指數(shù)，缺測項用其余部分均值插補，最終有75項大氣環(huán)流指數(shù)被保留，如表1所示，用于后續(xù)模型構建，選用時段為1953年1月至2020年12月。

表1 選用的75項大氣環(huán)流指數(shù)Tab.1 Information of selected 75 atmospheric indices.

（3）海溫數(shù)據(jù)選用Kaplan Extended SST V2海溫波動月值數(shù)據(jù)資料，下載自美國海洋與大氣管理局網(wǎng)站（http：//www.esrl.noaa.gov/psd/），該數(shù)據(jù)集是目前公認的可靠性較高的海溫數(shù)據(jù)集之一［20，21］，空間分辨率為5°×5°，選用的時間跨度為1953 年1 月至2020年12月。

2 研究方法

研究方法主要由①待選預報因子集生成；②預報因子篩選；和③多元線性回歸模型構建3 部分組成，以下對3 部分作簡要介紹。①由于大氣環(huán)流和海溫因子對徑流的影響往往在幾個月甚至更長的滯后時間之后才能反映出來，因此在長期徑流預報中，需要考慮各因子對徑流的影響滯后期。本研究對所有預見期均考慮預報發(fā)起點前推12個月的歷史徑流、大氣環(huán)流和海溫偶極作為待選預報因子，通過將不同的預報因子組合共生成7種待選預報因子集（即預報因子只考慮前期徑流、只考慮大氣環(huán)流、只考慮海溫偶極、考慮徑流和大氣環(huán)流、考慮徑流和海溫偶極、考慮大氣環(huán)流和海溫偶極、考慮徑流、大氣環(huán)流和海溫偶極）。②以待預報月徑流為因變量，以各預報因子為自變量，通過對自變量與因變量的相關性進行顯著性檢驗完成預報因子初選后，再通過逐步回歸法進一步篩選預報因子。③多元線性回歸模型構建過程中選用奇數(shù)年率定偶數(shù)年驗證的方法，對年內(nèi)各月各預見期分別構建7種因子組合的預報模型。下文對各部分作詳細介紹。

2.1 海溫偶極模型

研究表明，不同海域的表面海溫場（SST）存在著一定的空間關聯(lián)性，利用SST 的偶極特征構建的統(tǒng)計預報模型對降水和徑流的指示作用甚至優(yōu)于傳統(tǒng)的ENSO 事件［22］。SAJI 等［23］研究發(fā)現(xiàn)海溫場的空間分布具有偶極的特征，基于此，美國佐治亞理工學院的GEORGAKAKOS 教授團隊［22］提出了一種海溫偶極模型，通過在全球海域篩選與區(qū)域降水相關的海溫場以構建偶極模型預報未來降水，其中，海溫偶極的定義是特定大小和地理位置上兩個海洋區(qū)域的平均表面溫度異常的函數(shù)，表示成兩個海區(qū)表面溫度場的差或和。關于海溫偶極的篩選和模型構建的詳細過程可參考CHEN［22］和QIAN［24］等人的文章，以下就本研究對海溫偶極模型的使用作簡要說明。主要步驟為：①以Gerrity Skill Score（GSS）評分方法［25］為標準，從海溫網(wǎng)格大數(shù)據(jù)中識別顯著的海溫偶極，將其配對、擴展、篩選后，得到海溫偶極的最終形態(tài)；②重復步驟①后得到多對海溫偶極，為年內(nèi)各月分別挑選與徑流序列相關性最為顯著的海溫偶極；③以所有顯著海溫偶極區(qū)域海溫差為預報因子，以徑流為預報對象，對年內(nèi)各月分別構建線性回歸模型，用留一交叉驗證方法保留平均絕對誤差最小的10個顯著海溫偶極區(qū)域海溫差（即海溫偶極因子）；④對預報發(fā)起點前推12 個月的海溫場重復步驟①至③，獲得12×10 個海溫偶極，對每個影響滯時均保留10 個海溫偶極因子中與奇數(shù)年實測徑流序列皮爾遜相關系數(shù)最大的1個，作為長期徑流預報模型的備選預報因子。

2.2 預報因子篩選

2.2.1 相關性分析

考慮到本研究的預報因子集可能過大，故先分別計算各因子和徑流序列的皮爾遜相關系數(shù)，篩選出集合內(nèi)與徑流序列顯著相關的因子。在本研究中，率定期（奇數(shù)年）年內(nèi)各月樣本數(shù)為33，顯著性水平取0.05，則因子與徑流序列相關系數(shù)大于0.344者通過相關性挑選。

2.2.2 逐步回歸分析

經(jīng)2.2.1 步驟篩選出的多個因子可能存在兩兩之間顯著相關的情況，這會導致多元回歸出現(xiàn)多重共線性的問題，從而影響預報效果，因此需要采用逐步回歸法進一步篩選因子。逐步回歸法是按自變量和因變量的顯著性大小逐個引入變量，所以逐步回歸會產(chǎn)生與引入步驟數(shù)量同樣多的因子組合，對應多個多元線性回歸模型。為得到因子數(shù)量合理、涵蓋主要信息且率定期擬合效果合格的模型，在逐步回歸分析中以率定期的修正復相關系數(shù)為標準確定回歸的最終結果。修正的復相關系數(shù)矯正了復相關系數(shù)擬合優(yōu)度對自由度的依賴關系，使得只有當自變量確實對因變量有所作用時值才會增加。在本研究中，率定期修正復相關系數(shù)達到0.8 時，對應的最少因子個數(shù)組合即為逐步回歸最終結果；若所有因子組合在率定期擬合效果都無法達到修正復相關系數(shù)0.8，則選擇達到0.75 對應的最少因子個數(shù)組合；達不到0.75 則選用0.65；率定期修正復相關系數(shù)小于0.65判定為無法有效擬合。

2.3 多元線性回歸模型

多元回歸模型能夠比較全面地綜合各個預報因子的作用。本研究將由逐步回歸法挑選出的預報因子按最小二乘原理確定各因子的回歸系數(shù)，建立多元線性回歸方程，公式為：

式中：x（t）i為預報因子序列；i為預報因子序號；n為預報因子總個數(shù)；bc，i為回歸系數(shù)；fc為常數(shù)項。

多元線性回歸模型中的回歸系數(shù)和常數(shù)項采用最小二乘法進行估計，最小二乘法的目標是使得預報結果與觀測值之間的離差平方和最小，離差平方和的計算公式如下所示：

式中：m為樣本數(shù)；j為樣本點序號。

2.4 預報效果評價指標

長期徑流預報的精度的計算公式參考《水文情報預報規(guī)范》（GB/T22482-2008）［26］，具體如下：

式中：Qoi和Qsi分別代表實測和預報月流量值。其中，預報值大于實際值2倍以上時，預報精度按0處理。

3 研究結果與討論

圖2展示了5月份預見期為6個月和9月份預見期為3個月時，篩選出的對天一水庫入庫徑流影響最為顯著的海溫偶極，考慮文章篇幅問題，對其他月份各預見期對應的海溫偶極篩選結果不一一展示。由圖可知，5 月份6 個月預見期對應的正偶極子主要分布在太平洋和印度洋低緯度地區(qū)，負偶極子分布在東經(jīng)40°至65°之間（印度洋）和東經(jīng)180°（東太平洋）附近。9月份3 個月預見期對應的偶極在太平洋海域分布差別較大，沒有明顯的區(qū)域聚集規(guī)律。對于5 月份6 個月預見期而言，挑選出的偶極與Nino 34、Nino 4 海域有所重疊；對于9 月份3 個月預見期而言，北太平洋白令海峽附近海溫異?，F(xiàn)象重復出現(xiàn)，兩種情況下對應的最顯著的偶極均主要分布在印度洋海域內(nèi)。由于海洋變化、大氣環(huán)流及其交互作用十分復雜，不同海域的溫度場之間的關聯(lián)性和空間關系難以歸納出一般規(guī)律，我們較難對特定研究區(qū)域解釋選取的正負偶極呈現(xiàn)的空間規(guī)律與關系的物理成因。本研究的重點不在此，但相關解釋在未來的研究中應完善。

圖2 對預見期為6個月時的5月份和預見期為3個月時的9月份篩選出的偶極位置分布圖Fig.2 Selected SST dipoles for May at the 6-month lead time and for September at the 3-month lead time

考慮到天氣系統(tǒng)的混沌性和隨機性，引入不同來源的因子可能會產(chǎn)生噪聲干擾從而影響預報精度，本研究向多元線性回歸模型中逐步引入前期徑流、大氣環(huán)流和海溫偶極因子，得到多種因子組合方案，如圖3所示。圖中以預報精度為評價指標，預報精度越高，顏色越藍，反之則越紅。由于僅考慮前期徑流進行預報時只有個別月份能達到率定期修正復相關系數(shù)大于0.65 的要求，故在圖中不予展示。圖3（a）和（b）分別為預報因子只考慮大氣環(huán)流和同時考慮徑流和大氣環(huán)流因子的情況。圖3（c）到（f）為預報因子含有海溫偶極的情況。由圖可知，預報因子只考慮大氣環(huán)流或考慮徑流和大氣環(huán)流因子時，在1 月份至4 月份的12 個預見期內(nèi)仍可以達到較高的預報精度。以預報因子考慮徑流和大氣環(huán)流因子為例，1 月至4 月份12 個預見期的平均值分別為71.7、75.3%、65.8%和70.2%；而在5 月份、9月份和11 月的預報精度降至47.4%、47.9%和44.0%，特別是9月份和11月份，有多個預見期的預報精度低于40%。當預報因子考慮海溫偶極因子、徑流和海溫偶極因子、大氣和海溫偶極因子以及徑流、大氣和海溫偶極因子時［對應圖3（c）到（f）］，9月份12 個月預見期的平均精度提升至了54.3%、55.0%、54.9%和55.6%，11 月份12 個月預見期的平均精度提升至了53.3%、53.3%、49.7%和49.7%。

圖3 多元線性回歸模型在預見期1-12個月的預報精度Fig.3 Comparison of prediction accuracy at lead times of 1 to 12 months among different predictor combinations

圖4展示了6種預報因子組合在驗證期12個月預見期的平均徑流預報精度。預見期1 到3 個月時，預報效果最佳的因子組合為考慮徑流和大氣環(huán)流或考慮徑流和海溫偶極，3 個預見期2 種組合的預報精度分別為66.7%、66.8%，63.4%、63.9%，62.5%、62.3%。由此可見，盡管前期徑流因子不能單獨用作預報，但它的加入能提高大氣環(huán)流因子和海溫偶極因子單獨使用時的預報精度，主要體現(xiàn)在預見期較短時。而含有海溫偶極因子的預報因子組合，在預見期較長時預報效果較優(yōu)。僅含有海溫偶極因子或海溫偶極與前期徑流因子的組合，在預見期4 到11個月預報精度最優(yōu)。當同時考慮徑流、大氣環(huán)流和海溫偶極因子時，在各預見期的預報效果均不是最佳，可能是由于過擬合所導致的。

圖4 各種預報因子組合在驗證期的12個月預見期的平均預報精度Fig.4 Prediction accuracy of all predictor combinations at 1- to 12-month lead times in the validation period.

圖5（a）和圖5（b）分別展示了不同預報因子組合方式（與圖3順序一致）在天一水庫汛期（5月至10月）和非汛期（11月至次年4 月）的徑流預報效果。如圖所示，在汛期，海溫偶極對徑流預測的指示性比大氣環(huán)流因子強。過往研究表明中國夏季降水主要受夏季風異常影響，該夏季風來自印度洋，強度受印度洋海溫異常影響［27-30］，可認為中國大部分流域汛期徑流的變化與印度洋海溫變化息息相關。圖2 表明影響5 月和9 月的徑流變化的最顯著的偶極均主要分布在印度洋海域內(nèi)，與相關文獻一致。對于非汛期而言，在預見期4個月以內(nèi)時，大氣環(huán)流因子對徑流預測的指示性遠強于海溫偶極因子，而當預見期超過4個月后，海溫的指示作用超過大氣環(huán)流。其主要原因是海洋記憶的時間尺度較長，海洋的溫度變化遠慢于大氣，故海溫異常對徑流的作用的滯后時間要長于大氣環(huán)流因子。從圖中還可知，當預見期較短時（汛期預見期2 個月，非汛期預見期4 個月），引入流域的前期徑流作為預報因子之一可以有效地提升預報效果，而當預見期較長時，引入前期徑流對徑流預測的指示作用消失。主要原因是通常當滯時超過3 個月時，月徑流的偏自相關性就不再顯著，可認為2 到4 個月前的徑流對當前月份徑流的影響很小，引入其作為預報因子無法提升預報效果。

圖5 各種預報因子組合在汛期和非汛期的預報精度Fig.5 Prediction accuracy of all predictor combinations at 1 to 12 month lead times in the flood season and non-flood season respectively

圖6 展示了由2012 年12 月發(fā)起的2013 年全年的作業(yè)預報和由2017 年12 月發(fā)起的2018 年全年的作業(yè)預報徑流過程線圖，其中，黑色實線代表實測徑流過程線。由圖可知，各種因子組合預報徑流量區(qū)別最大的月份為7 月至11 月。在預報2013年和2018年徑流時，是否考慮徑流作為預報因子對預報結果的影響較小?？傮w而言，預報因子考慮徑流和海溫偶極與僅考慮海溫偶極的預報效果最好，模擬出了2018年汛期徑流的雙峰型特征。

圖6 由2012年12月和2017年12月發(fā)起的12個月預見期的徑流預報過程線圖Fig.6 Streamflow predictions for 2013 and 2018 issued in December 2012 and December 2017， respectively

4 結 論

研究以天一水庫為研究對象，基于多元線性回歸方法在年內(nèi)各月分別構建了預見期為1～12 個月的徑流預報模型。在關鍵因子挑選部分以海溫偶極取代了傳統(tǒng)的固定海域的海洋表面溫度場，并綜合考慮了水庫前期徑流和大氣環(huán)流對月徑流預報的影響，主要結論如下。

（1）僅考慮前期徑流作為預報因子時，月徑流預報精度較差，但它的加入能提高大氣環(huán)流因子和海溫偶極因子單獨使用時的預報精度，主要體現(xiàn)在預見期為1～3個月時。

（2）含有海溫偶極因子的預報因子組合，在預見期較長時預報效果優(yōu)于預報因子只考慮大氣環(huán)流因子和考慮徑流和大氣環(huán)流因子的模型，其中以徑流和海溫偶極為預報因子的模型表現(xiàn)效果最佳，其在預見期為4～9 個月的平均精度較前述模型分別提升了2.4%和2.1%。此外，預報精度提升效果最為顯著的月份為9月和11月，該模型對這兩個月份在預見期1～12個月的平均精度較前述模型分別提升了8.2%、7.1%和12.6%、9.3%。