雨量预报数学模型及求解代码

1、仅供个人学习参考赛区评阅编号（由赛区组委会填写）：2015高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了全国大学生数学建模竞赛章程和全国大学生数学建模竞赛参赛规则（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛章程和参赛规则，以

2、保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。我们参赛选择的题号（从A/B/C/D中选择一项填写）：我们的报名参赛队号（12位数字全国统一编号）：参赛学校（完整的学校全称，不含院系名）：参赛队员（打印并签名）：1张伟尧2熊霜3韩倩指导教师或指导教师组负责人（打印并签名）：董银丽日期：2015年9月6日赛区评阅编号（由赛区组委会填写）：-2015高教社杯全国大学生数学建模竞赛编号专用页赛区评阅记录（可供赛区评阅时使用）：全国

3、评阅随机编号（由全国组委会填写）：雨量预报方法的评价模型摘要本文建立了一个关于雨量预报方法的评价模型。首先，通过给定的大量数据（预测数据和实测数据）利用Matlab图形处理功能的基本绘图命令plot画出2491个网格点和91个观测点的位置。在可接受范围内，计算各观测站点和等距网格点之间的距离，并按升序排序。取其前5个到观测站点距离最小的等距网格点，再根据欧拉公式距离倒数加权的方法对它们赋权重，取其前5个网格点的雨量，分别乘以它们各自对应的权并求和，就是相对应观测站点的雨量。分别对两种方法预报的41天每天4个时段各网格点的雨量处理，最后对两种方法得到站点的预测雨量进行分析。问题一：利用上述得到的

4、站点预测雨量，计算某个观测点在某时段的预报偏差率，并对其预报偏差率求和再算术平方根得出标准率，作为评价准确性高低的指标，从而得到第一种雨量预报方法的准确率为S二104.5506，第二种雨量预报方法准确率为1S=104.6044。因此得出第一种方法比第二种方法的准确性高。2问题二：分别将两种方法的站点预测雨量与实测雨量按照题目给定的等级划分，得到等级矩阵，然后分别用两种方法的预测等级矩阵和实测等级矩阵作差，统计同一级别、相差1级、相差2级相差6级数据的频数和频率。等级差越小，相对报错率越低，公众越满意。因此得出第一种方法比第二种方法的公众满意度高。关键词：欧拉公式权标准率等级差一、问题重述雨量预

5、报对农业生产和城市工作和生活有重要作用，但如何准确、及时的对雨量做出预报是一个十分困难的问题。本题由气象部门提出，希望建立一种科学评价预报方法好坏的数学模型与方法。我国某地气象台和气象研究所正在研究6小时雨量预报方法，即每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度，北纬32度附近的53x47的等距网格点上。同时设立91个分布不均匀的观测站点实测这些时段的实际雨量。气象部门提供了41天用两种不同方法的预报数据和相应的实测数据，预测数据在文件夹FORECAST中，实测数据在文件夹MEASURIN

6、G中，其中的文件都可以用Windows系统的“写字板”程序打开阅读。经（lon.dat）纬（lat.dat）度也分别包括在文件夹FORECAST中，其余文件名为f日期i_dis1和f日期i_dis2表示该日期用的第一或第二种方法。题目中有如下假设：1、雨量用毫米作单位，小于0.1毫米的视为无雨。2、气象部门将6小时降雨量分为6等：0.12.5毫米为小雨，2.66毫米为中雨，6.112毫米为大雨，12.125毫米为暴雨，25.160毫米为大暴雨，大于60.1毫米为特大暴雨。题目的问题如下：问题一：请建立数学模型来评价两种6小时雨量预报方法的准确性；问题二：若按此分级向公众预报，如何在评价方法中考

7、虑公众的感受？二、问题分析我们从题目中了解分析得到：气象台每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度，北纬32度附近的5347的等距网格点上。同时设立91个分布不均匀的观测站点实测这些时段的实际雨量。由于网格点比较多，且每个网格点的位置是以经度和纬度表示在一定的区域，所以我们将纬度看作X轴，经度看作y轴，采用Matlab图形处理功能的基本绘图命令plot画出散点图（图1），（程序见附录一）。图中绿色表示网格结点的分布，红色表示实测站点的分布。从图中可以分析看出，气象部门提供了2491个网格

