雨量预报分析的评价模型数学建模Word文档下载推荐.docx
《雨量预报分析的评价模型数学建模Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《雨量预报分析的评价模型数学建模Word文档下载推荐.docx（19页珍藏版）》请在冰豆网上搜索。
准确度越高，表明预报准确度越高，相应公众感受就越好。
预报1的预报准确度为%高于预报2的准确度%，公众更易接受第一种预报方法。
关键字：
散乱节点插值残差平方和平均误差预报等级矩阵预测准确度
二、问题重述
雨量预报对农业生产和城市工作和生活有重要作用，但准确、及时地对雨量作出预报是一个十分困难的问题，广受世界各国关注。
我国某地气象台和气象研究所正在研究6小时雨量预报方法，即每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度、北纬32度附近的53×
47的等距网格点上。
同时设立91个观测站点实测这些时段的实际雨量，由于各种条件的限制，站点的设置是不均匀的。
气象部门希望建立一种科学评价预报方法好坏的数学模型与方法。
气象部门提供了41天的用两种不同方法的预报数据和相应的实测数据。
预报数据在文件夹FORECAST中，实测数据在文件夹MEASURING中，其中的文件都可以用Windows系统的“写字板”程序打开阅读。
FORECAST中的文件和分别包含网格点的经纬度，其余文件名为<
f日期i>
_dis1和<
_dis2，例如f6181_dis1中包含2002年6月18日晚上20点采用第一种方法预报的第一时段数据（其2491个数据为该时段各网格点的雨量），而f6183_dis2中包含2002年6月18日晚上20点采用第二种方法预报的第三时段数据。
MEASURING中包含了41个名为<
日期>
.SIX的文件，如表示2002年6月18日晚上21点开始的连续4个时段各站点的实测数据（雨量），这些文件的数据格式是：
站号纬度经度第1段第2段第3段第4段
58138
58139
58141
58143
58146
……
雨量用毫米做单位，小于0.1毫米视为无雨。
（1）请建立数学模型来评价两种6小时雨量预报方法的准确性；
（2）气象部门将6小时降雨量分为6等：
—毫米为小雨，—6毫米为中雨，—12毫米为大雨，—25毫米为暴雨，—60毫米为大暴雨，大于60.1毫米为特大暴雨。
若按此分级向公众预报，如何在评价方法中考虑公众的感受
三、名词和符号说明
表示预测点的经度值
表示预测点的纬度值
表示预测点已知的预测值
表示观察站点的经度值
表示观察站点的纬度值
表示三次样条插值的参数选项
观察站点的预测值
表示两种预测值矩阵中第i个元素的测量值
表示实测矩阵中第i个元素的测量值
表示残差平方和
表示均值误差
表示预测准确度
预报正确次数
预测次数
四、模型假设
：
假设题目中全部数据真实可靠，忽略误差；
假设观测站所在位置的经纬度准确无误；
假设天气预报针对的位置在所给网格点附近；
假设雨量在各网点之间的变动是连续的；
五、问题分析
针对问题1，我们将两种预测方法的所有预测值构造成两个以有序时间段对应的预测值为列,以网格点的个数为行的2491×
164矩阵，对于91观测站点41天的实测值做同样的处理，构造成91×
164的矩阵。
这样，繁琐的数据经过预处理后就整理成了三个矩阵。
由于观测站点相应位置没有两种预测方法对应的预测值，无法直接进行评价，我们采用了三次样条插值的方法进行插值预处理，到了91个观测站点两种预测方法的相应时刻的预测值，然后将两种预测方法雨量预测值与雨量实测值进行比较，从而判断出两种预测方法的准确性。
针对问题2，我们根据要求的雨量分级方式来考虑观众的感受。
我们将问题1中91个观测站点预测处理后雨量预测值构成的两个91×
164矩阵和实际雨量观测值构成的91×
164这个三个矩阵分别采用雨量等级记法构造出三个新的矩阵，然后分别把两个预测值构成的降雨量等级矩阵和观测值构成的等级矩阵对应元素相减并取绝对值，并进行等级统计，再利用预测准确度检验法进行判断，准确度越高说明我们预报的误差越小，表明预测方法更准确。
六、模型建立
1、数据预处理
（1）针对问题1
根据上面的分析，我们先对数据进行预处理。
处理方法为：
把FORECAST文件夹中的第一种和第二种预测方式得到的数据分开两个文件夹，分别以记事本格式按照日期的先后顺序有序的导入Matlab的workspace工作空间中，然后建立m文件编辑公式将两部分的数据导入53×
164的三维矩阵rain1和rain2中；
在命令窗中直接导入预测点的经度和纬度存入矩阵lon和lat中，导入实测点的经度和纬度存入矩阵lon1和lat1中，并对其作图，如图5-1。
实现的matlab语句已呈现在附录中。
图5-1预测点（彩色实线）与实测点（蓝色孤点）
由于三维矩阵无法用表格的形式呈现，我们分别截取了rain1和rain2矩阵的第一层呈现在下表5-1和5-2中，rain3是二维矩阵，将其数据呈现在表5-3中：
表5-1预测方法一在6月18日第一个时间段的预测值构造成的矩阵
预测值
B1
B2
B3
……….
B47
A1
A2
A3
A52
A53
表5-2预测方法二在6月18日第一个时间段的预测值构造成的矩阵
表5-391个站点在41天共164个时段的雨量实测值
D1
D2
D3
D164
P1
P2
P3
P90
P91
（2）针对问题2
由问题1我们可以整理得到91个观测点41天的预测值和测量值对应的两个91×
164矩阵，再根据问题2中气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵，如下表5-4和5-5：
表5-4预测方法一在6月18日第一个时间段的预报等级构造成的矩阵
预报等级
表5-5预测方法二在6月18日第一个时间段的预报等级构造成的矩阵
2、模型建立与求解
由于91个观测站点没有相应的预测值，因此不能够直接对实测值进行评价，属于离散的散乱节点，我们利用插值基点为散乱节点的插值函数
[1]在Matlab中进行三次样条插值处理，插值函数
为：
[1]
其中
表示预测点的经度值，
表示预测点的纬度值，
表示预测点的已知的预测值，
表示观察站点的纬度值，
表示观察站点的经度值，
表示三次样条插值的参数选项，
观察站点的预测值。
将91个观测站点41天164个时段的雨量情况进行预测之后，我们可以建立模型来评价这两种预测方法。
这里我们利用残差平方和
[2]以及平均误差
[3]来作为评价的标准：
[2]
[3]
最后，我们根据
和
的值进行评价，取值越大，表明预报的值准确性越低。
因此，残差平方和
利用公式
（1），我们在Matlab中应用编程求解，程序代码见附录。
求解之后得到91个观测点41天164个时段的预测值，整理成91×
164矩阵，然后把预测矩阵和实测矩阵对应元素值相减取平方作残差平方和，再作平均误差，最后结果如下表5-6：
表5-6雨量预测计算结果数据表
结果
平均误差
预测方法一
预测方法二
由此可知，雨量预报方法一的准确性更高一些。
164矩阵，再根据问题2中气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵，继续利用残差平方和以及平均误差来作为评价的标准。
表5-7雨量等级计算结果数据表
2774
2806
由结果可知，预测方法一的准确度更高。
七、模型优化
针对问题2，由于残差平方和与平均误差难以反映真实汇报的准确度，我们将模型改进优化。
164矩阵，再根据问题2中气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，得到两个预报等级矩阵，将两个预报等级矩阵与实测等级矩阵做差值运算，得到两个等级差矩阵，对等级差作绝对值处理后，我们就可以从中进行等级差统计。
[4]
准确度越高，表明预报越准确度越高，相应公众感受就越好。
雨量分级与统计程序见附录4，数据统计结果见下表7-1：
表7-1两种预报等级与实测等级的级差与比例表
雨量等级差
预报1
等级差个数
等级差个数比例
预报2
12425
12404
1
2122
2128
2
55
52
3
10
4
5
6
由表7-1可见预报1的预报准确度为%高于预报2的准确度%。
预报1与实测雨量的符合度高，那么两种预报方法中，公众更易接受第一种预报方法。
八、模型的评价与推广
优点：
本文思路自然流畅，较多地使用了原始数据，对91个站点的预测值比较准确。
缺点：
认为降雨量在不同的地区之间的变化有连续性，过于理想化。
推广：
本文对降雨量的预测方法也可以运用在工农业生产中其他方面，比如一段时期某地区的气温变化、全国粮食产量预测等。
本文应用的散乱节点插值和残差分析法可以应用在很多方面。
九、附录
附录1：
数据初始化
>
File_FORECAST1=dir（'
预测方法一'
）
File_FORECAST1=
166x1structarraywithfields:
name
date
bytes
isdir
File_FORECAST2=dir（'
预测方法二'
File_FORECAST2=
isdir
File_MEASURING=dir（'
MEASURING'
File_MEASURING=
43x1structarraywithfields:
【】
rain1=zeros（53,47,164）;
forn=1:
164
filename=File_FORECAST1（n+2）.name;
rain1（:
:
n）=importdata（filename）;
end
rain2=zeros（53,47,164）;
filename=File_FORECAST2（n+2）.name;
rain2（:
temp=zeros（91,7,41）;
41
filename=File_MEASURING（n+2）.name;
temp（:
rain3=zeros（91,164）;
i=1;
forj=1:
fork=4:
7
rain3（:
i）=temp（:
k,j）;
i=i+1;
end
附录2：
问题1模型的建立与求解
calc1
pjwc=
ccpfh=
+005
calc2
lat;
lon;
lat1;
lon1;
ccpfh=0;
pjwc=0;
fori=1:
wea=rain1（:
i）;
wear=rain3（:
weap=griddata（lat,lon,wea,lat1,lon1）;
chazhi=wear-weap;
pingfang=chazhi.^2;
juedui=abs（chazhi）;
ccpfh=ccpfh+sum（pingfang,1）;
pjwc=pjwc+sum（juedui）;
pjwc=pjwc/（91*164）;
pjwc
ccpfh
wea=rain2（:
附录3：
问题2模型的建立与求解
calc3
pjwc1=
ccpfh1=
2774
pjwc2=
ccpfh2=
2806
rank1=zeros（91,164）;
rank2=zeros（91,164）;
rank=zeros（91,164）;
ccpfh1=0;
pjwc1=0;
ccpfh2=0;
pjwc2=0;
wea1=rain1（:
wea2=rain2（:
weap1（:
i）=griddata（lat,lon,wea1,lat1,lon1,'
cubic'
）;
weap2（:
i）=griddata（lat,lon,wea2,lat1,lon1,'
forj=1:
91
ifweap1（j,i）<
rank1（j,i）=0;
elseifweap1（j,i）<
=
rank1（j,i）=1;
rank1（j,i）=2;
rank1（j,i）=3;
rank1（j,i）=4;
rank1（j,i）=5;
else
rank1（j,i）=6;
ifweap2（j,i）<
rank2（j,i）=0;
elseifweap2（j,i）<
rank2（j,i）=1;
rank2（j,i）=2;
rank2（j,i）=3;
rank2（j,i）=4;
rank2（j,i）=5;
rank2（j,i）=6;
ifrain3（j,i）<
rank（j,i）=0;
elseifrain3（j,i）<
rank（j,i）=1;
rank（j,i）=2;
rank（j,i）=3;
rank（j,i）=4;
rank（j,i）=5;
rank（j,i）=6;
chazhi1=rank1-rank;
chazhi2=rank2-rank;
pingfang1=chazhi1.^2;
pingfang2=chazhi2.^2;
juedui1=abs（chazhi1）;
juedui2=abs（chazhi2）;
ccpfh1=ccpfh1+sum（sum（pingfang1,1）,2）;
ccpfh2=ccpfh2+sum（sum（pingfang2,1）,2）;
pjwc1=pjwc1+sum（sum（juedui1）,2）;
pjwc2=pjwc2+sum（sum（juedui2）,2）;
pjwc1=pjwc1/（91*164）;
pjwc2=pjwc2/（91*164）;
pjwc1
ccpfh1
pjwc2
ccpfh2
附录4：
问题2的模型优化求解
calc4
count=
1242512404
24282451
6056
1011
12
00
rate=
w1=abs（rank1-rank）;
w2=abs（rank2-rank）;
count=zeros（7,2）;
rate=zeros（7,2）;
ifw1（j,i）==0
count（1,1）=count（1,1）+1;
elseifw1（j,i）==1
count（2,1）=count（2,1）+1;
elseifw1（j,i）==2
count（3,1）=count（3,1）+1;
elseifw1（j,i）==3
count（4,1）=count（4,1）+1;
elseifw1（j,i）==4
count（5,1）=count（5,1）+1;
elseifw1（j,i）==5
count（6,1）=count（6,1）+1;
elseifw1（j,i）==6
count（7,1）=count（7,1）+1;
ifw2（j,i）==0
count（1,2）=count（1,2）+1;
elseifw2（j,i）==1
count（2,2）=count（2,2）+1;
elseifw2（j,i）==2
count（3,2）=count（3,2）+1;
elseifw2（j,i）==3
count（4,2）
升级会员 