雨量预报分析的评价模型数学建模.docx

资源描述

雨量预报分析的评价模型数学建模.docx

《雨量预报分析的评价模型数学建模.docx》由会员分享，可在线阅读，更多相关《雨量预报分析的评价模型数学建模.docx（19页珍藏版）》请在冰豆网上搜索。

雨量预报分析的评价模型数学建模.docx

雨量预报分析的评价模型数学建模

雨量预报分析的评价模型

一、摘要

我们将FORECAST文件夹中的数据按日期先后顺序导入Matlab，建立53×47×164的三维矩阵rain1和rain2；把MEASURING文件夹中的数据以同样方法导入91×7×41的三维矩阵temp中，然后建立循环将temp矩阵中每一层的后4列提取，另存入一个91×164的rain3矩阵；在命令窗中直接导入预测点的经度和纬度存入矩阵lon和lat中，导入实测点的经度和纬度存入矩阵lon1和lat1中，并对其作图，得到实测点和预测点的经纬度图。

整理得到91个观测点41天的预测值和测量值对应的两个91×164矩阵，根据气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵。

针对问题

（1），利用插值基点为散乱节点的插值函数

[1]在Matlab中进行三次样条插值处理，将91个观测站点41天164个时段的雨量情况进行预测。

利用残差平方和

以及平均误差

来作为评价的标准。

残差平方和

与平均误差

值较小的一种预测方法作为较好的预报方法。

残差平方和以及平均误差数值越小，表明预报越准确度越高。

预测方法一的残差平方和为，平均误差为。

预测方法二的残差平方和为，平均误差为。

雨量预报方法一的准确性更高一些。

针对问题

（2），两个预报等级矩阵，继续利用残差平方和以及平均误差来作为评价的标准。

残差平方和以及平均误差数值越小，表明预报越准确度越高，相应公众感受就越好。

预测方法一的残差平方和为2774，平均误差为。

预测方法二的残差平方和为2806，平均误差为。

雨量预报方法一的准确性更高一些。

由于残差平方和与平均误差难以反映真实汇报的准确度，我们将模型改进优化。

把矩阵中相应的降雨量值转化为其所对应等级值，得到两个预报等级矩阵，将两个预报等级矩阵与实测等级矩阵做差值运算，得到两个等级差矩阵，对等级差作绝对值处理，进行等级差统计。

我们利用预测准确度检验法对两种预报进行评价。

预测准确度（

）等于预报正确次数（

）（即运算之差为0的情况）和预测次数（

）之比，即

。

准确度越高，表明预报准确度越高，相应公众感受就越好。

预报1的预报准确度为%高于预报2的准确度%，公众更易接受第一种预报方法。

关键字：

散乱节点插值残差平方和平均误差预报等级矩阵预测准确度

二、问题重述

雨量预报对农业生产和城市工作和生活有重要作用，但准确、及时地对雨量作出预报是一个十分困难的问题，广受世界各国关注。

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

同时设立91个观测站点实测这些时段的实际雨量，由于各种条件的限制，站点的设置是不均匀的。

气象部门希望建立一种科学评价预报方法好坏的数学模型与方法。

气象部门提供了41天的用两种不同方法的预报数据和相应的实测数据。

预报数据在文件夹FORECAST中，实测数据在文件夹MEASURING中，其中的文件都可以用Windows系统的“写字板”程序打开阅读。

FORECAST中的文件和分别包含网格点的经纬度，其余文件名为_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毫米视为无雨。

请建立数学模型来评价两种6小时雨量预报方法的准确性；

气象部门将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×47×164的三维矩阵rain1和rain2中；把MEASURING文件夹中的数据以同样方法导入91×7×41的三维矩阵temp中，然后建立循环将temp矩阵中每一层的后4列提取，另存入一个91×164的rain3矩阵；在命令窗中直接导入预测点的经度和纬度存入矩阵lon和lat中，导入实测点的经度和纬度存入矩阵lon1和lat1中，并对其作图，如图5-1。

实现的matlab语句已呈现在附录中。

图5-1预测点（彩色实线）与实测点（蓝色孤点）

由于三维矩阵无法用表格的形式呈现，我们分别截取了rain1和rain2矩阵的第一层呈现在下表5-1和5-2中，rain3是二维矩阵，将其数据呈现在表5-3中：

