潘家口水库上游面雨量气候特征及入库水量预测模型探究

王朋朋，薛思嘉，易永力，瞿卫东，张海英

（1.承德市气象局，河北承德067000；2.引滦工程管理局水文水质监测中心，河北迁西063009）

潘家口水库位于河北省宽城满族自治县、兴隆 县与迁西县交界处，于1975年动工，至1979年竣工。潘家口水库为引滦入津的主体工程，以向天津供水为主，兼顾发电与防洪功能。水库流域面积为44 600 km2，其中坝址以上流域面积33 700 km2，占滦河流域面积的75%。滦河流域水系发达，比较大的支流有闪电河、蚂蚁吐河、伊逊河、柳河、瀑河等。流域地形差异比较大，总体变化趋势是西北高东南低，海拔高度为1 300～1 800 m。上游区域多年年平均降水量572.2 mm，平均气温为8.4℃，降水量与气温分布的特征是随纬度的增高而递减。

流域面雨量是指一定区域内的实际平均降水量，能客观地描述该区域实际降水资源状况，它作为洪水与水库调度中非常重要的参数，是各级政府组织防汛抗洪和水库调度等决策的重要依据，对流域面雨量进行气候特征分析及实际业务预测也是气象部门拓展服务领域的新举措。

面雨量是传统洪水预报非常重要的输入强迫因子之一。在洪水预报实时作业中，面雨量估测与预报的精确度在较大程度上决定了洪水预报的精度。因此，为了能更好地为各级政府等部门提供防汛抗洪、水库调度等信息，面雨量的研究工作就显得尤为重要，同时这也是气象部门和水利部门加强合作的重要途径。

本研究利用潘家口水库提供的长时间序列的入库水量数据，结合水库上游流域6个国家气象站的气象数据，通过分析潘家口水库上游流域气候特征、潘家口水库多年面雨量时空分布特征，为指导业务人员更好地开展水库气象服务工作提供科学依据，为更有效地提高河北省防汛抗旱和水资源调度气象服务能力提供决策参考。

1 数据与方法

1.1 数据来源

本文常规气象数据来源于水库上游区域隆化县、滦平县、承德市、承德县、兴隆县和宽城县6个国家气象站数据，主要包括逐日气温、逐日降水量（20—20时），资料长度为1965—2019年。潘家口入库水量数据来源于海委引滦工程管理局潘家口，其中1980—2017年为年入库数据、2017—2019年部分数据为逐日入库水量。

1.2 数据处理方法

1.2.1 水库上游流域气候特征变化分析方法

气象数据的预处理在Excel中完成。采用Mann-Kendall（以下简称M-K）检验法检测各气象要素序列的变化趋势[1]，通过检验可得统计量Z值，Z为正值表示增加趋势、负值表示减少趋势。Z的绝对值在大于等于1.28、1.64、2.32时表示分别通过了信度90%、95%、99%显著性检验。水库上游极端降水概率的检验，采用皮尔逊-Ⅲ型分布模型[2]。

除了使用几种常用的面雨量计算方法以外，近些年，计算面雨量新的方法不断被提出，比如网格法、基于数值预报产品法、基于函数插值的最优插值和薄板平滑插值法等[3-6]。

本文采用基于Gis软件的泰森多边形法（Thiessen Polygons）来计算潘家口水库流域的面雨量[7]，计算公式如下：

式中：P为流域平均面雨量（mm）；fi为各雨量站用多边形面积计算的权重数；Pi为各测站同时期降雨量（mm）。水库入库水量的回归方程采用逐步回归法，利用R2.3中MASS包中的“stepAIC（）”函数实现逐步回归模型。

2 结果与分析

2.1 水库上游流域气候特征变化

水库上游流域近55 a年均降雨量，如图1（a）所示。水库上游年均降水量为571.7 mm，年均最小降水量为381.2 mm，出现在2002年；年均最大降水量为826.0 mm，出现在1978年。其中，夏季（6—8月）降水量占年总降水量的68.7%，主汛期（7—8月）降水量占年总降水量的52.9%，降水集中。从图1（a）可知，水库上游流域年均降水量呈现出逐年弱减少的趋势。以M-K检验法检测水库上游降水统计量Z=-0.75，表明水库上游近55 a降雨量下降趋势不明显。

