最大降水量多年一遇计算方法及Matlab实现

N=length(li); M=1:1:N; pm=1:1:1000; x=li; temp=sort(x); [p,pci]=gamfit(temp); ffit=gamcdf(temp,p(1),p(2)); ave_x = mean(x); s_x =std(x); y=-log(-log(1-M/(N+1))); ave_y=mean(y);
文 章 编 号 :1007-0907(2012)01-0075-02
当暴雨洪涝等灾情出现时，“几年一遇” 经常作为标志灾害
强度出现在媒体口中。 那么如何计算和快速计算某一特定气象
要素值究竟是几年一遇呢？ 本文通过 Matlab 实例对比分析三种
常见的多年一遇的计算方法。 “几年一遇”在数学上实际上就是
概率与数理统计中极值分布问题。这类问题的解决，在城市建设
灾害防御规划中，有重要的社会和经济意义。
本 文 使 用 的 Matlab 版 本 为 7.10.0 (R2010a)， 以 阿 里 河 镇
1961-2009 年 49 年最大日降水量为例，分别通过指数分布法、耿
贝尔分布法及皮尔逊 III 分布法讨论及编程计算。 具体数据如下：
内蒙古农业科技 2012（1）：75~76 Inner Mongolia Agricultural Science And Technology
最大降水量多年一遇计算方法及 Matlab 实现
王 涵， 李 玲 （鄂伦春旗气象局，内蒙古 鄂伦春 165450）
中图分类号:S161.6 文献标识码:A 10.3969/j.jssn.1007-0907.2012.01.042
型的，指数分布形式如下：
x=algT+b
根据求解过程，关键部分编程如下：
N=length(li);
laver=mean(li);
for i=1:N
Ti(i)=(N+1)/i;
end
Taver=mean(Ti);
for i= 1:N
收 稿 日 期 :2011-12-28
a1=a1+abs(li(i)-laver)*abs(log(Ti(i))-log(Taver)); a2=a2+abs(log(Ti(i))-log(Taver))^2; end a=a1/a2; b=laver-a*log(Taver); latT=a*log(T)+b; end
ylabel(' 降水量(mm)'); title(' 阿里河镇日最大降水量多年一遇极值曲线 ');
先通过上文中三种分布 法 计 算 函 数 得 出 Y1、Y2、Y3， 并 运 行该绘图函数，得到以下对比曲线（图 4）。
需要指出的是， 无论用经验分布函数还是极值的理论分布 函数研究实际极值问题， 其效果优劣主要取决于所采用的分布 函数对实际资料的拟合程度， 本文没有通过实际数据对预测结 果进行拟合。 另外，要使研究的问题获得满意的结果，所采用的 资料样本数必须足够多， 各次观测必须相互独立且服从同一分 布。 对于不同的气象要素， 还有其他一些分布方法， 比如韦伯 （Weibull）分布、对数正态分布 等 ，本 文 没 有 一 一 测 试 ，对 于 气 候 极值的观测而言， 这些分布函数一般都能满足相互独立和服从 同一分布的要求。 参考文献： [1] 王 丽 文.珠 海 市 年 雨 量 和 年 最 大 日 雨 量 多 年 一 遇 的 极 值 计 算[J].广
东 气 象 ,2006,(1):25-28. [2] 杨 娟. 贵州年降水量和年最大月降水量多年一遇的 极 值 计 算[J].
贵 州 气 象 ,2008,(6):10-12. [3] 赵佩红，林国生，聂燕红.新会年和月最大日降水量多年一遇的极值
计 算 [J].气 象 ,2010,(2):88-89. [4] 李明杰，齐 鹏 ，侯 一 筠.山 东 沿 岸 多 年 一 遇 最 高 水 位 计 算[J].气 象,
分布变量的分布计算。 皮尔逊 III 分布形式如下：
軃
軃軃α=[
軃 軃
2 cs
軃
],β=[
2 x軃·cv·cs
]
軃
軃
軃軃軃y=x軃·（1-2
軃 軃
cv cs
）
其 中 ，x 为 样 本 ，y 为 重 现 期 值 ，cs 为 偏 态 系 数 ，cv 为 变 差 系 数，α 为形状参数，β 为尺度参数。 根据求解过程，关键部分编程
通过计算重现期为 1～1000 的重现期值得到如下曲线（图 1）。
2 耿贝尔分布法 耿 贝 尔 (Gumbel)分 布 是 一 个 适 用 较 广 的 分 布 理 论 , 同 样 也 属
于指数型分布。 耿贝尔分布形式如下： x=u-a*lgT(x) 其中 u 表示一个与降水序列有关的函数。 根据求解过程，关键部分编程如下：
latT=gaminv((1-1/T),Alpha,1/Beta)+X0;
end
通过计算重现期为 1～1000 的重现期值得到如下曲线（图 3）。
4 绘图及结论
编制绘图函数如下：
function createfigure(Y1,Y2,Y3)
figure1 = figure;
axes1 = axes('Parent',figure1);
如下：
N=length(li);
laver=mean(li);
lstdevp=std(li);
lskew=skewness(li);
Cs=lskew*(N-1)*(N-2)/N^2;
Cv=lstdevp/laver;
Alpha=4/Cs^2;
Beta=2/(Istdevp*Cs);
X0=laver*(1-2*Cv/Cs);
8.2,33.2,63.8,53.6,30.2。
需要说明的是，本文程序使用的数据为列矩阵，所以如果输
入行矩阵，需要做转置操作。
极值统计的根本目的是准确的推断极值序列的重现期值或
某一极值平均可能在多少年内出现一的事件，平均在 T 年内出现 1 次时，
76
内蒙古农业科技
No.1
s_y=std(y); a = s_x/ s_y; u = ave_x -s_x/s_y*ave_y; xp_Genbel = u - a*log(-log(1-1./pm) );
通过计算重现期为 1～1000 的重现期值得到如下曲线（图 2）。
3 皮尔逊Ⅲ分布法
皮 尔 逊Ⅲ（Pearson-III）原 理 一 般 适 用 于 统 计 降 水 量 等 左 偏
60.8,34,40,36.6,38.8,52.2,48.4,53.2,73.2,51.6,37.3,74.3,43.2,51,40.7,
36.7,49,55.4,70.7,42.1,106.8,72.9,32.2,83.1,47.5,58,46.3,29.6,48.8,7
1.2,66,29.7,53.2,51.9,47.7,43.4,38,26.6,58.2,37.7,58.6,34.8,42,31.8,3
2009,(11):78-81. [5] 王 涛，朱敏嘉，崔爱萍，等.山西近 45 年降水变化对旱涝的影响分
析[J].山 西 农 业 科 学 ，2010，38（12）：61-63. [6] 王 健，孙 明，薛明霞，等.井 灌 区 降 水 量 的 特 征 分 析 及 灰 色 预 测
[J].山 西 农 业 科 学 ，2010，37（11）：51-56. （责任编辑 侯旭光）
box(axes1,'on');
hold(axes1,'all');
plot(Y1,'b',Y2,' -.',Y3,' r:');
legend (' 指 数 分 布 法 ',' 皮 尔 逊 III 分 布 法 ',' 耿 贝 尔 分 布 法 ',
'Location','SouthEast');
xlabel(' 重现期(年)');
1-F(x)
F(x)
这就是说，重现期就是右侧概率或左侧概率的倒数。而重现期内
的
极大
值
或极
小
值可
分
别由
F(x)=1-
1 T(x)
和
F(x)=
1 T(x)
解
出
x。
即
重现期和重现期值互为函数关系，可以互相解出，只要知道函数
F(x)即 可 。
1 指数分布法
对于极大值问题而言， 很多气候要素的极值分布都属指数
则把这个 T 叫做 X 的特定值 x 的重现期，而在 T 年内平均出现
1 次的这个特定值 x 叫做重现期值。 如果令 F(x)为 X 的 分 布 函
数，x 代表 X 的年最大值或最小值 ，则根据概率分布 ，最 大 值 或
最 小 值 的 再 现 期 T(x)可 分 布 由 ：T(x)= 1 和 T(x)= 1 给 出 。