洪峰流量预测

承 诺 书

我们仔细阅读了数学建模选拔赛的规则.

我们完全明白，在做题期间不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人研究、讨论与选拔题有关的问题。

我们知道，抄袭别人的成果是违反选拔规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守选拔规则，以保证选拔的公正、公平性。如有违反选拔规则的行为，我们将受到严肃处理。

我们选择的题号是（从A/B/C中选择一项填写）： A 队员签名 ：

日期：2013年8月23日

2013年河南科技大学数学建模竞赛选拔

编 号 专 用 页

评阅人 评分 备注 评阅编号（评阅前进行编号）：

评阅记录（评阅时使用）：

洪峰流量预测

摘要

本文讨论了关于百年一遇的定义、百年中最大洪峰流量预测、怎样根据现有的年最大洪峰样本来预测未来三年的洪峰值等问题。

首先对于问题一，通过查阅资料并客观的分析，可认为“百年一遇”不能简单机械的理解为“一百年只出现一次”。百年一遇应该理解为数学统计中的“概率”，指的是随机事件发生的概率，表现在时间上就是我们平常所说的“时间周期”。百年一遇，就是指超过一定标准的洪峰流量超过一定标准发生的概率是1%。

对于问题二，考虑到洪峰流量的分布，在较大值和较小值分布较少，处于一般情况下的分布最多，因此，考虑正态分布建立模型。我们通过spss进行正态分布拟合的探索。然后我们对拟合的正态分布的正态性进行检验，结果发现正态性良好，我们建立的正态分布模型成立。然后我们求得此正态分布概率密度函数，然后对正态分布概率积分函数进行求解，得到概率积分为0.99的临界点x0781.14，作为“百年一遇”的观测值，即超过此值，我们就称为“百年一遇”。此观测点x0781.14为百年一遇的洪峰流量值。 对于问题三，题中要求对后三年的洪峰量进行预测，再分析了题目后我们选择了用灰色理论对数据进行预测，建立灰色模型求解预测结果三年内结果均无明显变化，对模型检验发现模型并不合理，后来我们又对模型进行了改进，用马尔科夫与灰色理论混合模型能够比较好的通过检验并预测了后三年的洪峰量为： 年份 2013 2014 2015 洪峰值 583.67 516.67 513.95

关键字:正态分布 马尔科夫预测 灰色预测

1

问题重述

外界对三峡工程“万年一遇”“千年一遇”“百年一遇”等防洪标准说法不一提出质疑，并将相关报道整理如下：

现假设附件中是某水文站的每年的最大洪峰流量观测值。

1、有人说，百年一遇就是一百年内只出现过一次，能否这样理解？给出你对“百年一遇”的定义。

2、 能否计算出附件中水文站百年一遇的洪峰流量值；如能，给出具体的计算过程及结果。

3、 预测该水文站将来3年内每一年的最大洪峰流量。

2

问题分析

对于问题一，如何定义“百年一遇”，首先否定了“百年一遇就是一百年内只出现过一次”的理解。经过查相关资料所得：所谓“百年一遇”这是一个关于频率的概念。它表示在许多次试验中某一事件重复出现的时间间隔的平均数。需要特别指出的是所谓“重现期”并不是说正好多少年中出现一次，它带有统计平均的意义，说得更确切一点是表示某种水文变量大于或等于某一指定值，每出现一次平均所需的时间间隔数。 对此，可以通过随机模拟来对此进行解释，通过MATLAB随机产生一百个随机数来代表这一百年的洪峰，用其中的最大洪峰值来代表百年一遇的最大洪峰值，再用随机抽样的方法确定最大洪峰值发生的概率。从而可以确定上述说法的错误在于把平均发生的次数认为了为一定会发生。准确定义应用概率来定义为：在一百年中每天最大洪峰发生的概率为1%。

