潘家口水庫上游面雨量氣候特征及入庫水量預(yù)測模型探究

王朋朋，薛思嘉，易永力，瞿衛(wèi)東，張海英

（1.承德市氣象局，河北承德067000；2.引灤工程管理局水文水質(zhì)監(jiān)測中心，河北遷西063009）

潘家口水庫位于河北省寬城滿族自治縣、興隆 縣與遷西縣交界處，于1975年動工，至1979年竣工。潘家口水庫為引灤入津的主體工程，以向天津供水為主，兼顧發(fā)電與防洪功能。水庫流域面積為44 600 km2，其中壩址以上流域面積33 700 km2，占灤河流域面積的75%。灤河流域水系發(fā)達(dá)，比較大的支流有閃電河、螞蟻吐河、伊遜河、柳河、瀑河等。流域地形差異比較大，總體變化趨勢是西北高東南低，海拔高度為1 300～1 800 m。上游區(qū)域多年年平均降水量572.2 mm，平均氣溫為8.4℃，降水量與氣溫分布的特征是隨緯度的增高而遞減。

流域面雨量是指一定區(qū)域內(nèi)的實際平均降水量，能客觀地描述該區(qū)域?qū)嶋H降水資源狀況，它作為洪水與水庫調(diào)度中非常重要的參數(shù)，是各級政府組織防汛抗洪和水庫調(diào)度等決策的重要依據(jù)，對流域面雨量進(jìn)行氣候特征分析及實際業(yè)務(wù)預(yù)測也是氣象部門拓展服務(wù)領(lǐng)域的新舉措。

面雨量是傳統(tǒng)洪水預(yù)報非常重要的輸入強迫因子之一。在洪水預(yù)報實時作業(yè)中，面雨量估測與預(yù)報的精確度在較大程度上決定了洪水預(yù)報的精度。因此，為了能更好地為各級政府等部門提供防汛抗洪、水庫調(diào)度等信息，面雨量的研究工作就顯得尤為重要，同時這也是氣象部門和水利部門加強合作的重要途徑。

本研究利用潘家口水庫提供的長時間序列的入庫水量數(shù)據(jù)，結(jié)合水庫上游流域6個國家氣象站的氣象數(shù)據(jù)，通過分析潘家口水庫上游流域氣候特征、潘家口水庫多年面雨量時空分布特征，為指導(dǎo)業(yè)務(wù)人員更好地開展水庫氣象服務(wù)工作提供科學(xué)依據(jù)，為更有效地提高河北省防汛抗旱和水資源調(diào)度氣象服務(wù)能力提供決策參考。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)來源

本文常規(guī)氣象數(shù)據(jù)來源于水庫上游區(qū)域隆化縣、灤平縣、承德市、承德縣、興隆縣和寬城縣6個國家氣象站數(shù)據(jù)，主要包括逐日氣溫、逐日降水量（20—20時），資料長度為1965—2019年。潘家口入庫水量數(shù)據(jù)來源于海委引灤工程管理局潘家口，其中1980—2017年為年入庫數(shù)據(jù)、2017—2019年部分?jǐn)?shù)據(jù)為逐日入庫水量。

1.2 數(shù)據(jù)處理方法

1.2.1 水庫上游流域氣候特征變化分析方法

氣象數(shù)據(jù)的預(yù)處理在Excel中完成。采用Mann-Kendall（以下簡稱M-K）檢驗法檢測各氣象要素序列的變化趨勢[1]，通過檢驗可得統(tǒng)計量Z值，Z為正值表示增加趨勢、負(fù)值表示減少趨勢。Z的絕對值在大于等于1.28、1.64、2.32時表示分別通過了信度90%、95%、99%顯著性檢驗。水庫上游極端降水概率的檢驗，采用皮爾遜-Ⅲ型分布模型[2]。

除了使用幾種常用的面雨量計算方法以外，近些年，計算面雨量新的方法不斷被提出，比如網(wǎng)格法、基于數(shù)值預(yù)報產(chǎn)品法、基于函數(shù)插值的最優(yōu)插值和薄板平滑插值法等[3-6]。

本文采用基于Gis軟件的泰森多邊形法（Thiessen Polygons）來計算潘家口水庫流域的面雨量[7]，計算公式如下：