表5-1预测方法一在6月18日第一个时间段的预测值构造成的矩阵

预测值

……….

B47

……….

A52

……….

A53

……….

表5-2预测方法二在6月18日第一个时间段的预测值构造成的矩阵

预测值

……….

B47

……….

A52

……….

A53

……….

表5-391个站点在41天共164个时段的雨量实测值

预测值

……….

D164

……….

P90

……….

P91

……….

（2）针对问题2

由问题1我们可以整理得到91个观测点41天的预测值和测量值对应的两个91×164矩阵，再根据问题2中气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵，如下表5-4和5-5：

表5-4预测方法一在6月18日第一个时间段的预报等级构造成的矩阵

预报等级

……….

B47

……….

A52

……….

A53

……….

表5-5预测方法二在6月18日第一个时间段的预报等级构造成的矩阵

预报等级

……….

B47

……….

A52

……….

A53

……….

2、模型建立与求解

（1）针对问题1

由于91个观测站点没有相应的预测值，因此不能够直接对实测值进行评价，属于离散的散乱节点，我们利用插值基点为散乱节点的插值函数

[1]在Matlab中进行三次样条插值处理，插值函数

为：

[1]

其中

表示预测点的经度值，

表示预测点的纬度值，

表示预测点的已知的预测值，

表示观察站点的纬度值，

表示观察站点的经度值，

表示三次样条插值的参数选项，

观察站点的预测值。

将91个观测站点41天164个时段的雨量情况进行预测之后，我们可以建立模型来评价这两种预测方法。

这里我们利用残差平方和

[2]以及平均误差

[3]来作为评价的标准：

[2]

[3]

最后，我们根据

和

的值进行评价，取值越大，表明预报的值准确性越低。

因此，残差平方和

与平均误差

值较小的一种预测方法作为较好的预报方法。

利用公式

（1），我们在Matlab中应用编程求解，程序代码见附录。

求解之后得到91个观测点41天164个时段的预测值，整理成91×164矩阵，然后把预测矩阵和实测矩阵对应元素值相减取平方作残差平方和，再作平均误差，最后结果如下表5-6：

表5-6雨量预测计算结果数据表

结果

残差平方和

平均误差

预测方法一

预测方法二

由此可知，雨量预报方法一的准确性更高一些。

（2）针对问题2

由问题1我们可以整理得到91个观测点41天的预测值和测量值对应的两个91×164矩阵，再根据问题2中气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵，继续利用残差平方和以及平均误差来作为评价的标准。

残差平方和以及平均误差数值越小，表明预报越准确度越高，相应公众感受就越好。

表5-7雨量等级计算结果数据表

结果

残差平方和

平均误差

预测方法一

2774

预测方法二

2806

由结果可知，预测方法一的准确度更高。

七、模型优化

针对问题2，由于残差平方和与平均误差难以反映真实汇报的准确度，我们将模型改进优化。

由问题1我们可以整理得到91个观测点41天的预测值和测量值对应的两个91×164矩阵，再根据问题2中气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，得到两个预报等级矩阵，将两个预报等级矩阵与实测等级矩阵做差值运算，得到两个等级差矩阵，对等级差作绝对值处理后，我们就可以从中进行等级差统计。

我们利用预测准确度检验法对两种预报进行评价。

预测准确度（

）等于预报正确次数（

）（即运算之差为0的情况）和预测次数（

）之比，即

[4]

准确度越高，表明预报越准确度越高，相应公众感受就越好。

雨量分级与统计程序见附录4，数据统计结果见下表7-1：

表7-1两种预报等级与实测等级的级差与比例表

雨量等级差

预报1

等级差个数

预报1

等级差个数比例

预报2

等级差个数

预报2

等级差个数比例

12425

12404

2122

2128

由表7-1可见预报1的预报准确度为%高于预报2的准确度%。

预报1与实测雨量的符合度高，那么两种预报方法中，公众更易接受第一种预报方法。

八、模型的评价与推广

优点：

本文思路自然流畅，较多地使用了原始数据，对91个站点的预测值比较准确。

缺点：

认为降雨量在不同的地区之间的变化有连续性，过于理想化。

