雨量预测方法

承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

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

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。

如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛的题目是：C题雨量预报方法的评价
我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）：江西师范大学科学技术学院
参赛队员(打印并签名)：1.熊军军
2.许谞
3.许盛敏
指导教师或指导教师组负责人(打印并签名)：温利民
日期：2005年09月19日赛区评阅编号（由赛区组委会评阅前进行编号）：
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：
雨量预报方法的评价模型
摘要
本文建立了一个关于雨量预报方法的评估模型。

首先，通过对给定的大量数据（预报数据和实测数据）进行统计画图分析，得出了散点图。

然后分别对两种不同方法预报的41天中每天4个时段各等距网格点的雨量数据进行处理和分析。

在可接受的度数差范围内搜索与各个观测站点距离最近的网格点，按从小到大排序后取其最小的4个网格点，再根据欧氏距离倒数加权的方法对它们赋权重，取出4个网格点对应的雨量，分别与各自的权重相乘，累加得到的值来预测相对应观测站点的雨量。

对得到的观测站点的预测雨量进行两种方法的分析，方法一：将预测雨量与实测雨量求偏差率，并对所有偏差率求出一个偏差率的算术平方根，作为评价准确性的指数，从而得到第一种雨量预报方法的准确性的指数为102.8755，第二种雨量预报方法的准确性的指数为726.6841；方法二：将预测雨量与实测雨量分别转化为对应的级别（如雨量在区间0.1——2.5为1级），用同级率比较法将它们作比较，从而得到第一种雨量预报方法的同级率为73.9346%，第二种雨量预报方法的同级率为70.9662%。

本文利用数学软件Matlab很好地实现了编程模拟计算，并结合实际测得的数据得出了雨量预报方法的同级率，很好地指导了人们的生活与工作。

关键词：（预报、实测、网格点、同级率）
一、问题的重述与分析
1、问题的重述
随着气象事业现代化建设的快速发展，雨量预报对指导农业生产和城市工作和生活有重要作用，但如何准确、及时地对雨量作出预报是一个十分困难的问题，近年来，随着社会经济的不断发展，预报方法对于提高气象服务水平，增强防灾减灾能力具有重要意义，因此，广受世界各国关注。

我国某地气象台和气象研究所正在研究6小时雨量预报方法，即每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度、北纬32度附近的53×47的等距网格点上。

同时设立91个分布不均匀的观测站点实测这些时段的实际雨量。

气象部门提供了41天的用两种不同方法的预报数据和相应的实测数据（预报数据在文件夹FORECAST中，实测数据在文件夹MEASURING中）。

现在我们所关心的问题就是：
（1）对气象部门提供的大量数据（预报数据和实测数据），怎样进行合理、有效地分析，进而建立数学模型，来评价两种6小时雨量预报方法的准确性；
（2）气象部门将6小时降雨量分为6等：0.1—2.5毫米为小雨，2.6—6毫米为中雨，
6.1—12毫米为大雨，12.1—25毫米为暴雨，25.1—60毫米为大暴雨，大于60.1
毫米为特大暴雨。

所以，若按此分级向公众预报，如何在评价方法中考虑公
众的感受？
2、问题的分析
我们从题目中了解分析到：气象台每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度、北纬32度附近的53×47的2491个等距网格点上。

同时设立91个分布不均匀的观测站点实测这些时段的实际雨量。

由于网格点比较多，且每个网格点的位置是以经度和纬度表示处在一定的区域，所以我们把纬度看作x轴，经度看作y轴，采用Matlab图形处理功能的基本绘图命令plot 画出散点图（图一），程序见附录一。

从图中可以分析看出，气象部门提供了在2491个网格点上41天4个时间段的大量预报数据（雨量），并且同样给出了91个观测站点的实测数据（雨量）。

所以我们想通过网格点上的预报数据来预测实测站点的数据。

然而，观测站点集中在所有网格点的中央部分，而四周是大量的距离比较远的网格点。

因此，通过搜索出2491个网格点中对站点影响比较大的几个网格点，再用搜索出来的几个网格点的预测数据加权求出一个预测数据（雨量），进而和该站点实测数据进行比较，来评价两种6小时雨量预报方法的准确性。

在向公众预报时，采取一种合理、准确的预测方法，增加雨量分等级预报的同级率，能对公众起到良好的出行指导作用，使人们对雨量预报有更深的理解，更多的关注。