式中：P為流域平均面雨量（mm）；fi為各雨量站用多邊形面積計算的權(quán)重數(shù)；Pi為各測站同時期降雨量（mm）。水庫入庫水量的回歸方程采用逐步回歸法，利用R2.3中MASS包中的“stepAIC（）”函數(shù)實現(xiàn)逐步回歸模型。

2 結(jié)果與分析

2.1 水庫上游流域氣候特征變化

水庫上游流域近55 a年均降雨量，如圖1（a）所示。水庫上游年均降水量為571.7 mm，年均最小降水量為381.2 mm，出現(xiàn)在2002年；年均最大降水量為826.0 mm，出現(xiàn)在1978年。其中，夏季（6—8月）降水量占年總降水量的68.7%，主汛期（7—8月）降水量占年總降水量的52.9%，降水集中。從圖1（a）可知，水庫上游流域年均降水量呈現(xiàn)出逐年弱減少的趨勢。以M-K檢驗法檢測水庫上游降水統(tǒng)計量Z=-0.75，表明水庫上游近55 a降雨量下降趨勢不明顯。

水庫上游區(qū)域雨日強度計算方法為上游區(qū)域年降水總量與降水日數(shù)之比，單位為mm/d。水庫上游雨日強度在5.3～10.0 mm/d，近55 a平均雨日強度為7.3 mm/d，最大、最小雨日強度分別出現(xiàn)在1994、1980年。由圖1（b）可得，近55 a水庫上游雨日強度呈現(xiàn)上升趨勢，經(jīng)M-K趨勢檢驗得Z=1.99＞1.64，即通過了95%顯著性檢驗，說明水庫上游近55 a雨日強度上升趨勢明顯。降水強度在增加，也預(yù)示著極端降水事件發(fā)生的概率在增加。

當(dāng)日20時至次日20時降水量大于25、50 mm時，分別記為1 d大雨日、暴雨日。水庫上游年平均大雨日數(shù)為4.2 d，年平均暴雨日數(shù)為1.9 d。其中，平均大雨日數(shù)最大、最小值分別出現(xiàn)在1995、2002年，分別為7、1.6 d；平均暴雨日數(shù)最大、最小值分別出現(xiàn)在2018、2002年，分別為3.8、1 d。年平均大雨、暴雨日數(shù)都有趨于減少的趨勢，經(jīng)M-K檢驗可得大雨日數(shù)、暴雨日數(shù)減少趨勢不明顯，如圖1（c）所示。

③明挖段隧道為出地面結(jié)構(gòu)，場地所在橋下現(xiàn)狀地面為一南北走向斜坡，斜坡高差達(dá)20m，樁基托換在地面上實施，施工作業(yè)場地困難。

水庫上游年平均氣溫為8.47℃，最大、最小值分別出現(xiàn)在1998、1969年，分別為9.56℃、6.59℃，其余年份在平均值上下波動。從圖1（d）可以看出，水庫上游年均氣溫呈現(xiàn)出上升趨勢，經(jīng)M-K趨勢檢驗可得Z=3.83＞2.32，表明通過了信度99%顯著性檢驗，因此可以說水庫上游平均氣溫上升趨勢明顯。水庫區(qū)域氣溫變化與氣候變暖的背景相一致，氣溫變暖導(dǎo)致蒸發(fā)量的增加，會進(jìn)一步影響入庫水量。

圖1 1965—2019年水庫上游流域氣候特征變化

2.2 水庫上游流域極端降水概率分布

通過皮爾遜-Ⅲ分布函數(shù)計算，得出水庫上游年降水量出現(xiàn)各種概率的重現(xiàn)水平及不同概率級別出現(xiàn)的年份，其中10、20、50、100 a一遇分別指的是強降水概率為10%、5%、2%、1%的可能最大降水量。

水庫上游年降水量極值重現(xiàn)水平及出現(xiàn)時間，詳見表1。由表1可知，近55 a水庫上游年降水量未出現(xiàn)100 a一遇的降水；50 a一遇的降水出現(xiàn)在1978年，降水量為826.0 mm；20 a一遇的降水出現(xiàn)在1979年，降水量為748.6 mm；10 a一遇的降水出現(xiàn)4次，分別在1969、1973、1994、1998年。極端降水事件多出現(xiàn)在20世紀(jì)70—90年代，近10 a未出現(xiàn)極端降水事件。