推广：

本文对降雨量的预测方法也可以运用在工农业生产中其他方面，比如一段时期某地区的气温变化、全国粮食产量预测等。

本文应用的散乱节点插值和残差分析法可以应用在很多方面。

九、附录

附录1：

数据初始化

>>File_FORECAST1=dir（'预测方法一'）

File_FORECAST1=

166x1structarraywithfields:

name

date

bytes

isdir

>>File_FORECAST2=dir（'预测方法二'）

File_FORECAST2=

166x1structarraywithfields:

name

date

bytes

isdir

>>File_MEASURING=dir（'MEASURING'）

File_MEASURING=

43x1structarraywithfields:

name

date

bytes

isdir

【】

rain1=zeros（53,47,164）;

forn=1:

164

filename=File_FORECAST1（n+2）.name;

rain1（:

n）=importdata（filename）;

end

【】

rain2=zeros（53,47,164）;

forn=1:

164

filename=File_FORECAST2（n+2）.name;

rain2（:

n）=importdata（filename）;

end

【】

temp=zeros（91,7,41）;

forn=1:

filename=File_MEASURING（n+2）.name;

temp（:

n）=importdata（filename）;

end

rain3=zeros（91,164）;

i=1;

forj=1:

fork=4:

rain3（:

i）=temp（:

k,j）;

i=i+1;

end

附录2：

问题1模型的建立与求解

>>calc1

pjwc=

ccpfh=

+005

>>calc2

pjwc=

ccpfh=

+005

【】

lat;

lon;

lat1;

lon1;

ccpfh=0;

pjwc=0;

fori=1:

164

wea=rain1（:

i）;

wear=rain3（:

i）;

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）;

end

pjwc=pjwc/（91*164）;

pjwc

ccpfh

【】

lat;

lon;

lat1;

lon1;

ccpfh=0;

pjwc=0;

fori=1:

164

wea=rain2（:

i）;

wear=rain3（:

i）;

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）;

end

pjwc=pjwc/（91*164）;

pjwc

ccpfh

附录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;

fori=1:

164

wea1=rain1（:

i）;

wea2=rain2（:

i）;

wear=rain3（:

i）;

weap1（:

i）=griddata（lat,lon,wea1,lat1,lon1,'cubic'）;

weap2（:

i）=griddata（lat,lon,wea2,lat1,lon1,'cubic'）;

end

fori=1:

164

forj=1:

ifweap1（j,i）<

rank1（j,i）=0;

elseifweap1（j,i）<=

rank1（j,i）=1;

elseifweap1（j,i）<=

rank1（j,i）=2;

elseifweap1（j,i）<=

rank1（j,i）=3;

elseifweap1（j,i）<=

rank1（j,i）=4;

elseifweap1（j,i）<=

rank1（j,i）=5;

else

rank1（j,i）=6;

end

fori=1:

164

forj=1:

ifweap2（j,i）<

rank2（j,i）=0;

elseifweap2（j,i）<=

rank2（j,i）=1;

elseifweap2（j,i）<=

rank2（j,i）=2;

elseifweap2（j,i）<=

rank2（j,i）=3;

elseifweap2（j,i）<=

rank2（j,i）=4;

elseifweap2（j,i）<=

rank2（j,i）=5;

else

rank2（j,i）=6;

end

fori=1:

164

forj=1:

ifrain3（j,i）<

rank（j,i）=0;

elseifrain3（j,i）<=

rank（j,i）=1;

elseifrain3（j,i）<=

rank（j,i）=2;

elseifrain3（j,i）<=

rank（j,i）=3;

elseifrain3（j,i）<=

rank（j,i）=4;

elseifrain3（j,i）<=

rank（j,i）=5;

else

rank（j,i）=6;

end

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

rate=

【】

w1=abs（rank1-rank）;

w2=abs（rank2-rank）;

count=zeros（7,2）;

rate=zeros（7,2）;

fori=1:

164

forj=1:

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;

end

fori=1:

164

forj=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）=count（4,2）+1;

elseifw2（j,i）==4

count（5,2）=count（5,2）+1;

elseifw2（j,i

展开阅读全文

雨量预报分析的评价模型 数学建模.docx

雨量预报分析的评价模型数学建模.docx