对于问题二求最大洪峰值,可以根据现有的某水文站的数据，将其分成间距相同的组，计算出各个组中的频数，通过计算出这些样本的平均值和方差来估计整体的均值和方差，从而确定洪峰值分布的概率分布曲线。通过确定洪峰值大于等于百年洪峰值得概率为1%所对应于的洪峰值作为百年一遇的最大洪峰值。

对于问题三，可以采用灰色理论预测模型，预测该水文站将来3年内每一年的最大洪峰流量。

基本假设

1、所给数据真实有效，水文站的检测结果可靠无差。

2、每年最大洪峰流量只考虑时间的影响，附表没给出影响因素，不加考虑。 3、不考虑环境剧烈变化所引起的对预测结果的影响。

符号说明

Xi第i个洪峰流量值Fx洪峰流量的概率密度函数Xi的均值Xi的标准差g1偏度g2峰度1

3

g偏度的标准误g峰度的标准误2

ttx

带估参数

x0概率积分为0.99的观测点 Xi第i年份的最大洪峰值

0Xi第i个叠加序列的值

1模型的建立与求解

对于问题一：

经查阅相关资料可以得到关于百年一遇的理解: 所谓“百年一遇”这是一个关于频率的概念。往往用“重现期”来替代“频率”，它表示在许多次试验中某一事件重复出现的时间间隔的平均数。需要特别指出的是所谓“重现期”并不是说正好多少年中出现一次，它带有统计平均的意义，说得更确切一点是表示某种水文变量大于或等于某一指定值，每出现一次平均所需的时间间隔数。经分析可以用随机模拟的方法来解释，先用MATLAB随机产生一百个数如附录：

用随机抽样的方法对数据进行单次抽样，进行一百次，分别代表一百年，每个数分别代表每一年的洪峰值，最大的数可以作为百年一遇的最大洪峰值。用excel进行一百次的随机抽样并统计不同阶段的的概率得：

流0-100 100-200 200-300 300-400 400-500 500-600 600-700 700-800 800-900 量 频0.01 0.02 0.15 0.43 0.31 0.07 0.01 0 0 率 并得出其直方图：

4

可以得出最大流量值得概率为：P=0.01。则发生次数均值为EX1，百年中发生

次数为1仅仅是在大量实测统计的情况下，百年发生的均值。但是在实际中每一年都有0.01的概率发生。所以不能说是一百年中一定会发生。

我们期望找到一个的洪峰流量值，作为一个观测点，即衡量的标准，超过此洪峰流量，我们就称之为“百年一遇”。此观测点，应能保证有99%每年最大洪峰流量都小于此观测点，即大于此观测点的每年最大洪峰流量出现的概率为1%，引入重现期，即重现期为100年。 对于问题二：

考虑到每年的洪峰流量的值，分布的范围在，很大或很小的情况很少，而大多数都分布在洪峰流量既不很大，也不很小的范围内，由此我们联系到所学的正态分布，可以考虑建立起一个关于洪峰流量的正态分布的模型。

1、模型的建立

建立以下正态分布的模型：

Fx1e2(x)222

2、模型求解 Xi X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 洪峰638658.流.10 99 量 XiX12 X13 577.11 502.42 384.55 455.04 405.57 555.06 510.22 578.33 470.90 X14 X15 X16 X17 X18 X19 X20 X21 X22 洪峰834.流00 量

703.92 597.79 521.35 535.82 348.53 463.45 502.04 597.32 472.09 433.74 Xi X23 X24 X25 X26 X27 X28 最大洪610.67 512.23 465.35 555.56 485.52 672.23 峰流量 通过spss中分析探索描述性统计，可得到正态分布的，。

5

以下是spss处理的结果：

可得到=537.43，=104.597。

正态性检验：经查询资料，得知，正态分布的正态性的两个评判标准为偏度g1和峰度g2，以及Q-Q图。

（1）、偏度，峰度检验方法

通过计算g1和g2及其标准误g1及g2然后作U检验。两种检验同时得出

