衢州水文站中长期年均流量预报

2018-05-25冯晓蕾何锡君胡孙坚

浙江水利科技 2018年3期

冯晓蕾，何锡君，胡孙坚

（1.浙江水文新技术开发经营公司，浙江杭州 310009；2.浙江省水文局，浙江杭州 310009；3.浙江省工业环保设计研究院有限公司，浙江杭州 310000）

1 问题的提出

水文预报就是根据前期或现时已出现的水文、气象等信息，运用水文学、气象学、水力学的原理和方法，对河流、湖泊等水体未来一定时段内的水文情势做出定量或定性的预报[1]。按其预见期的长短，可分为短期水文预报与中长期水文预报[2]。中长期水文预报的定义为[3]：根据前期水文气象要素，用成因分析与数理统计的方法，对未来较长时间的水文要素进行科学的预测，其预见期一般为3 d ～ l a。

本文以衢州水文站的年平均径流量为预报对象，22项气象环流因子为预报因子，分析对衢州水文站年平均径流量有重要影响的因子。使用逐步回归分析模型，对各年的平均流量分别进行拟合、预报，以期为流域水资源的合理调度提供重要依据。

2 预报模型基本原理

回归分析法是从分析影响预报对象的因子中挑选出一批预报因子，建立其统计规律作为预报的依据。为了预报某对象（如某一水文要素）未来时刻的变化，选择预报量前期已发生的多个有关的水文、气象要素或其他物理要素，即预报因子为研究对象，利用回归方法去分析多个预报因子与这个预报量之间的关系，建立它们统计关系的方程式。从而，对预报对象做出估计，其主要步骤为：① 确定预报量并选择适当的因子；②根据数据计算回归系数标准方程组所包含的有关统计量；③解方程组，定出回归系数；④建立回归方程进行统计显著性检验；⑤利用已出现的因子值代入回归方程做出预报量的估计。

2.1 预报因子分析

利用统计方法来挑选预报因子，统计方法中最常用的就是相关分析，通过计算2个序列y和x之间的线性相关系数，从而判断它们之间的线性相关程度。预报因子与预报对象的相关性计算有很多种方法，本文中采用单相关计算的方法。相关系数的计算公式为：

式中：xi、yi分别为影响因子和待预报水文要素的实测值；x、y 分别为影响因子和待预报水文要素的多年平均值；r为样本的相关系数，其显著程度常用检验来说明：

式中：r为样本的相关系数；n为样本容量（个）。

选定信度α后可从t分布表中查出相应的tα，当t≥tα时，则认为在这一信度下两者是线性相关的，否则认为是不相关的。

2.2 逐步回归模型原理

2.2.1 基本思路

先定义一个衡量因子对预报对象重要性的指标，以便从中挑选出对其影响显著的因子。因子的挑选应逐步进行，在建立回归方程的过程中，每一步只挑选出一个因子，要求当步选出的因子能满足残差平方和下降最多并且能通过指定的信度检验。直至还未引入回归方程的因子中不存在对作用显著的因子为止。随着因子的引入，由于因子之间的相互配合关系，可能产生当后面的因子引入以后会引起前面的因子对的作用显著变小，甚至不显著，要将不显著的因子剔除。因此，在逐步回归中每一步都要做剔除和引进因子的检验。直至方程中既不能引入也不能剔除因子为止。最后的方程中只包含了对预报对象影响显著的因子。而没有进入方程的因子，添加任何一个，都不会对回归效果有显著的改进。

2.2.2 计算步骤

2.2.2.1 建立标准化的正规方程组

假设预报因子X是由p个组成，每个因子的长度为n；预报对象y，其长度为n。在逐步回归中采用的是标准化的正规回归方程组：

其系数与常数项组成的零步增广矩阵为：

记作R(0)= Rij(0)。

2.2.2.2 因子的剔除和引入

假设逐步回归已经进行了l步，回归方程引入了l个预报因子，则第l + 1步的计算步骤如下：

（1）计算全部待选因子的方差贡献：

式中：为第l步时矩阵R(l)中的元素，表示已引入因子的方差贡献，表示未引入因子的方差贡献。

（2）剔除因子。在已引入的l个因子中挑选出方差贡献最小者，记为min{Vk(l)}，其对应的因子为xk，定义方差比为：

给定显著性水平α，在F分布表中查出自由度为1，（n - l - 1）时的上分位点Fα，若Fk≤Fα，则认为因子xk的方差贡献不显著，可剔除；否则不剔除。

（3）引入因子。在为引进方程的因子中挑选出方差贡献最大者，记为max{}，其对应的因子为xk，定义方差比为：

式中：n为样本容量；l为当前因子数。

给定显著性水平α，在F分布表中查出自由度为1，（n - l - 2）时的上分位点Fα，若Fk≥Fα，则引进因子xk，否则不引进。

在逐步回归建模过程中，因子xk的剔除和引进都需要进行消去变换，消去变换的方法为：

式中： i表示矩阵的行数；j表示矩阵的列数，k是被选入的因子的序号。

2.2.2.3 计算结果分析

循环因子引进和剔除步骤，若可供挑选的因子全被引进，或者既无因子可引进又无因子可剔除时，则终止逐步回归分析。

假设最后共引入了m个因子进入回归方程，则应用逐步回归分析法建立的标准化回归方程：

根据标准化前后数据之间的关系，即：

则标准化回归方程可转化为下式：

式中：m为因子数（个）；b0，b1，…bm为回归系数。

3 预报模型的建立

3.1 衢州水文站概况