表1 潘家口水庫上游流域年降水量極值重現(xiàn)水平及出現(xiàn)時間

2.3 水庫上游流域面雨量年際變化特征

水庫上游流域近55 a面雨量及其距平歷年變化，如圖2所示。通過計算知道，近55 a水庫上游流域年均面雨量為585.9 mm，面雨量最大值出現(xiàn)在1978年，為842.3 mm；最小值出現(xiàn)在2002年，為387.7 mm。近55 a面雨量呈現(xiàn)波動下降趨勢，下降速率為8.7 mm/10 a。用降水距平百分率來衡量旱澇的發(fā)生發(fā)展情況，定義降水量距平百分率≤－25%為旱、≥25%為澇。由圖2可以看出，近55 a水庫上游流域發(fā)生旱年份略低于澇年份，達(dá)到干旱指標(biāo)的有3 a，多集中在20世紀(jì)70年代和21世紀(jì)初，其中2002年干旱程度最明顯，降水距平百分率達(dá)-38%；達(dá)到水澇指標(biāo)的有5 a，多集中在20世紀(jì)70年代末和90年代，其中1978年的降水距平百分率高達(dá)44%，這是流域自有氣象記錄以來降水量最多的一年。

圖2 潘家口水庫上游流域面雨量年變化

2.4 水庫上游流域面雨量月際變化特征

水庫上游流域近55 a面雨量月變化特征分布，如圖3所示。水庫上游流域夏、秋季6—9月各月面雨量都在50 mm以上，遠(yuǎn)高于其他月份，為流域的多雨季。多雨季占常年平均雨量的78%。從各月降水變化趨勢來看，除了1—2、7—8月有逐年減少的趨勢外，其他月份均有逐年增加的趨勢。經(jīng)過M-K趨勢檢驗分析可知，8月Z值為2.34＞2.32，即8月面雨量減少趨勢顯著，通過了信度為99%的檢驗；10月面雨量降水呈顯著增加趨勢，通過了信度95%的檢驗。7—8月雖然處于水庫主降水期，但面雨量降水有逐年減少的趨勢，主汛期前后的時間段（4—6、9—11月）面雨量降水有逐年增加的趨勢，仍要注意做好這個時間段的水庫入庫水量預(yù)測工作，為水庫安全運行提供氣象服務(wù)保障。

圖3 潘家口水庫上游流域面雨量月變化

2.5 水庫上游流域不同量級面雨量分布特征

在以上研究的基礎(chǔ)上，本文利用近55 a水庫上游流域逐日面雨量資料，并參考我國江河面雨量等級劃分標(biāo)準(zhǔn)[8]（小雨0.5～5.9 mm、中雨6.0～14.9 mm、大雨15.0～29.9 mm、暴雨30～59.9 mm、大暴雨60～149.9 mm、特大暴雨≥150 mm）統(tǒng)計分析1965—2019年水庫上游流域不同量級降水日數(shù)及降水量的分布情況，發(fā)現(xiàn)近55 a流域年均雨日為76 d，其中小雨日49 d，占總雨日的64%；中雨日16.5 d，占總雨日的22%；大雨日7 d，占總雨日的9%；暴雨日3 d；大暴雨日0.5 d；無特大暴雨日。

1965—2019 年水庫上游流域不同量級降水日數(shù)在各個月份的分布情況，詳見表2。從各月雨日來看，5—9月雨日居多，最多的7月平均雨日達(dá)15.8 d，6月為13.3 d，最少的1月不足1 d。從不同量級降水在各月的分布情況看，小雨和中雨量級降水多出現(xiàn)在6、7、8月，大雨及以上量級降水集中在7、8月。各降水量級最多的天數(shù)均出現(xiàn)在7月，因此7月前后兩個月即6—9月是水庫入庫水量預(yù)測的重點時間段。

2.6 水庫上游流域3 d以上集中強降水特征