水库上游区域雨日强度计算方法为上游区域年降水总量与降水日数之比，单位为mm/d。水库上游雨日强度在5.3～10.0 mm/d，近55 a平均雨日强度为7.3 mm/d，最大、最小雨日强度分别出现在1994、1980年。由图1（b）可得，近55 a水库上游雨日强度呈现上升趋势，经M-K趋势检验得Z=1.99＞1.64，即通过了95%显著性检验，说明水库上游近55 a雨日强度上升趋势明显。降水强度在增加，也预示着极端降水事件发生的概率在增加。

当日20时至次日20时降水量大于25、50 mm时，分别记为1 d大雨日、暴雨日。水库上游年平均大雨日数为4.2 d，年平均暴雨日数为1.9 d。其中，平均大雨日数最大、最小值分别出现在1995、2002年，分别为7、1.6 d；平均暴雨日数最大、最小值分别出现在2018、2002年，分别为3.8、1 d。年平均大雨、暴雨日数都有趋于减少的趋势，经M-K检验可得大雨日数、暴雨日数减少趋势不明显，如图1（c）所示。

③明挖段隧道为出地面结构，场地所在桥下现状地面为一南北走向斜坡，斜坡高差达20m，桩基托换在地面上实施，施工作业场地困难。

水库上游年平均气温为8.47℃，最大、最小值分别出现在1998、1969年，分别为9.56℃、6.59℃，其余年份在平均值上下波动。从图1（d）可以看出，水库上游年均气温呈现出上升趋势，经M-K趋势检验可得Z=3.83＞2.32，表明通过了信度99%显著性检验，因此可以说水库上游平均气温上升趋势明显。水库区域气温变化与气候变暖的背景相一致，气温变暖导致蒸发量的增加，会进一步影响入库水量。

图1 1965—2019年水库上游流域气候特征变化

2.2 水库上游流域极端降水概率分布

通过皮尔逊-Ⅲ分布函数计算，得出水库上游年降水量出现各种概率的重现水平及不同概率级别出现的年份，其中10、20、50、100 a一遇分别指的是强降水概率为10%、5%、2%、1%的可能最大降水量。

水库上游年降水量极值重现水平及出现时间，详见表1。由表1可知，近55 a水库上游年降水量未出现100 a一遇的降水；50 a一遇的降水出现在1978年，降水量为826.0 mm；20 a一遇的降水出现在1979年，降水量为748.6 mm；10 a一遇的降水出现4次，分别在1969、1973、1994、1998年。极端降水事件多出现在20世纪70—90年代，近10 a未出现极端降水事件。

表1 潘家口水库上游流域年降水量极值重现水平及出现时间

2.3 水库上游流域面雨量年际变化特征

水库上游流域近55 a面雨量及其距平历年变化，如图2所示。通过计算知道，近55 a水库上游流域年均面雨量为585.9 mm，面雨量最大值出现在1978年，为842.3 mm；最小值出现在2002年，为387.7 mm。近55 a面雨量呈现波动下降趋势，下降速率为8.7 mm/10 a。用降水距平百分率来衡量旱涝的发生发展情况，定义降水量距平百分率≤－25%为旱、≥25%为涝。由图2可以看出，近55 a水库上游流域发生旱年份略低于涝年份，达到干旱指标的有3 a，多集中在20世纪70年代和21世纪初，其中2002年干旱程度最明显，降水距平百分率达-38%；达到水涝指标的有5 a，多集中在20世纪70年代末和90年代，其中1978年的降水距平百分率高达44%，这是流域自有气象记录以来降水量最多的一年。

图2 潘家口水库上游流域面雨量年变化

2.4 水库上游流域面雨量月际变化特征