UU0.051.96，即在95%的置信概率下，可以认为正态性良好。

Ug1g

1Ug2g

2下面给出spss做出的偏度、峰度检验的结果 ：

6

描述 均值 均值的 95% 置信区间 下限 上限 统计量 496.867 577.983 532.936 516.790 10940.440 104.5966 348.5 834.0 485.5 130.9 .733 1.122 标准误 .441 .858 537.425 19.7669 5% 修整均值 中值 方差 最大洪峰流量 标准差 极小值 极大值 范围 四分位距 偏度 峰度

可得g10.733 g10.441；g21.122 g20.858 计算得：偏度检验：Ug10.7331.66211.96 0.4411.1221.30771.96 0.858g1 峰度检验：Ug2g2由检验结果可知，在95%的置信概率下，偏度、峰度检验均合格，即在偏度、峰度检验下，正态分布的正态性良好。

（2）、 Q-Q图检验

以样本的分位数作为横坐标，以按照正态分布计算的相应分位点作为纵坐标，把样本表现为指教坐标系的散点。如果资料服从正态分布，则样本点应该呈一条围绕第一象限对角线的直线。

我们根据spss作出Q-Q图如下：

7

由图可知，28个观测点中，由27个都是很好的围绕（几乎完全分布）在第一象限对角线上，只有第28个观测点，偏离稍远。因此，我们完全可以认为Q-Q图检验的结果为，其探索的正态分布的正态性良好。

（4）、正态分布模型

Fx1e2(x)222

=537.43

=104.597

即：

N,2x0

我们此时是由样本来估计总体，由样本求出的=537.43，均是对总体，=104.597，

的无偏估计。故我们建立的这一正态分布模型可以推广到总体更一般的情形。

则，总体最大洪峰流量的概率密度函数为：

Fx1e104.5972(x535.43)22(104.597)2

下面是我们用matlab做出的总体正态分布曲线图

8

3、对“百年一遇”的最大洪峰流量值得计算 （1）、Px概率积分函数 t概率积分函数

PxFxdx

引入新的变量t，tx

对于一般的正态分布N,2,Fx的值可以利用以下变换公式计算：

uu11x222Fxedtedu 22注：此变换公式为摘录，其中变量与我们建立模型所用变量含义可能不同，仅供参考。

故我们做类似变换得到：

x(t)2tx2Pxt

不同的t值，可由正态分布积分表查出。

由我们在第一问给出的定义，“百年一遇”即大于此洪峰流量的概率小于1%，即“百

9

年一遇”的临界值x0，小于此洪峰流量值的概率积分为99%

Px99% t99%

查表得出t2.33，得出了“百年一遇”的最大洪峰流量值临界点x0

x0t

代入数据：104.975 537.43 得：

x02.33104.597537.43781.14

由此，我们得出了总体的“百年一遇”的洪峰流量值的临界点x0781.14，即当某

年的最大洪峰流量大于此临界点，我们就称为百年一遇。

对于问题三

（1）模型建立

通过对问题的分析,得出对洪峰值的预测可以建立灰色理论预测模型。建立某水文站从1985到2012年最大洪峰值的原始时间序列为：

X(0)X(0)(1),X(0)(2),X(0)(3),...,X(0)(28)=638.10，658.99，，485.52，672.23

将原始数据按如下方式累加：

X1X11=638.1010010X2X1X2X1X2=1297.09XX10

27i12827X0iX126X027=14375.6728X0iX127+X028=15047.9i1111得到新的序列：

XX11,X2,X3,...,X28=638.10，1297.09，，14375.67,15047.9

11令Z1为X1的紧邻均值生成序列：

Zk1111Xk1Xk 210

Z=Z2，Z3,,Z28967.595,1585.645,,14711.785

1111则GM(1,1)的灰微分方程为：

X0kZk

1其中：称为发展系数；称为内生控制灰数。

设为带估参数，



可利用最小二乘法求解。解得：

BTBBTY 其中