根據(jù)盧珊[9]等研究成果，定義水庫流域連續(xù)3 d以上日雨量大于等于20 mm的連續(xù)降水過程為集中強降水過程。表3統(tǒng)計了1965—2019年水庫上游流域各月集中強降水過程的各參數(shù)，近55 a水庫上游流域共出現(xiàn)8次強降水過程。從各月份分布來看，除7—9月外其他月份均未出現(xiàn)集中強降水過程，7月為強降水發(fā)生最多的月份（4次），平均過程雨量為95.2 mm；8月發(fā)生集中強降水過程的次數(shù)僅次于7月，但強度較7月更大，最大過程雨量達(dá)到151.3 mm。9月也有集中強降水過程的發(fā)生，但強度較夏季有明顯的減弱，進(jìn)入9月一般來說降雨量會減少，但仍要防范強降水對水庫蓄水的影響。

表3 1965—2019年潘家口水庫上游流域各月集中強降水過程參數(shù)統(tǒng)計

2.7 水庫上游流域面雨量與入庫水量的關(guān)系

根據(jù)上述可知，水庫上游流域夏季降水量占年總降水量的68.7%，降水主要集中在此期間，冬季降水量占比最少，只占到1.5%。因此，為研究水庫入庫水量與上游流域降水量的關(guān)系，篩選出2017—2020年夏季上游區(qū)域6個國家站的日降水量（20—20時），然后篩選出中雨以上的降水?dāng)?shù)據(jù)，即6個國家站中有一個站日降水量大于10 mm就記為一列有效數(shù)據(jù)，篩選出同一天的入庫水量匯總成數(shù)據(jù)集。面雨量與入庫水量的關(guān)系如圖4所示，兩者有很強的正相關(guān)關(guān)系，相關(guān)系數(shù)為0.67，即隨著面雨量的增大，入庫水量隨之增大。

圖4 潘家口水庫上游流域面雨量與入庫水量的關(guān)系

2.8 水庫上游流域入庫水量研究

針對潘家口水庫上游區(qū)域，從影響入庫徑流的氣象因子入手，通過深入分析入庫水量與各因子的響應(yīng)關(guān)系，以期建立基于氣象因子的入庫水量預(yù)測模型，從而與中長期天氣預(yù)報和短期氣候預(yù)測結(jié)合，為水庫提供專業(yè)氣象服務(wù)。

由圖5可以看出，1980—2017年水庫入庫水量是以2.5萬m3/10 a的速率在減少，這與水庫上游流域年降水量變化趨勢相一致。從表4可以看出，入庫水量與上游年均降水量、大雨日數(shù)、暴雨日數(shù)、雨日強度呈正相關(guān)性，且均通過了顯著性檢驗（P＜0.01），入庫水量與上游年均氣溫呈負(fù)相關(guān)關(guān)系但未通過顯著性檢驗。

圖5 潘家口水庫近55 a入庫水量變化

表4 入庫水量與相關(guān)因子間相關(guān)系數(shù)及顯著性

篩選出2017—2020年（7—9月）上游流域6個國家站的日降水量（20—20時），然后篩選出中雨以上的降水?dāng)?shù)據(jù)，即6個國家站中有一個站日降水量大于10 mm就記為一列有效數(shù)據(jù)，篩選出同一天的入庫水量匯總成數(shù)據(jù)集。選擇2017—2019年的數(shù)據(jù)進(jìn)行建模，以2020年的數(shù)據(jù)來驗證模型的準(zhǔn)確性。

因為距離的問題，考慮到并不是每個國家站當(dāng)日的降水量都會對當(dāng)日的入庫水量有所影響，例如距離水庫較遠(yuǎn)的隆化縣，如果當(dāng)日有降水，可能會在2 d或者3 d之后對入庫水量有影響，即入庫水量和降水量存在滯后性，將水庫上游6個站的降水量來構(gòu)建入庫水量預(yù)測模型，計算滯后1 d時用當(dāng)日的降水量作為自變量因子、后一天的入庫水量作為因變量，滯后2、3 d計算方法同滯后1 d。采用逐步回歸的方法逐步剔除不顯著的因子，從而得到入庫水量的最優(yōu)預(yù)報回歸模型。

利用AIC（赤池信息準(zhǔn)則）來比較預(yù)測模型，它考慮到了模型的統(tǒng)計擬合度以及用來擬合的參數(shù)數(shù)目。AIC越小的模型越優(yōu)先選擇。根據(jù)表5可得，選擇當(dāng)日降水量與當(dāng)日入庫水量時AIC值最小，因此以此建模較好。