二、模型的基本假设和符号说明
1、模型假设
（1）
观测站点的设置是不均匀的；（2）
题中网格是等距的正方形网格（所谓“正方形网格”是指每个格子都是正方形的网格；网络点是指纵线和横线的交叉点）；（3）一个x 轴、y 轴分别为纬度和经度的坐标，通过把点的纬度和经度分别看作
横坐标和纵坐标，用欧氏距离计算公式来作为两点
22)()(j i j i b b a a d −+−=之间的距离。

（4）
点到观测站点的距离越短，则对观测站点的雨量影响越大；（5）
单个网格点到观测站点距离倒数与所取的4个网格点到观测站点倒数之和的比为它的权值；（6）雨量用毫米做单位，小于0.1毫米视为无雨；
2、符号说明
网格点及其对应的纬度和经度:
)2491,,2,1)(,(⋯=i n m X i i i 个第i 观测站点及其对应的纬度和经度
:)91,,2,1)(,(⋯=i b a P i i i 个第i 可接受度数差:
ε
与观测站点的距离最小的前4个网格点的对
:)4,3,2,1,91,,2,1()(===×n i q Q n i in ⋯个第i 应权矩阵与观测站点的距离最小的前4个网格点的距
:)4,3,2,1,91,,2,1()(===×n i d D n i in ⋯个第i 离矩阵欧氏距离的计算公式
22)()(j i j i y y x x d −+−=时段里，与观测站点的距离最小的前4)
4,3,2,1,91,...,2,1()(===×n i f Y n i in j 个第j 个第i 个网格点的预测降雨量矩阵
观测站点的时段预测降雨量矩阵)164,...,2,1,91,...,2,1()(===×j i y M j i ij 个第i 个第j 观测站点的时段实测降雨量矩阵)164,...,2,1,91,...,2,1()(===×j i s I j i ij 个第i 个第j 预报偏差率:
σ预报偏差率的算术平方根（准确性指数）:
S 统计方法的预报数据与实测数据处在同一级6210,...,,i i i i J J J J 种第i 别、相差1级、相差2级、、、、相差6级的频数
三、模型的建立及求解
一、问题（1）及其求解
算法算法：：
1．
根据题意，气象部门提供了41天用两种不同方法的预报数据和相应的实测数据，每种预报方法都有大量的预测数据。

为了评价两种6小时雨量预报方法的准确性，我们采用网格点上的预报数据来预测观测站点的数据，再来和实际测得的数据相比，判断其准确性。

以坐标为基准点，给定一个可接受度数差（在求解中取),(i i i b a p )91,,2,1(⋯=i ε，可搜索得到9至19个网格点），对任意的，搜索其任一个25.0=ε)91,,2,1(⋯=i p i 观测站点在纬度和经度都上下增加的正方形内的所有等距网格点。

若网格点ε()n m X ,纬度和经度在同时满足和时，即认为该网格点是可接εε+≤≤−i i a m a εε+≤≤−i i b n b 受范围内的网格点。

2．
找到可接受范围内的网格点后，我们计算网格点和这些观测站点间的距离。

再得到各观测站点和等距网格点之间的距离后，将各观测站点按距离从小到大排序后保存，我们用Matlab 编程求得结果（见附录二）。

3．
各观测站点和等距网格点之间的距离从小到大排序后，为了更好地用网格点的预报雨量来预测观测站点的雨量，我们取前4个到观测站点距离最小的等距网格点。

根据欧氏距离的倒数加权的方法，先算出前4个网格点到观测点的距离，再分别对它们求倒，则4个网格点分别到观测站点的权重为它们之间距离的倒数。

权的计算公式为：)4,3,2,1,91,...2,1(,11
41===∑
=n i d d q j ij
in in 其中为了预测各观测站点在某月某日某个时段的雨量值，我们采用距离的倒数加权的方法，取出4个等距网格点分别在某月某日某个时段的雨量值，然后分别乘以它们各自对观测站点的权重，再求和就为预测降雨量。

预测降雨量的计算公式：∑=×=4
1n in
in ij f q y 每个时段中，每个观测站点对应有4个网格点预测雨量值，可计算出1个观测站点预测雨量值，91个观测站点，164个时段就可计算出一个观测站点的时段预测降个第i 个第j 雨量矩阵)
164,...,2,1,91,...,2,1()(===×j i y M j i ij 这里我们用Matlab 编程求得两个方法对应的预测降雨量矩阵（见附录三）。

4．
将两个方法对应的预测降雨量矩阵分别与实测降雨量矩阵进行比较，分析出哪一个的准确性高。

这里我们用计算出预报偏差率的算术平方根作为一个准确性指数，来辨别准确性的高低。

取预测降雨量矩阵和实测降雨量矩阵)164,...,2,1,91,...,2,1()(===×j i y M j i ij 。

则有：
)164,...,2,1,91,...,2,1()(===×j i s I j i ij 预报偏差率计算公式：)
164,...,2,1,91,...,,2,1(==−=j i y s y ij ij
ij ij σ取中的元素，计算预报偏差率的算术平方根()j i ij ×σ()∑∑===91116412
i j ij
S σ当S 越小，准确性越高。

这里我们用Matlab 编程求得两个方法（见附录四）求得：预报偏差率的算术平方根（即准确性指数）：
方法一的为：8755
.1021=S 方法二的为：7522
.7262=S 所以可以反映出第一种方法比第二种方法准确性高。

二、问题问题（（2）及其求解
由题意可得：气象部门将6小时降雨量分为6等：0.1—2.5毫米为小雨，2.6—6毫米为中雨，6.1—12毫米为大雨，12.1—25毫米为暴雨，25.1—60毫米为大暴雨，大于60.1毫米为特大暴雨。

为了比较91个观测点的预报数据与实测数据之间量级上的差别，分别将降雨量的预报值和实测值按大小划分成7个级别后，分别记为：[0，0.1]——0，[0.1，2.5]——1，[2.6，6]——2，[6.1，12]——3，[12.1，25]——4，[25.1，60]——5，[60.1，——6。

然后分别统计预报数据与实测数据处在同一级别、相差1级、相差2级、相]∞差3级、相差4级、相差5级、相差6级的频数，并计算出对应频率。

把两个方法对应的预测降雨量矩阵（题一中已计算出）与实测降雨量矩阵转化为对应的等级矩阵，将两个方法对应的预测等级矩阵分别与实测等级矩阵进行比较，分别统计预报数据与实测数据处在同一级别、相差1级、相差2级、相差3级、相差4级、相差5级、相差6级的频数。

我们通过Matlab 编程求解（见附录四），得出结果：
方法一的为：方法二的为：0,
0,
2,
17,113,
3758,1103416151413121110=======J J J J J J J 0
,
11,19,56,
145,
4102,
1059126252423222120=======J J J J J J J 并计算出频率：
表一
由表一中和图二可以看出方法一的报错率低，相对报错等级差小，公众应该更满意使用方法一进行降雨量的预报。

方法一频数方法二频数方法一频率
方法二频率11034105910.73934602
0.70966229375841020.25180917
0.274859291131450.0075717
0.0097158917560.0011391
0.003752352190.00013401
0.001273120110
0.0007370700
00
四、模型的误差分析
1、在可接受度数差范围内搜索与各个观测站点距离最近的网格点取的较少，如1个点
时，则不能很好地测出预测雨量，具有一定的偶然性。

2、当与各个观测站点距离最近的网格点取的比较多，如15个点时，则每个点都有一定
的权重，具有大面积的网格点来测雨量，这对与观测站点较近的网格点来说是很不公平的。

所以，我们为了得到较为合理的结果，在本模型中取与各个观测站点距离最近的4个网格点。

五、模型的评价
1、本模型利用欧氏距离计算方法，对等距网格点和观测站点之间的距离进行了求解，
简单明了，通俗易懂。

2、本模型综合运用了多种方法对问题进行求解。

其中欧氏距离倒数加权的方法对预测
观测站点的雨量起到了极为重要的作用。

3、利用算术平方根方法来反映预报数据的偏差率，以及在数据的准确性方面有着非常
高的精度。

六、模型的推广
此模型是针对均匀网格的雨量预报方法的评价，而大多数数值天气预报模式都是采用这种单一均匀的网格来进行计算。

从60年代开始逐步研究使用不同格距的网格作预报，也就是变网格预报模式。

变网格(又称其为非均匀网格)能减轻网格之间差异带来的问题，而且以平缓渐变的网格为优。

设计有限区域变网格模式的计算方案，可采用坐标变换使不均匀网格变成均匀网格再进行计算；或利用人工构造变网格，用有限元法求解等。

对于有限区域，如重点考虑某一局地天气，在模式中使用嵌套网格技术是最常见的方法，它是将更细的网格嵌套在所关心区域的粗网格上，以此来提高局地空间分辨率。

但是嵌套网格技术存在一个缺陷，那就是在内边界附近，由于网格距的突变而导致空间截断误差的突变，造成波动由粗网格进入细网格(或相反)时产生较严重的折射、反射和寄生波现象。

虽然现在对变网格研究较多，但投入业务使用的只有加拿大气象中心和法国气象局。

本文是将原来均匀网格的MM4模式改造成非均匀变网格模式，然后建立一套投入业务使用的客观分析系统、变网格数值预报系统、预报资料图形显示系统和网络自动化传输资料系统等。

下图三给出了变网格和均匀网格500hPa高度场预报的均方差随预报时间的演变图。

从图中可看出，两种网格模式预报的均方差都是随着预报时间的增加而增加的，变网格预报的均方差小于均匀网格预报的均方差，这说明变网格数值模式在天气形势预报上也是优于均匀网格的。

图三1999年8月8日08时均匀网格和变网格
所以本系统经过实际数值预报检验，得到在基本没有增加计算量的情况下，变网格数值预报模式对所关心区域的环流形势和局部天气系统的预报明显优于传统的均匀网格预报，证明了变网格数值预报是提高数值预报准确率的一种非常有效的手段，开展变网格数值预报的研究和应用是非常有意义的。

参考文献
[1]王沫然编，《Matlab与科学计算，第二版》，电子工业出版社，北京，2003。

[2]姜启源等编，《数学模型》，高等教育出版社，2003。

[3]李火林等编，《数学模型及方法》，江西高校出版社，1997。

[4]毛恒青等编，《简化典型相关预测模型面雨量中期预报试验》，《气象》。

[5]毛恒青，李小泉.我国夏季降水与前期太平洋海温场关系的典型相关分析.南京气象学院学报，1998，21(1)：130～137.
[6]曾庆存，李荣凤.不等距差分格式的计算紊乱问题.大气科学，1982，6(4)：345～354.
附录
***********注意：运行时需按步骤进行，中途不可清除Matlab中的变量.********* 1。

运行程序二。

（搜索实测站点周围在可接受度数的网格点，然后根据距离长短排序并保存）
2。

选择对应编号的文件，导入Y=[降雨量]53*47的数据，
修改编号:
M1(k,编号)=sum;%M1在第一个程序中定义了，算第二个方案时M1换为M2
运行程序三。

（根据模型利用距离的倒数加权计算每天每个时段预测的结果）
3。

修改:
M1算第二个方案时M1换为M2
运行程序四
（预测矩阵与实测矩阵的数据分析）
1、附录一、等距网格点和观测站点的散点分布图程序
x=[纬度];%x是由文件夹FORECAST中lat文件的数据生成的一个纬度值矩阵,x 大小为53*47
y=[经度]；%y是由文件夹FORECAST中lon文件的数据生成的一个经度值矩阵,y 大小为53*47
lat=[32.983333.300033.666733.800033.483333.033333.2333 33.766733.383333.200032.100032.300032.000032.683332.8000 32.933332.416732.333332.200032.866732.183332.533332.3833 32.333332.066731.800031.950031.333331.566731.700031.0833 31.983331.750031.766731.950031.433331.366731.266731.8833 31.650031.583331.416731.066731.150031.900031.100031.4000 31.366731.616731.200031.050031.233331.466730.300030.8500 30.683330.933330.300030.616730.083330.883331.133330.9833 30.966730.233330.050030.850030.850030.783330.000030.5333 30.516730.633330.233330.200030.883331.133331.000030.9333 30.616730.266730.066730.733330.033330.250029.866729.7167 29.783329.816729.700029.9667];
lon=[118.5167118.8500119.2667119.8000119.8167119.0333119.3000 120.2500120.1500120.4833118.2667118.3000118.8000119.0167119.4500 119.8333119.4167119.9333120.0000120.3167119.4667120.4500120.5667 121.1833121.6000121.6667118.8500118.3833118.5000118.5167118.1833 119.5833119.5500119.9333119.1667119.4833119.8167120.6333120.2667 120.7333120.3167120.9500120.4333120.6333121.2000121.3667121.4833 121.2500121.4500121.4333121.7833121.5333121.1000118.1333118.3167 118.4000118.7500118.5333118.9833118.5833119.4167119.1833119.8833 119.6833119.7000119.9500120.0833120.9000120.7333120.6333120.0667 120.6833120.5333120.1667120.3167121.1667121.1167121.2500121.4833 121.0833121.2167121.1500122.4500122.1000122.1833118.4333118.2833 118.1833119.6833120.2500121.7500];
plot(lat,lon,'dk',x,y,'.g');
xlabel('纬度','FontWeight','bold');
ylabel('经度','FontWeight','bold');
legend('实测站点','网格结点')
axis([2736117125])
3、附录二搜索观察站点周围最近的结点以及根据距离长短排序
a=[经度]；%a是由文件夹FORECAST中lon文件的数据生成的一个经度值矩阵,a大小为53*47
b=[纬度]；%b是由文件夹FORECAST中lat文件的数据生成的一个纬度值矩阵,b大小为53*47
c=[经度，纬度];%c为91个实测站点对应的经度和纬度,c大小为91*2
g=zeros(91,20);
e=0.25;%e为可接受度数差
for k=1:91
m=1;
for i=1:53
for j=1:47
if abs(b(i,j)-c(k,1))<=e&abs(a(i,j)-c(k,2))<=e%在可接受度数差范围内，搜索几个网格点
d1(k,m)=i;
d2(k,m)=j;
s=(c(k,1)-b(i,j))^2+(c(k,2)-a(i,j))^2;
g(k,m)=sqrt(s);
m=m+1;
end
end
end
N(k)=m-1;
p=[1:20;g(k,:)];
for c1=1:N(k)-1%计算所搜索的网格点与实测站点的距离，并排序后保存
for c2=1:N(k)-c1
if p(2,c2)>p(2,c2+1)
c3=p(2,c2);
p(2,c2)=p(2,c2+1);
p(2,c2+1)=c3;
c4=p(1,c2);
p(1,c2)=p(1,c2+1);
p(1,c2+1)=c4;
end
end
end
for c5=1:N(k)
J1(k,c5)=d1(k,p(1,c5));
J2(k,c5)=d2(k,p(1,c5));
G(k,c5)=g(k,p(1,c5));
end
end
d1=d1;d2=d2;g=g;G=G;J1=J1;J2=J2;N=N;
4、附录三根据模型利用距离的倒数加权计算每天每个时段预测的结果
Y=[降雨量];%Y大小53*47,为某个时段的对应降雨量矩阵,其数据取于FORECAST 中该时段降雨量文件，
如:矩阵调用f6194_dis2中的53*47个数据,按方案中测量时间先后编号为8.
（矩阵共两种方案164个=2*164个）
for k=1:91
R=0;
for t1=1:4
if G(k,t1)==0
break;
else R=1/G(k,t1)+R;%R为权的分母
end
end
sum=0;
for t2=1:4
y=Y(J1(k,t2),J2(k,t2));
if G(k,t1)==0%若当某站点正好位于网格点上，则该网格点预测数据为该
站点预测数据
sum=y;
break;
else w=1/G(k,t2);
sum=(w/R)*y+sum;%对所取数据加距离权处理
end
end
M1(k,填入编号)=sum;%M1=[预测降雨量]91*164;算第二个方案时M1换为M2
end
5、附录四预测方法
m0=0;m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;
I=[实际降雨量];%I=[实际降雨量]91*164，实际降雨量矩阵
for s1=1:2
if s1==1
D=M1;%M1=[预测降雨量]91*164;当取第二个方案时，改M1为M2。

else D=I;
end%D=[降雨量]91*164，暂时保存降雨量的矩阵DJ=zeros(91,164);%DJ=[等级]91*164，暂时保存对应等级的矩阵
for i=1:91
for j=1:164
D(i,j)=round(D(i,j)*10);
if D(i,j)<1
DJ(i,j)=0;
elseif D(i,j)>=1&D(i,j)<=25
DJ(i,j)=1;
elseif D(i,j)>=26&D(i,j)<=60
DJ(i,j)=2;
elseif D(i,j)>=61&D(i,j)<=120
DJ(i,j)=3;
elseif D(i,j)>=121&D(i,j)<=250
DJ(i,j)=4;
elseif D(i,j)>=251&D(i,j)<=600
DJ(i,j)=5;
else
DJ(i,j)=6;
end
end
end
if s1==1
MM1=DJ;
else
II=DJ;
end
end
E1=abs(M2-I);%算出两矩阵对应元素的偏差值，存于E1中
sum=0;
for r1=1:91
for r2=1:164
r3=M2(r1,r2);
if r3==0
r3=0.01;
end
sum=(E1(r1,r2)/r3)^2+sum;
x=abs(MM1(r1,r2)-II(r1,r2));
if x==0
m0=m0+1;
elseif x==1
m1=m1+1;
elseif x==2
m2=m2+1;
elseif x==3
m3=m3+1;
elseif x==4
m4=m4+1;
elseif x==5
m5=m5+1;
else
m6=m6+1;
end
end
end
sum=sqrt(sum)
m0=m0
m1=m1
m2=m2
m3=m3
m4=m4
m5=m5
m6=m6。