1Z121X0201Z31X3 BY

Z1281X028用MATLAB解得:0.0001401315909 534.7959806170588 可得:

01Xk1X1ek 带入参数得:

534.79598061705880.0001401315909k534.79598061705881Xk1638.10e 0.00014013159090.0001401315909还原值可得：X0k1Xk1X1k

X1则可通过MATLAB求得2013、2014、2015年最大洪峰值:

2013532.57530X2014532.5007 0X2015532.4261对灰色模型的检验： 年份 1985 1986 1987

0误差 0.00% 23.25% 7.95% 11

年份 1999 2000 2001 误差 2.31% 0.41% 34.69% 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 6.01% 28.05% 14.85% 24.09% 3.9% 4.48% 8.29% 11.82% 56.20% 31.86% 11.99% 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 13.13% 5.88% 12.00% 11.47% 18.65% 14.55% 3.9% 12.68% 4.26% 8.87% 26.19%

（2）模型反思与改进

我们通过灰色理论预测2013,2014,2015的最大洪峰流量：

X2013532.57530X2014532.5007 0X2015532.4261发现预测的结果，每年最大洪峰流量的值几乎不变，这与我们实际每年最大洪峰流量的波动不符，并且误差较大，所以我们决定对这一模型，进行改进。

改进思路：考虑到每一年的最大洪峰流量和下一年的最大洪峰流量之间有较大的相关性，所以我们采用马尔科夫状态转移矩阵来改进。

首先进行数据处理，把洪峰流量进行一个分级： 800-900 一级；700-800 二级；600-700 三级； 500-600 四级；400-500 五级；300-400 六级。

我们可以由附表得到每一年的最大洪峰流量等级和下一年的最大洪峰流量等级之间存在一定的联系，所以建立由该水文站上一年的最大洪峰值转移到下年最大洪峰值时转移概率矩阵如下：

次年 一级 二级 三级 四级 五级 六级 当年 一级 0 1 0 0 0 0 二级 0 0 0 1 0 0 三级 0 0 1/3 2/3 0 0 四级 0 0 0 1/2 1/3 1/6 五级 1/8 0 1/4 3/8 1/4 0 六级 0 0 0 0 1 0 012

01000000010012000033因此可得到马尔科夫状态转移矩阵为：R111

0002361131008484000010

由附表我们知道：

2012年最大洪峰流量值为672,23，属于三级，故对应的2013年最大洪峰流量值

12出现在三级的概率为，出现在四级的概率为，出现在其他等级的概率为0。

3312故预测2013年最大洪峰流量值为650550583.33（650为第三等级范围

33中值，550为第四范围等级的中值）

2013年预测的最大洪峰流量值为583.33，属于四级。 由MATLAB可得：

0001000000.50.33330.1667000.11110.55560.22220.11112R

0.041700.08330.3750.41670.08330.03130.1250.14580.44790.18750.06250.12500.250.3750.250故对应的2014年最大洪峰流量值出现在三级的概率为0.1111，出现在四级的概率为0.5556，出现在五级的概率为0.2222，出现在六级的概率为0.1111。

故预测2014年最大洪峰流量值为5500.555645670

0000.50.33330.16670.041700.08330.3750.41670.08330.027800.09260.43520.35190.09263R

0.05210.04170.13190.39930.31250.06250.02340.03130.09550.51650.25870.07470.03130.1250.14580.44790.18750.0625故对应的2015年最大洪峰流量值出现在一级的概率为0.0278，出现在三级的概率为0.096，出现在四级的概率为0.4352，出现在五级的概率为0.3519，出现在六级的概率为0.0926。

故预测2015年最大洪峰流量值为

8500.02786500.09265500.43524500.35193500.0926513.95

13

模型评价

模型优点：

1、本文对百年一遇的解释巧妙地利用了随机模拟的模型，充分的解释了百年一遇是一种关于概率的概念。