表5 入庫水量與降水量滯后性模型AIC值

經(jīng)過逐步回歸分析，程序自動剔除了承德市、灤平、隆化、承德縣4個站，得到的水庫入庫水量與上游流域降水量的模擬方程為：

式中：Q入庫量為入庫水量；xl，kc分別代表興隆，寬城國家站的降水量數(shù)據(jù)。為了驗證方程的模擬效果，利用2020年汛期降水日的降水?dāng)?shù)據(jù)，結(jié)合同時期的入庫水量數(shù)據(jù)，擬合出了2020年7—9月潘家口水庫入庫水量模擬值。實測值和模擬值的變化曲線如圖6所示，由此可以看出實測值和模擬值曲線的變化趨勢基本相同，實測值和模擬值大部分能重合。這說明利用興隆和寬城的降水?dāng)?shù)據(jù)能較好地模擬出水庫入庫水量的變化趨勢，擬合效果較好。

圖6 2020年7—9月入庫水量實測值和模擬值變化曲線

為了進(jìn)一步驗證上述回歸方程模型的實用性、檢驗2020年7—9月入庫水量實測值和模擬值之間的誤差，采用均方根誤差（RMSE）、平均絕對誤差（MAE）和Theil系數(shù)等來檢驗?zāi)Ｐ偷男Ч蛯崪y與模擬之間的誤差。其計算公式如下：

式中：ai表示實測值；mi表示模擬值。RMSE反映了模擬值偏離實測值的程度，RMSE越小表示模擬精度越高，因此可用RMSE作為評定預(yù)測模型精度的標(biāo)準(zhǔn)。Theil不等系數(shù)能描述兩組數(shù)據(jù)的偏差，不等系數(shù)越小表明實測值和模擬值越接近。

由表6可知，根據(jù)入庫水量的預(yù)測模型得出的均方根誤差（RMSE）為206.13、不等系數(shù)為0.23。模擬值與實測值偏差不大，方程的模擬效果較好。

表6 2020年7—9月入庫水量模擬值與實測值偏差

3 結(jié)論

（1）潘家口水庫上游流域近55 a年均降水量、大雨、暴雨日數(shù)呈現(xiàn)逐年略減少的趨勢；雨日強度呈現(xiàn)逐年上升的趨勢顯著，降水強度的增加也預(yù)示著極端降水事件發(fā)生的概率在增加；水庫上游流域氣溫升高趨勢顯著，可能導(dǎo)致蒸發(fā)、蒸散的增強，進(jìn)而導(dǎo)致水庫蓄水量的減少。7—8月雖然處于水庫主汛期，但面雨量有逐年減少的趨勢，主汛期前后的時間段（4—6、9—11月）面雨量降水有逐年增加的趨勢，要做好這個時間段的水庫入庫水量預(yù)測工作，為水庫安全運行提供氣象服務(wù)保障。7月為強降水過程發(fā)生次數(shù)最多的月份，8月集中強降水過程的次數(shù)僅次于7月，但強度較7月更大，9月也有集中強降水過程的發(fā)生，但次數(shù)和強度較7—8月有明顯的減弱，進(jìn)入9月一般來說降雨量會減少，但仍要防范強降水過程對水庫蓄水的影響。

（2）本文采用6個國家站的降水量來計算流域面雨量，有研究表明插值雨量站點越多降雨插值結(jié)果的精度越高，因此在以后的研究中應(yīng)增加觀測站的數(shù)量。本文僅是基于數(shù)學(xué)統(tǒng)計方法對水庫入庫水量和上游降水量進(jìn)行簡單的數(shù)學(xué)建模，但影響入庫水量的機理性研究不深入，例如蒸發(fā)、蒸散、地表徑流等都可能對入庫水量產(chǎn)生影響。

（3）潘家口水庫上游流域有26個水文站，觀測時間序列較長，積累了長序列觀測資料，可以研究其氣候變化特征。要加快氣象、水文部門信息資料共享進(jìn)程，建立流域信息匯集和共享的機制和渠道，可以通過共同建設(shè)實時數(shù)據(jù)庫、資料存儲和發(fā)布系統(tǒng)，實現(xiàn)流域氣象、水文信息共享。