8、点上41天4个时间段的大量预报数据（雨量），并且同样给出了91个观测站点的实测数据（雨量）。所以我们想通过网格上的预报数据来预测实测站点的数据。然而，观测站点集中在所有网格的中央部分，而四周是大量的距离比较远的网格点。因此，通过搜索出2491个网格点中对站点影响比较大的几个网格点，再用搜索出来的几个网格点的预测数据加权求出一个预测数据（雨量），进而和该站点实测数据进行比较，来评价两种6小时雨量预报方法的准确性。在向公众预报时，采用一种合理、准确的预测方法，增加雨量分等级预报的同级率，能对公众的出行起到良好的指导作用，使人们对雨量预报有更深的理解，更多的关注。125124123122度121经1

9、2011911811727282930313233343536纬度图1、网格结点与实测站点的散点分布图三、模型假设耀莓左；lymi拆1、观测站点的设置是不均匀的；2、题中的网格式等间距的正方形网格(所谓“正方形网格”是指每个格子都是正方形的网格；网络点是指纵线和横线的交叉点)；3、一个x轴、y轴分别是纬度和经度的坐标，通过把点的纬度和经度分别看作横坐标和纵坐标，用欧拉公式计算d=.;(a-a)2+(b-b)2来作为两点之间的距离。ijij4、点到观测站点的距离越短，则对观测站点的雨量影响越大；5、单个网格点到观测站点倒数与所取的5个网格点到观测点倒数之和的比为它的权值；6、雨量用毫米做单位，小

10、于0.1毫米视为无雨；四、符号说明Q二(q)(i=1,2.,91,n二1,2,3,4,5):与第i个观测站点的距离最小的前5个网格点ini.n的对应权矩阵；D二(d)(i二1,2,.91,n二1,2,3,4,5):与第i个观测站点的距离最小的前5个网格点ini.n的距离矩阵；dx-x)2+(y-y)2():欧拉公式计算第i个观测站点到第j个实测站点的距vijij离计算；p(x,y)(i二1,2,.91):第i个观测站点的坐标点的纬度和经度；iii8:可接受度数差；qin第i个实测站点的距离最小的前5个网格点-!l(i二l,2,.91,n二1,2,3,4,5)：f丄djTij到第i个观测站点的权

11、重;y=fqf(iijininn=1=l,2,.,91,n=1,2,3,4,5):第i个观测站点在第j个时段的预测雨量之和；M=(y)(i=1,2,.91,j=12.164):第i个观测站点在第j个时段的预测雨量矩阵;ijijijI丄)(i=1,2,.91,j=1,2,.164):第i个观测站点在第j个时段的实际雨量矩阵；Q=ijysijyijj(i=12.91,j=1,2,.164):第i个观测站点在第j个时段的预报偏差率;164(q)(i=1,2,.,91,j=1,2,.,164):预报偏差率的算术平方根;i=1j=1J,J,.,J(c=1,2):第c种方法的预报数据与实测数据处在同一级别

12、，相差同一c0c1c6级别，相差1级相差6级的频数五、模型的建立与求解5、1问题一模型建立及其求解由于问题涉及的数据量较大，为了模型计算的简便，需要对数据进行处理。(一)数据处理：为了评价两种6小时雨量预报方法的准确性，我们采用网格点上的预报数据来预测观测站点的数据，再来和实际测得的数据相比，判断其准确性，首先从预测数据选取可接受范围内的网格点，其次选取实测点附近5个距离最小的预测点，计算5个网格点分别到观测站点的权重，并计算观测站点的预测雨量。(1)可接受范围内网格点的选取以坐标p(a,b)(i二1,2,.,91)为基准点，给定一个可接受度数差(在求解中取iii=0.2，对任意的p(i二1,

13、2,.,91),搜索其任一个观测站点在纬度和经度都上下增加的i正方形内所有等距网格点。若网格点x(m,n)纬度和经度在同时满足a-ma+和iib-nb+时，即认为该网格点时可接受范围内的网格点。ii(2)网格点到站点距离的计算找到可接受范围内的网格点后，我们计算网格点和这些观测站点间的距离。再得到个观测站点和等距网格点之间的距离后，将各观测站点按距离从小到大排序后保存，我们用Matlab编程求得结果(程序见附录二)。(3)权的计算各观测站点和等距网格点之间的距离从小到大排序后，为了更好地用网格点的预报雨量来预测观测站点的雨量，我们取前5个到观测站点距离最小的等距网格点。根据欧拉公式距离的倒数加

14、权的方法，先算出前5个网格点到网格点的距离，再分别对它们求倒，则5个网格点分别到观测站点的权重为它们之间距离的倒数。1权的计算公式为：qindj=1ij-in,(i二1,2,.,91,n二1,2,3,4,5)4)观测站点预测雨量的计算为了预测各观测站点在某月某日某个时段的雨量值，我们采用距离的倒数加权的方法，取出5个等距网格点分别在某月某日某个时段的雨量值，然后分别乘以它们各自对观测站点的权重，再求和就为预测降雨量。预测降雨量的计算公式：y=qf；ijininn=1每个时段中，每个观测点对应有5个网格点预测雨量值，可计算出1个观测站点预测雨量值，91个观测站点，164个时段就可计算出一个第i个