2、本文对问题二最大洪峰值得预测充分的利用了正态分布曲线的概率模型，巧妙地利用概率统计的方法得出最大洪峰值。 模型的不足：

1、对于用灰色模型来预测未来三年的洪峰值，经过验证并不理想。 3、对洪峰值得等级分配会带来一定的误差。

参考文献

[1]、党耀国,灰色预测与决策模型研究,,科学出版社,2009-12-01 [2]、彭放 杨瑞琰 罗文强 肖海军 何水明,数学建模方法,科学出版社 2007 [3]、费业泰,误差理论与数据处理,机械工业出版社,2010-6

14

附录

随机数：

625.33 591.82 313.76 506.70 841.59 444.54 535.33 449.61 448.73 623.74 547.89 533.79 679.11 420.26 681.66 480.47 453.94 425.28 801.64 426.74 569.17 643.98 637.94 710.59 488.41 474.63 523.49 550.40 569.59 508.93 398.10 293.73 385.85 572.31 578.36 584.67 500.41 451.28 372.47 590.56 566.92 540.93 395.00 585.02 448.65 502.40 595.24 646.11 539.82 510.02 354.36 507.55 450.47 435.00 517.01 514.66 505.72 539.84 542.79 623.83 588.67 462.69 687.70 614.76 678.78 332.32 716.48 513.91 516.75 517.12 475.81 411.09 用excel抽样数表：841.5939 638.2805 586.2628 710.5851 627.5899 585.0197 710.5851 627.5899 585.0197 697.1422 623.8287 584.6693 697.1422 620.2047 578.3594 681.6613 614.7602 578.3594 678.7759 610.8041 572.3092 678.7759 602.818 572.3092 646.1108 588.6739 572.3092 646.1108 586.2628 569.5923 统计表：

MATLAB程序：

B=[-967.595

-1585.645 -2125.41 -2568.895 -2988.69 -3418.995 -3899.31 -4431.95 -4976.225

405.93 652.32 446.90 508.36 518.96 610.80 620.20 322.81

569.5923 539.8355 569.5923 539.8355 556.6383 539.8355 556.6383 539.8192 556.6383 530.3265 547.8944 530.3265 542.7898 518.9601 540.9269 518.9601 540.9269 517.1203 540.9269 517.1203

流量 频率800-900

700-800 600-700 500-600 400-500 300-400 200-300 100-200 0-100 523.80 556.64 487.62 627.59 517.1203517.0054514.664514.664508.9255508.3604507.5468507.5468506.6993506.14360.01 0.02 0.15 0.43 0.31 0.07 0.01 0 01

1 1 1 1 1 1 1 1

15

397.93 655.36 574.05 506.14 506.1436481.6165481.6165480.4689475.8143474.6326474.6326453.9395453.9395451.2803420.52 669.29 606.47 530.33 451.2803451.2803451.2803450.4666450.4666449.6065448.7345448.7345448.7345448.6525416.47 481.62 327.96 638.28 446.8996444.5395444.5395444.5395435.003426.7425426.7425425.2808425.2808425.2808697.14 586.26 515.49 602.82

420.5233 405.9342 332.3207 332.3207 332.3207

322.8112 322.8112 322.8112 313.7595 293.7279

-5500.84 -6153.29 -6922.25 -7573.105 -8132.675 -8661.26 -9103.435 -9509.425 -9992.17 -10541.85 -11076.555 -11529.47 -12051.675 -12613.125 -13101.915 -13612.37 -14132.91 -14711.785 >> yN=[658.99 577.11 502.42 384.55 455.04 405.57 555.06 510.22 578.33 470.9 834 703.92 597.79 521.35 535.82 348.53 463.45 502.04 597.32 472.09 433.74 610.67 512.23 465.35 555.56 485.52

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];

16

672.23];

>> a=inv(B'*B)*B'*yN a =

0.0001 534.7960

>> digits(8) >> vpa(a) ans =

0.00014013159 534.79598

17