水库上游流域近55 a面雨量月变化特征分布，如图3所示。水库上游流域夏、秋季6—9月各月面雨量都在50 mm以上，远高于其他月份，为流域的多雨季。多雨季占常年平均雨量的78%。从各月降水变化趋势来看，除了1—2、7—8月有逐年减少的趋势外，其他月份均有逐年增加的趋势。经过M-K趋势检验分析可知，8月Z值为2.34＞2.32，即8月面雨量减少趋势显著，通过了信度为99%的检验；10月面雨量降水呈显著增加趋势，通过了信度95%的检验。7—8月虽然处于水库主降水期，但面雨量降水有逐年减少的趋势，主汛期前后的时间段（4—6、9—11月）面雨量降水有逐年增加的趋势，仍要注意做好这个时间段的水库入库水量预测工作，为水库安全运行提供气象服务保障。

图3 潘家口水库上游流域面雨量月变化

2.5 水库上游流域不同量级面雨量分布特征

在以上研究的基础上，本文利用近55 a水库上游流域逐日面雨量资料，并参考我国江河面雨量等级划分标准[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）统计分析1965—2019年水库上游流域不同量级降水日数及降水量的分布情况，发现近55 a流域年均雨日为76 d，其中小雨日49 d，占总雨日的64%；中雨日16.5 d，占总雨日的22%；大雨日7 d，占总雨日的9%；暴雨日3 d；大暴雨日0.5 d；无特大暴雨日。

1965—2019 年水库上游流域不同量级降水日数在各个月份的分布情况，详见表2。从各月雨日来看，5—9月雨日居多，最多的7月平均雨日达15.8 d，6月为13.3 d，最少的1月不足1 d。从不同量级降水在各月的分布情况看，小雨和中雨量级降水多出现在6、7、8月，大雨及以上量级降水集中在7、8月。各降水量级最多的天数均出现在7月，因此7月前后两个月即6—9月是水库入库水量预测的重点时间段。

2.6 水库上游流域3 d以上集中强降水特征

根据卢珊[9]等研究成果，定义水库流域连续3 d以上日雨量大于等于20 mm的连续降水过程为集中强降水过程。表3统计了1965—2019年水库上游流域各月集中强降水过程的各参数，近55 a水库上游流域共出现8次强降水过程。从各月份分布来看，除7—9月外其他月份均未出现集中强降水过程，7月为强降水发生最多的月份（4次），平均过程雨量为95.2 mm；8月发生集中强降水过程的次数仅次于7月，但强度较7月更大，最大过程雨量达到151.3 mm。9月也有集中强降水过程的发生，但强度较夏季有明显的减弱，进入9月一般来说降雨量会减少，但仍要防范强降水对水库蓄水的影响。

表3 1965—2019年潘家口水库上游流域各月集中强降水过程参数统计

2.7 水库上游流域面雨量与入库水量的关系

根据上述可知，水库上游流域夏季降水量占年总降水量的68.7%，降水主要集中在此期间，冬季降水量占比最少，只占到1.5%。因此，为研究水库入库水量与上游流域降水量的关系，筛选出2017—2020年夏季上游区域6个国家站的日降水量（20—20时），然后筛选出中雨以上的降水数据，即6个国家站中有一个站日降水量大于10 mm就记为一列有效数据，筛选出同一天的入库水量汇总成数据集。面雨量与入库水量的关系如图4所示，两者有很强的正相关关系，相关系数为0.67，即随着面雨量的增大，入库水量随之增大。

图4 潘家口水库上游流域面雨量与入库水量的关系

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

针对潘家口水库上游区域，从影响入库径流的气象因子入手，通过深入分析入库水量与各因子的响应关系，以期建立基于气象因子的入库水量预测模型，从而与中长期天气预报和短期气候预测结合，为水库提供专业气象服务。

由图5可以看出，1980—2017年水库入库水量是以2.5万m3/10 a的速率在减少，这与水库上游流域年降水量变化趋势相一致。从表4可以看出，入库水量与上游年均降水量、大雨日数、暴雨日数、雨日强度呈正相关性，且均通过了显著性检验（P＜0.01），入库水量与上游年均气温呈负相关关系但未通过显著性检验。

图5 潘家口水库近55 a入库水量变化

表4 入库水量与相关因子间相关系数及显著性

筛选出2017—2020年（7—9月）上游流域6个国家站的日降水量（20—20时），然后筛选出中雨以上的降水数据，即6个国家站中有一个站日降水量大于10 mm就记为一列有效数据，筛选出同一天的入库水量汇总成数据集。选择2017—2019年的数据进行建模，以2020年的数据来验证模型的准确性。

因为距离的问题，考虑到并不是每个国家站当日的降水量都会对当日的入库水量有所影响，例如距离水库较远的隆化县，如果当日有降水，可能会在2 d或者3 d之后对入库水量有影响，即入库水量和降水量存在滞后性，将水库上游6个站的降水量来构建入库水量预测模型，计算滞后1 d时用当日的降水量作为自变量因子、后一天的入库水量作为因变量，滞后2、3 d计算方法同滞后1 d。采用逐步回归的方法逐步剔除不显著的因子，从而得到入库水量的最优预报回归模型。

利用AIC（赤池信息准则）来比较预测模型，它考虑到了模型的统计拟合度以及用来拟合的参数数目。AIC越小的模型越优先选择。根据表5可得，选择当日降水量与当日入库水量时AIC值最小，因此以此建模较好。

表5 入库水量与降水量滞后性模型AIC值

经过逐步回归分析，程序自动剔除了承德市、滦平、隆化、承德县4个站，得到的水库入库水量与上游流域降水量的模拟方程为：

式中：Q入库量为入库水量；xl，kc分别代表兴隆，宽城国家站的降水量数据。为了验证方程的模拟效果，利用2020年汛期降水日的降水数据，结合同时期的入库水量数据，拟合出了2020年7—9月潘家口水库入库水量模拟值。实测值和模拟值的变化曲线如图6所示，由此可以看出实测值和模拟值曲线的变化趋势基本相同，实测值和模拟值大部分能重合。这说明利用兴隆和宽城的降水数据能较好地模拟出水库入库水量的变化趋势，拟合效果较好。

图6 2020年7—9月入库水量实测值和模拟值变化曲线

为了进一步验证上述回归方程模型的实用性、检验2020年7—9月入库水量实测值和模拟值之间的误差，采用均方根误差（RMSE）、平均绝对误差（MAE）和Theil系数等来检验模型的效果和实测与模拟之间的误差。其计算公式如下：

式中：ai表示实测值；mi表示模拟值。RMSE反映了模拟值偏离实测值的程度，RMSE越小表示模拟精度越高，因此可用RMSE作为评定预测模型精度的标准。Theil不等系数能描述两组数据的偏差，不等系数越小表明实测值和模拟值越接近。

由表6可知，根据入库水量的预测模型得出的均方根误差（RMSE）为206.13、不等系数为0.23。模拟值与实测值偏差不大，方程的模拟效果较好。

表6 2020年7—9月入库水量模拟值与实测值偏差

3 结论

（1）潘家口水库上游流域近55 a年均降水量、大雨、暴雨日数呈现逐年略减少的趋势；雨日强度呈现逐年上升的趋势显著，降水强度的增加也预示着极端降水事件发生的概率在增加；水库上游流域气温升高趋势显著，可能导致蒸发、蒸散的增强，进而导致水库蓄水量的减少。7—8月虽然处于水库主汛期，但面雨量有逐年减少的趋势，主汛期前后的时间段（4—6、9—11月）面雨量降水有逐年增加的趋势，要做好这个时间段的水库入库水量预测工作，为水库安全运行提供气象服务保障。7月为强降水过程发生次数最多的月份，8月集中强降水过程的次数仅次于7月，但强度较7月更大，9月也有集中强降水过程的发生，但次数和强度较7—8月有明显的减弱，进入9月一般来说降雨量会减少，但仍要防范强降水过程对水库蓄水的影响。

（2）本文采用6个国家站的降水量来计算流域面雨量，有研究表明插值雨量站点越多降雨插值结果的精度越高，因此在以后的研究中应增加观测站的数量。本文仅是基于数学统计方法对水库入库水量和上游降水量进行简单的数学建模，但影响入库水量的机理性研究不深入，例如蒸发、蒸散、地表径流等都可能对入库水量产生影响。

（3）潘家口水库上游流域有26个水文站，观测时间序列较长，积累了长序列观测资料，可以研究其气候变化特征。要加快气象、水文部门信息资料共享进程，建立流域信息汇集和共享的机制和渠道，可以通过共同建设实时数据库、资料存储和发布系统，实现流域气象、水文信息共享。