15、观测站点的第j个时段预测降雨量矩阵M=(y)(i=j=1,2,.,164)。iji.j这里我们用Matlab编程求得两个方法对应的预测降雨量矩阵(程序见附录三)。(二)模型建立及求解：将两个方法对应的预测降雨量矩阵分别与实测降雨量矩阵进行比较，分析出哪一个的准确性高。这里我们用计算出预报偏差率的算术平方根作为一个准确性指数，来辨别准确性的高低。预测降雨量的矩阵M=(y)(i=j=1,2,.,164)和实际降雨量矩阵iji.jI=(s)(i=1,2,.,91,j=1,2,.,164)。则有:iji.j预报偏差率计算公式：b=九為(i=1,2,.,91,j=1,2,.,164)；ijyij取二sh

16、ice_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:91forj=1:ld(i,j)二sqrt(xb4(i)-xb3(j)厂2+(yb4

17、(i)-yb3(j)厂2);%计算距离endend附录三观测站点预测降雨量的计算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=l:rK(i,l:l)=F(i,:)/max(max(F(i,:);%每个最小值点权的计算endH=K(:,1,2,3,4,5);H(find(isnan(H)=l)=l%将H中NaN替换为1HH二reshape(H,455,l);%前五个最小值的权并排序a=bl(b3(D);%返回91*5个点在x中

18、的矩阵位置%预测数据调用和降雨量的筛选%将CurrentDirectory路径改为数据所在位置%A1=;A2=;formonth1=6:7ifmonth1=6;fordata1=18:28fortime1=1:4str1二fint2str(month1)int2str(data1)int2str(time1)_dis1;str2=fint2str(month1)int2str(data1)int2str(time1)_dis2;tmp1=load(str1);tmp2=load(str2);A1=A1;tmp1(:);A2二A2;tmp2(:);endendelsefordata1=1:30i

19、fdata110strtmp=0int2str(data1);elsestrtmp=int2str(data1);endfortime1=1:4str1=fint2str(month1)strtmpint2str(time1)_dis1;str2=fint2str(month1)strtmpint2str(time1)_dis2;tmp1=load(str1);tmp2=load(str2);A1=A1;tmp1(:);A2=A2;tmp2(:);endendendend%A1为第一种方法统计的所有数据，A2为第二种方法统计的所有数据w二reshape(a,455,l);%将91*5个点转换成

20、455*1All二Al(:,w);%第一类164*455个降雨量A22=A2(:,w);%第二类164*455个降雨量gg=repmat(HH,164,1)%为使HH与All，A22大小相同，扩展平铺HH,gg大小为164*455A1_1=gg.*A11%第一类方法权与降雨量乘积A2_2=gg.*A22%第二类方法权与降雨量乘积大小为164*455%分别对A1_1，A2_2每隔5列求和并返回sum1,sum2%sum1=;fori=1:5:454sum1二sum1A1_1(:,i)+A1_1(:,i+1)+A1_1(:,i+2)+A1_1(:,i+3)+A1_1(:,i+4);endsum2=

21、;fori=1:5:454sum2二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:7ifmonth1=6;fordata1=18:28%str3=020int2str(month1)int2str(data1);str3=020int2str(month1)int2str(data1).SIX;tmp3=load(str3);B1=B1;tmp3(:,4:7);endelsefordata1=1:30ifdata110strtmp=0int2str(data1);elsestrtmp=int2str(data1);endstr3=020int2str(month1)strtmp.SIX;tmp3=load(str3);B1=B1;tmp3(:,4:7);endendendB_1二abs(suml-Bl)./suml;%第一种方法的偏差率B_2=abs(sum2-Bl)./sum2;%第二种方法的偏差率B_l(find(isnan(B_l)=l)=O;%将B_1中NaN替换成为0B_2(find(isnan(B_2)=l)=0;%将B_2中NaN替换成为0B_11