雨量预报数学模型及求解代码.docx
《雨量预报数学模型及求解代码.docx》由会员分享，可在线阅读，更多相关《雨量预报数学模型及求解代码.docx（17页珍藏版）》请在冰豆网上搜索。
雨量预报数学模型及求解代码
赛区评阅编号（由赛区组委会填写）：
2015高教社杯全国大学生数学建模竞赛
承诺书
我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。
如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。
我们参赛选择的题号（从A/B/C/D中选择一项填写）：
我们的报名参赛队号（12位数字全国统一编号）：
参赛学校（完整的学校全称，不含院系名）：
参赛队员（打印并签名）：
1.张伟尧
2.熊霜
3.韩倩
指导教师或指导教师组负责人（打印并签名）：
董银丽
日期：
2015年9月6日
赛区评阅编号（由赛区组委会填写）：
2015高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅记录（可供赛区评阅时使用）：
评
阅
人
备
注
送全国评阅统一编号（由赛区组委会填写）：
全国评阅随机编号（由全国组委会填写）：
雨量预报方法的评价模型
摘要
本文建立了一个关于雨量预报方法的评价模型。
首先，通过给定的大量数据（预测数据和实测数据）利用Matlab图形处理功能的基本绘图命令plot画出2491个网格点和91个观测点的位置。
在可接受范围内，计算各观测站点和等距网格点之间的距离，并按升序排序。
取其前5个到观测站点距离最小的等距网格点，再根据欧拉公式距离倒数加权的方法对它们赋权重，取其前5个网格点的雨量，分别乘以它们各自对应的权并求和，就是相对应观测站点的雨量。
分别对两种方法预报的41天每天4个时段各网格点的雨量处理，最后对两种方法得到站点的预测雨量进行分析。
问题一：
利用上述得到的站点预测雨量，计算某个观测点在某时段的预报偏差率，并对其预报偏差率求和再算术平方根得出标准率，作为评价准确性高低的指标，从而得到第一种雨量预报方法的准确率为
，第二种雨量预报方法准确率为
。
因此得出第一种方法比第二种方法的准确性高。
问题二：
分别将两种方法的站点预测雨量与实测雨量按照题目给定的等级划分，得到等级矩阵，然后分别用两种方法的预测等级矩阵和实测等级矩阵作差，统计同一级别、相差1级、相差2级……相差6级数据的频数和频率。
等级差越小，相对报错率越低，公众越满意。
因此得出第一种方法比第二种方法的公众满意度高。
关键词：
欧拉公式权标准率等级差
一、问题重述
雨量预报对农业生产和城市工作和生活有重要作用，但如何准确、及时的对雨量做出预报是一个十分困难的问题。
本题由气象部门提出，希望建立一种科学评价预报方法好坏的数学模型与方法。
我国某地气象台和气象研究所正在研究6小时雨量预报方法，即每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度，北纬32度附近的53
47的等距网格点上。
同时设立91个分布不均匀的观测站点实测这些时段的实际雨量。
气象部门提供了41天用两种不同方法的预报数据和相应的实测数据，预测数据在文件夹FORECAST中，实测数据在文件夹MEASURING中，其中的文件都可以用Windows系统的“写字板”程序打开阅读。
经（lon.dat）纬（lat.dat）度也分别包括在文件夹FORECAST中，其余文件名为_dis1和_dis2表示该日期用的第一或第二种方法。
题目中有如下假设：
1、雨量用毫米作单位，小于0.1毫米的视为无雨。
2、气象部门将6小时降雨量分为6等：
0.1—2.5毫米为小雨，2.6—6毫米为中雨，6.1—12毫米为大雨，12.1—25毫米为暴雨，25.1—60毫米为大暴雨，大于60.1毫米为特大暴雨。
题目的问题如下：
问题一：
请建立数学模型来评价两种6小时雨量预报方法的准确性；
问题二：
若按此分级向公众预报，如何在评价方法中考虑公众的感受
二、问题分析
我们从题目中了解分析得到：
气象台每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度，北纬32度附近的53
47的等距网格点上。
同时设立91个分布不均匀的观测站点实测这些时段的实际雨量。
由于网格点比较多，且每个网格点的位置是以经度和纬度表示在一定的区域，所以我们将纬度看作x轴，经度看作y轴，采用Matlab图形处理功能的基本绘图命令plot画出散点图（图1），（程序见附录一）。
图中绿色表示网格结点的分布，红色表示实测站点的分布。
从图中可以分析看出，气象部门提供了2491个网格点上41天4个时间段的大量预报数据（雨量），并且同样给出了91个观测站点的实测数据（雨量）。
所以我们想通过网格上的预报数据来预测实测站点的数据。
然而，观测站点集中在所有网格的中央部分，而四周是大量的距离比较远的网格点。
因此，通过搜索出2491个网格点中对站点影响比较大的几个网格点，再用搜索出来的几个网格点的预测数据加权求出一个预测数据（雨量），进而和该站点实测数据进行比较，来评价两种6小时雨量预报方法的准确性。
在向公众预报时，采用一种合理、准确的预测方法，增加雨量分等级预报的同级率，能对公众的出行起到良好的指导作用，使人们对雨量预报有更深的理解，更多的关注。
图1、网格结点与实测站点的散点分布图
三、模型假设
1、观测站点的设置是不均匀的；
2、题中的网格式等间距的正方形网格（所谓“正方形网格”是指每个格子都是正方形的网格；网络点是指纵线和横线的交叉点）；
3、一个x轴、y轴分别是纬度和经度的坐标，通过把点的纬度和经度分别看作横坐标和纵坐标，用欧拉公式计算
来作为两点之间的距离。
4、点到观测站点的距离越短，则对观测站点的雨量影响越大；
5、单个网格点到观测站点倒数与所取的5个网格点到观测点倒数之和的比为它的权值；
6、雨量用毫米做单位，小于0.1毫米视为无雨；
四、符号说明
：
与第
个观测站点的距离最小的前5个网格点的对应权矩阵；
：
与第
个观测站点的距离最小的前5个网格点的距离矩阵；
：
欧拉公式计算第
个观测站点到第
个实测站点的距离计算；
：
第
个观测站点的坐标点的纬度和经度；
：
可接受度数差；
：
第
个实测站点的距离最小的前5个网格点到第
个观测站点的权重；
：
第
个观测站点在第
个时段的预测雨量之和；
：
第
个观测站点在第
个时段的预测雨量矩阵；
：
第
个观测站点在第
个时段的实际雨量矩阵；
：
第
个观测站点在第
个时段的预报偏差率；
：
预报偏差率的算术平方根；
：
第
种方法的预报数据与实测数据处在同一级别，相差同一级别,相差1级……相差6级的频数
五、模型的建立与求解
5、1问题一模型建立及其求解
由于问题涉及的数据量较大，为了模型计算的简便，需要对数据进行处理。
（一）数据处理：
为了评价两种6小时雨量预报方法的准确性，我们采用网格点上的预报数据来预测观测站点的数据，再来和实际测得的数据相比，判断其准确性，首先从预测数据选取可接受范围内的网格点，其次选取实测点附近5个距离最小的预测点，计算5个网格点分别到观测站点的权重，并计算观测站点的预测雨量。
（1）可接受范围内网格点的选取
以坐标
为基准点，给定一个可接受度数差
（在求解中取
，对任意的
，搜索其任一个观测站点在纬度和经度都上下增加
的正方形内所有等距网格点。
若网格点
纬度和经度在同时满足
和
时，即认为该网格点时可接受范围内的网格点。
（2）网格点到站点距离的计算
找到可接受范围内的网格点后，我们计算网格点和这些观测站点间的距离。
再得到个观测站点和等距网格点之间的距离后，将各观测站点按距离从小到大排序后保存，我们用Matlab编程求得结果（程序见附录二）。
（3）权的计算
各观测站点和等距网格点之间的距离从小到大排序后，为了更好地用网格点的预报雨量来预测观测站点的雨量，我们取前5个到观测站点距离最小的等距网格点。
根据欧拉公式距离的倒数加权的方法，先算出前5个网格点到网格点的距离，再分别对它们求倒，则5个网格点分别到观测站点的权重为它们之间距离的倒数。
权的计算公式为：
（4）观测站点预测雨量的计算
为了预测各观测站点在某月某日某个时段的雨量值，我们采用距离的倒数加权的方法，取出5个等距网格点分别在某月某日某个时段的雨量值，然后分别乘以它们各自对观测站点的权重，再求和就为预测降雨量。
预测降雨量的计算公式：
；
每个时段中，每个观测点对应有5个网格点预测雨量值，可计算出1个观测站点预测雨量值，91个观测站点，164个时段就可计算出一个第
个观测站点的第
个时段预测降雨量矩阵
。
这里我们用Matlab编程求得两个方法对应的预测降雨量矩阵（程序见附录三）。
（二）模型建立及求解：
将两个方法对应的预测降雨量矩阵分别与实测降雨量矩阵进行比较，分析出哪一个的准确性高。
这里我们用计算出预报偏差率的算术平方根作为一个准确性指数，来辨别准确性的高低。
预测降雨量的矩阵
和实际降雨量矩阵
。
则有：
预报偏差率计算公式：
；
取
中的元素，计算预报偏差率的算术平方根
，
当
越小，准确性越高。
这里我们用Matlab编程求得两个方法（程序见附录四）求得：
预报偏差率的算术平方根（即准确性指数）：
方法一的为：
；
方法二的为：
；
所以可以反映出第一种方法比第二种方法准确性高。
5．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.5,6]——2；（6,12]——3；（12,25]——4；（25,60]——5；（60,
）——6.然后分别统计预报数据与实测数据处在同一级别、相差1级、相差2级、相差3级、相差4级、相差5级、相差6级的频数，并计算出对应频率。
把两个方法对应的预测降雨量矩阵（题一中已计算出）与实测雨量矩阵转化为对应的等级矩阵，将两个方法对应的预测等级矩阵进行比较，分别统计预报数据与实测数据处在同一级别、相差1级、相差2级、相差3级、相差4级、相差5级、相差6级的频数。
我们通过Matlab编程求解（程序见附录五），得出结果：
方法一的为：
方法二的为：
并计算出频率（如表一）：
方法一频数
方法二频数
方法一频率
方法二频率
11638
11626
77.9818%
77.9014%
3127
3139
20.9528%
21.0332%
136
136
0.9113%
0.9113%
20
20
0.1340%
0.1340%
3
3
0.0201%
0.0201%
0
0
0.0000%
0.0000%
0
0
0.0000%
0.0000%
表一：
频率
图2、等级差频率分布饼状图
由表一中和图2可以看出方法一的报错率低，相对报错等级差小，公众应该更满意使用方法一进行雨量预报。
六、模型的误差分析
1、在可接受范围内搜索与各个观测站点距离最近的网格点取的较少，如1个点时，则不能很好地测出预测雨量，具有一定的偶然性。
2、当与各个观测站点距离最近的网格点取的比较多，如15个点时，则每个点都有一定的权重，具有大面积的网格点来测雨量，这对于观测站点较近的网格点来说是很不公平的。
3、在可接受范围点中有三个网格点的坐标与三个站点的坐标一致，对于这三个站点的预测降雨量，可直接调用该坐标处网格点位置的降雨量。
所以，我们为了得到较为合理的结果，在本模型中取与各个观测站点距离最近的5个网格点。
七、模型的评价
1、本模型利用欧拉公式距离计算方法，对等距网格点和观测站点之间的距离进行求解，简明易了，通俗易懂。
2、本模型综合运用了多种方法对问题进行求解。
其中欧拉公式距离倒数加权的方法对预测观测站点的雨量起到了极为重要的作用。
3、利用算术平方根方法来反映预报数据的偏差率，以及在数据的准确性方面有着非常高的精度。
参考文献
[1]刘来福，曾文艺，《数学模型与数学建模》，北京，北京师范大学出版社，1997年；
[2]H.P.Williams,《数学规划模型建立与计算机应用》，北京，国防工业出版社，1991年；
[3]卓金武，《MATLAB在数学建模中的应用》，北京，北京航空航天大学出版社，2011年；
[4]周品、赵新芬，《MATLAB数理统计分析》，北京，国防工业出版社，2009年；
[5]，气象数据matlab处理_图文_XX文库；；
附录
附录一
等距网格点和观测站点的散点分布图程序（在文件FORECAT1中把lat,lon文件导入到Matlab中，在文件MEASURING把A020618.SIX导入到Matlab中）
x=lat;
y=lon;
xi=A020618（:
2）;
yi=A020618（:
3）;
plot（x,y,'.g',xi,yi,'.r'）%2491个预测点和91个观测点分布图
axis（[27.335.3116.8125]）;
附录二
可接受范围内的网格点到站点距离的计算
xi=A020618（:
2）;
yi=A020618（:
3）;
plot（x,y,'.g',xi,yi,'.r'）
axis（[27.335.3116.8125]）;
shice_x1=max（xi）+0.2;
shice_x2=min（xi）-0.2;
shice_y1=max（yi）+0.2;
shice_y2=min（yi）-0.2;
b1=find（x>=shice_x2&x<=shice_x1）;%满足条件的点的x坐标位置
xb2=x（b1）;
yb2=y（b1）;
plot（xb2,yb2,'.g',xi,yi,'.r'）;
b3=find（yb2>=shice_y2&yb2<=shice_y1）;%筛选后满足条件的点的y坐标位置
xb3=xb2（b3）;%有效的点
yb3=yb2（b3）;
plot（xb3,yb3,'.g',xi,yi,'.r'）;
axis（[29.334.3117.8122.8]）;
xb4=xi;
yb4=yi;
l=length（xb3）;
d=zeros（91,l）;
fori=1:
91
forj=1:
l
d（i,j）=sqrt（（xb4（i）-xb3（j））^2+（yb4（i）-yb3（j））^2）;%计算距离
end
end
附录三
观测站点预测降雨量的计算
[B,C]=sort（d,2）;%对d从小到大排序，并记录在原矩阵中的坐标
D=C（:
[1,2,3,4,5]）;
E=B（:
[1,2,3,4,5]）;
F=1./E;
F（:
6）=sum（F,2）;
[r,l]=size（F）;
fori=1:
r
K（i,1:
l）=F（i,:
）/max（max（F（i,:
）））;%每个最小值点权的计算
end
H=K（:
[1,2,3,4,5]）;
H（find（isnan（H）==1））=1%将H中NaN替换为1
HH=reshape（H',455,1）';%前五个最小值的权并排序
a=b1（b3（D））;%返回91*5个点在x中的矩阵位置
%%%%%%%%%%%%%%%%%%%%%预测数据调用和降雨量的筛选%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%将CurrentDirectory路径改为数据所在位置%%%%%%%%%%%%%
A1=[];A2=[];
formonth1=6:
7
ifmonth1==6;
fordata1=18:
28
fortime1=1:
4
str1=['f'int2str（month1）int2str（data1）int2str（time1）'_dis1'];
str2=['f'int2str（month1）int2str（data1）int2str（time1）'_dis2'];
tmp1=load（str1）;tmp2=load（str2）;
A1=[A1;tmp1（:
）'];A2=[A2;tmp2（:
）'];
end
end
else
fordata1=1:
30
ifdata1<10
strtmp=['0'int2str（data1）];
else
strtmp=int2str（data1）;
end
fortime1=1:
4
str1=['f'int2str（month1）strtmpint2str（time1）'_dis1'];
str2=['f'int2str（month1）strtmpint2str（time1）'_dis2'];
tmp1=load（str1）;tmp2=load（str2）;
A1=[A1;tmp1（:
）'];A2=[A2;tmp2（:
）'];
end
end
end
end
%A1为第一种方法统计的所有数据，A2为第二种方法统计的所有数据
w=reshape（a',455,1）';%将91*5个点转换成455*1
A11=A1（:
w）;%第一类164*455个降雨量
A22=A2（:
w）;%第二类164*455个降雨量
gg=repmat（HH,164,1）%为使HH与A11，A22大小相同，扩展平铺HH，gg大小为164*455
A1_1=gg.*A11%第一类方法权与降雨量乘积
A2_2=gg.*A22%第二类方法权与降雨量乘积大小为164*455
%%%%%%%%分别对A1_1，A2_2每隔5列求和并返回sum1,sum2%%%%%%%%%
sum1=[];
fori=1:
5:
454
sum1=[sum1A1_1（:
i）+A1_1（:
i+1）+A1_1（:
i+2）+A1_1（:
i+3）+A1_1（:
i+4）];
end
sum2=[];
fori=1:
5:
454
sum2=[sum2A2_2（:
i）+A2_2（:
i+1）+A2_2（:
i+2）+A2_2（:
i+3）+A2_2（:
i+4）];
end
附录四
模型的求解
%%%%%%%%%%%%%%%91个站实测数据的导入%%%%%%%%%%%%%%%%%%%%%%
%将CurrentDirectory路径改为数据所在位置
B1=[];
formonth1=6:
7
ifmonth1==6;
fordata1=18:
28
%str3=['020'int2str（month1）int2str（data1）];
str3=['020'int2str（month1）int2str（data1）'.SIX'];
tmp3=load（str3）;
B1=[B1;tmp3（:
4:
7）'];
end
else
fordata1=1:
30
ifdata1<10
strtmp=['0'int2str（data1）];
else
strtmp=int2str（data1）;
end
str3=['020'int2str（month1）strtmp'.SIX'];
tmp3=load（str3）;
B1=[B1;tmp3（:
4:
7）'];
end
end
end
B_1=abs（sum1-B1）./sum1;%第一种方法的偏差率
B_2=abs（sum2-B1）./sum2;%第二种方法的偏差率
B_1（find（isnan（B_1）==1））=0;%将B_1中NaN替换成为0
B_2（find（isnan（B_2）==1））=0;%将B_2中NaN替换成为0
B_11=B_1.^2
S1=sqrt（sum（B_11（:
）））%第一种方法的算术平方根，即标准率
B_22=B_2.^2
S2=sqrt（sum（B_22（:
）））%第二种方法的算术平方根，即标准率
附录五
问题二的求解
%%%%%将sum1，sum2，B1降雨量数据分等级，形成164*91大小的等级矩阵
sizesum1=size（sum1）;
z1=zeros（sizesum1）;
z2=z1;
z3=z1;
z4=z1;
z5=z1;
z6=z1;
z7=z1;
z1（sum1<0.1）=0;
z2（sum1>=0.1&sum1<=2.5）=1;
z3（sum1>2.5&sum1<=6）=2;
z4（sum1>6&sum1<=12）=3;
z5（sum1>12&sum1<=25）=4;
z6（sum1>25&sum1<=60）=5;
z7（sum1>60）=6;
sum11=z1+z2+z3+z4+z5+z6+z7;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sizesum2=size（sum2）;
r1=zeros（sizesum2）;
r2=r1;
r3=r1;
r4=r1;
r5=r1;
r6=r1;
r7=r1;
r1（sum2<0.1）=0;
r2（sum2>=0.1&sum2<=2.5）=1;
r3（sum2>2.5&sum2<=6）=2;
r4（sum2>6&sum2<=12）=3;
r5（sum2>12&sum2<=25）=4;
r6（sum2>25&sum2<=60）=5;
r7（sum2>60）=6;
sum22=r1+r2+r3+r4+r5+r6+r7;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sizeB1=size（B1）;
t1=zeros（sizeB1）;
t2=t1;
t3=t1;
t4=t1;
t5=t1;
t6=t1;
t7=t1;
t1（B1<0.1）=0;
t2（B1>=0.1&B1<=2.5）=1;
t3（B1>2.5&B1<=6）=2;
t4（B1>6&B1<=12）=3;
t5（B1>12&B1<=25）=4;
t6（B1>25&B1<=60）=5;
升级会员 