衢州水文站属国家级重点水文站，是钱塘江上游的一个重要控制站，流域属山地地形，多山溪河流，坡陡流急，容易诱发洪涝灾害。流域内降雨充沛，多年平均降雨量1 815.2 mm，但年内分配极不均匀，降雨主要集中在3 —6月，占全年总量的57%，特别是梅雨季节，常常出现大暴雨、特大暴雨天气。据实际资料统计，流域平均汇流时间为8 h左右。衢州水文站的水文预报，不仅关系到浙赣铁路的安全，还关系到碗窑、白水坑、湖南镇、铜山源、黄坛口等大型水库的错峰调度，在钱塘江流域防洪体系中，处于非常重要的地位。

3.2 衢州水文站流量资料

以衢州水文站年平均流量为预报对象，模型编制依据的水文资料为衢州水文站1951 — 1999年的年平均流量。

根据各参数物理特性，对原始资料及衢州水文站资料进行合理性、一致性分析，数据合理，资料系列完整，无需插补延长，可直接用于计算。由于衢州水文站从1951年起有完整的水文资料，故资料系列取1951 — 1994年共44 a同步资料进行预报分析及拟合检验；1995 — 1999年共5 a同步资料进行预报方程预报检验。

3.3 预报因子挑选

根据前述预报因子分析方法，采用单相关系数法计算衢州水文站年平均流量与22个预报因子之间的相关系数，并对求得的相关系数进行t检验，筛选出与预报对象相关性强的因子。

具体计算时，取衢州水文站1952 — 1994年的长系列年平均流量资料，与某一预报因子1951 — 1993年的系列资料分别计算相关系数，即在提前1 a的范围内进行挑选，在满足r＞rα的基础上，选择相关系数最大的系列为该因子相关月份系列选取；如果预报因子1 — 12月系列与预报对象相关系数都小于rα，则表明该因子与预报对象线性相关不显著，不予引进；依次计算该预报对象与22个预报因子1 — 12月系列的相关系数，即可挑选出相应的预报因子。

由可信度（α = 0.05）及样本数（n = 43），查表得最低相关系数rα = 0.301；由可信度（α = 0.01）及样本数（n = 43），查表得最低相关系数rα = 0.389，具体预报因子挑选结果见表1。

表1 预报因子及相关系数表

3.4 建立预报模型

首先，依次挑选出该预报对象与22个预报因子系列中相关系数最大的系列，再利用预报因子及衢州水文站年平均流量实测资料，根据上文逐步回归模型原理，通过程序计算方程的各项系数。

逐步回归每一步的回归方程系数，过程一共进行了11步，从最后第11步的计算结果可知：22个因子最终有9个因子进入方程中，分别是因子7、因子2、因子11、因子4、因子10、因子19、因子21、因子20及因子16。“非标准化回归系数”中的“B”列系数代入，得衢州水文站年平均流量的预报方程：

y = 1 049.377 + 1.221X1- 0.727X2- 2.117X3- 4.439X4+0.451X5+ 1.618X6- 17.557X7+ 5.894X8- 4.013X9（11）

式中：y为衢州水文站年平均流量（m3/s）；X1为上一年 12 月 H 500 mb（25 ～ 35 N、110 ～ 130 E）长江中下游区7点合计；X2为上一年3月 H 500 mb（50 ～ 55 N、70 ～ 90 E + 40 ～ 45 N + 65 ～ 85 E）巴尔喀什湖区；X3为上一年10月H 500 mb（30 ～ 40 N、80 ～ 90 E）西安高原子6点合计；X4为上一年1月H 500 mb（120 E、20 ～ 40 N）高度差（沿120 E线20 ～ 40 N）；X5为上一年4月H 500 mb（50 ～ 60 N、100 ～ 120 E）贝加尔湖区8点合计；X6为上一年10月鄂海平均高度H（135 ～ 150 E、45 ～ 60 N）；X7为上一年1月C102 102站西风指数（115 E、25 ～ 30 N）；X8为上一年1月西风风速V27.5 N（m/s）（105 E、25 ～ 35 N）；X9为上一年8月西风风速V42.5 N（m/s）（105 E、40.5 N）。

3.5 模型效果检验

（1）复相关系数（R）。复相关系数为：

式中：R为复相关系数；U为回归平方和；Q为残差平方和；Syy为总和平方和。

（2）剩余标准差（Sy）。剩余标准差为：

式中：Sy为剩余标准差；n为资料年限（a）；Q为残差平方和； m为挑选的因子（个）。

（3）回归效果的F检验。回归效果的F检验为：

由计算程序求得上述回归方程的复相关系数R、剩余标准差Sy和F检验值。

根据本例样本数n = 43，因子个数m = 9，在信度α =0.01的条件下，查复相关系数表有临界值Rα= 0.671，而衢州水文站逐步回归方程的R = 0.852，R＞Rα。因此在α =0.01的水平下回归的效果是显著的。该方程剩余标准差Sy=34.04 m3/s。

根据本例fU = 9，fQ = 33，在信度α = 0.01的条件下，查F分布表得临界值Fα= 3.01，而衢州水文站逐步回归方程的F = 9.72，F＞Fα。因此在α = 0.01的水平下回归效果是显著的。

4 模型方程拟合及预报结果

4.1 模型方程的拟合

根据已建立的衢州水文站年均流量预报模型所用的1952 — 1994年共43 a的实测年平均流量对预报模型进行历史拟合检验，计算相应年份的预报误差及许可误差，其中许可误差采用实测流量的20%来计算，如预报误差＜许可误差，则该年预报值合格，否则为不合格。该预报方程在43 a的平均流量历史拟合中预报合格次数为37次，合格率为86.0%，根据GB 22482 — 2008 — T《水文情报预报规范》规定，中长期预报合格率≥85%时，预报方程等级为甲等。具体拟合结果见图1和表2。