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

1、雨量预报分析的评价模型数学建模雨量预报分析的评价模型一、摘要我们将FORECAST文件夹中的数据按日期先后顺序导入Matlab，建立5347164的三维矩阵rain1和rain2；把MEASURING文件夹中的数据以同样方法导入91741的三维矩阵temp中，然后建立循环将temp矩阵中每一层的后4列提取，另存入一个91164的rain3矩阵；在命令窗中直接导入预测点的经度和纬度存入矩阵lon和lat中，导入实测点的经度和纬度存入矩阵lon1和lat1中，并对其作图，得到实测点和预测点的经纬度图。整理得到91个观测点41天的预测值和测量值对应的两个91164矩阵，根据气象部门将降雨的等级分为6

2、个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵。针对问题（1），利用插值基点为散乱节点的插值函数 1在Matlab中进行三次样条插值处理，将91个观测站点41天164个时段的雨量情况进行预测。利用残差平方和以及平均误差来作为评价的标准。残差平方和和平均误差值较小的一种预测方法作为较好的预报方法。残差平方和以及平均误差数值越小，表明预报越准确度越高。预测方法一的残差平方和为174290.00，平均误差为0.4553。预测方法二的残差平方和为195580.00，平均误差为0.4753。雨量预报方法一的准确性更高一些。针对问题（2），两个预

3、报等级矩阵，继续利用残差平方和以及平均误差来作为评价的标准。残差平方和以及平均误差数值越小，表明预报越准确度越高，相应公众感受就越好。预测方法一的残差平方和为2774，平均误差为0.1730。预测方法二的残差平方和为2806，平均误差为0.1745。雨量预报方法一的准确性更高一些。由于残差平方和和平均误差难以反映真实汇报的准确度，我们将模型改进优化。把矩阵中相应的降雨量值转化为其所对应等级值，得到两个预报等级矩阵，将两个预报等级矩阵和实测等级矩阵做差值运算，得到两个等级差矩阵，对等级差作绝对值处理，进行等级差统计。我们利用预测准确度检验法对两种预报进行评价。预测准确度（）等于预报正确次数（）（

4、即运算之差为0的情况）和预测次数（）之比，即。准确度越高，表明预报准确度越高，相应公众感受就越好。预报1的预报准确度为83.26%高于预报2的准确度83.11%，公众更易接受第一种预报方法。关键字：散乱节点插值 残差平方和 平均误差 预报等级矩阵 预测准确度二、问题重述雨量预报对农业生产和城市工作和生活有重要作用，但准确、及时地对雨量作出预报是一个十分困难的问题，广受世界各国关注。我国某地气象台和气象研究所正在研究6小时雨量预报方法，即每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度、北纬32度

5、附近的5347的等距网格点上。同时设立91个观测站点实测这些时段的实际雨量，由于各种条件的限制，站点的设置是不均匀的。气象部门希望建立一种科学评价预报方法好坏的数学模型和方法。气象部门提供了41天的用两种不同方法的预报数据和相应的实测数据。预报数据在文件夹FORECAST中，实测数据在文件夹MEASURING中，其中的文件都可以用Windows系统的“写字板”程序打开阅读。FORECAST中的文件lon.dat和lat.dat分别包含网格点的经纬度，其余文件名为_dis1和_dis2，例如f6181_dis1中包含2002年6月18日晚上20点采用第一种方法预报的第一时段数据（其2491个数据

6、为该时段各网格点的雨量），而f6183_dis2中包含2002年6月18日晚上20点采用第二种方法预报的第三时段数据。MEASURING中包含了41个名为.SIX的文件，如020618.SIX表示2002年6月18日晚上21点开始的连续4个时段各站点的实测数据（雨量），这些文件的数据格式是：站号 纬度 经度 第1段 第2段 第3段 第4段 58138 32.9833 118.5167 0.0000 0.2000 10.1000 3.100058139 33.3000 118.8500 0.0000 0.0000 4.6000 7.400058141 33.6667 119.2667 0.000

7、0 0.0000 1.1000 1.400058143 33.8000 119.8000 0.0000 0.0000 0.0000 1.800058146 33.4833 119.8167 0.0000 0.0000 1.5000 1.9000雨量用毫米做单位，小于0.1毫米视为无雨。(1)请建立数学模型来评价两种6小时雨量预报方法的准确性；(2)气象部门将6小时降雨量分为6等：0.12.5毫米为小雨，2.66毫米为中雨，6.112毫米为大雨，12.125毫米为暴雨，25.160毫米为大暴雨，大于60.1毫米为特大暴雨。若按此分级向公众预报，如何在评价方法中考虑公众的感受？三、名词和符号说明表

8、示预测点的经度值表示预测点的纬度值表示预测点已知的预测值表示观察站点的经度值表示观察站点的纬度值表示三次样条插值的参数选项观察站点的预测值表示两种预测值矩阵中第i个元素的测量值表示实测矩阵中第i个元素的测量值表示残差平方和表示均值误差 表示预测准确度预报正确次数预测次数四、模型假设：假设题目中全部数据真实可靠，忽略误差；：假设观测站所在位置的经纬度准确无误；：假设天气预报针对的位置在所给网格点附近；：假设雨量在各网点之间的变动是连续的；五、问题分析针对问题1，我们将两种预测方法的所有预测值构造成两个以有序时间段对应的预测值为列,以网格点的个数为行的2491164矩阵，对于91观测站点41天的实

9、测值做同样的处理，构造成91164的矩阵。这样，繁琐的数据经过预处理后就整理成了三个矩阵。由于观测站点相应位置没有两种预测方法对应的预测值，无法直接进行评价，我们采用了三次样条插值的方法进行插值预处理，到了91个观测站点两种预测方法的相应时刻的预测值，然后将两种预测方法雨量预测值和雨量实测值进行比较，从而判断出两种预测方法的准确性。针对问题2，我们根据要求的雨量分级方式来考虑观众的感受。我们将问题1中91个观测站点预测处理后雨量预测值构成的两个91164矩阵和实际雨量观测值构成的91164这个三个矩阵分别采用雨量等级记法构造出三个新的矩阵，然后分别把两个预测值构成的降雨量等级矩阵和观测值构成的

10、等级矩阵对应元素相减并取绝对值，并进行等级统计，再利用预测准确度检验法进行判断，准确度越高说明我们预报的误差越小，表明预测方法更准确。六、模型建立1、数据预处理（1）针对问题1根据上面的分析，我们先对数据进行预处理。处理方法为：把FORECAST文件夹中的第一种和第二种预测方式得到的数据分开两个文件夹，分别以记事本格式按照日期的先后顺序有序的导入Matlab的workspace工作空间中，然后建立m文件编辑公式将两部分的数据导入5347164的三维矩阵rain1和rain2中；把MEASURING文件夹中的数据以同样方法导入91741的三维矩阵temp中，然后建立循环将temp矩阵中每一层的后

11、4列提取，另存入一个91164的rain3矩阵；在命令窗中直接导入预测点的经度和纬度存入矩阵lon和lat中，导入实测点的经度和纬度存入矩阵lon1和lat1中，并对其作图，如图5-1。实现的matlab语句已呈现在附录2.1中。图5-1 预测点（彩色实线）和实测点（蓝色孤点）由于三维矩阵无法用表格的形式呈现，我们分别截取了rain1和rain2矩阵的第一层呈现在下表5-1和5-2中，rain3是二维矩阵，将其数据呈现在表5-3中：表5-1 预测方法一在6月18日第一个时间段的预测值构造成的矩阵预测值B1B2B3.B47A10.03100.03210.0313.0.0081A20.03020.

12、03100.0315.0.0078A30.03400.03410.0337.0.0079.A520.00880.00890.0084.0.0061A530.00890.00860.0086.0.0062表5-2 预测方法二在6月18日第一个时间段的预测值构造成的矩阵预测值B1B2B3.B47A10.03190.02930.0322.0.0077A20.02880.03110.0303.0.0085A30.03250.03230.0313.0.0082.A520.00950.00890.0079.0.0056A530.00850.00940.0093.0.0057表5-3 91个站点在41天共1

13、64个时段的雨量实测值预测值D1D2D3.D164P1000.0P2000.0P3000.0.P90004.1.0P91000.2.0（2）针对问题2由问题1我们可以整理得到91个观测点41天的预测值和测量值对应的两个91164矩阵，再根据问题2中气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵，如下表5-4和5-5：表5-4 预测方法一在6月18日第一个时间段的预报等级构造成的矩阵预报等级B1B2B3.B47A1000.0A2000.0A3000.0.A52000.0A53000.0表5-5 预测方法二在6月

14、18日第一个时间段的预报等级构造成的矩阵预报等级B1B2B3.B47A1000.0A2000.0A3000.0.A52000.0A53000.02、模型建立和求解（1）针对问题1由于91个观测站点没有相应的预测值，因此不能够直接对实测值进行评价，属于离散的散乱节点，我们利用插值基点为散乱节点的插值函数1在Matlab中进行三次样条插值处理，插值函数为： 1 其中表示预测点的经度值，表示预测点的纬度值，表示预测点的已知的预测值，表示观察站点的纬度值，表示观察站点的经度值，表示三次样条插值的参数选项，观察站点的预测值。将91个观测站点41天164个时段的雨量情况进行预测之后，我们可以建立模型来评价

15、这两种预测方法。这里我们利用残差平方和2以及平均误差3来作为评价的标准： 2 3最后，我们根据和的值进行评价，取值越大，表明预报的值准确性越低。因此，残差平方和和平均误差值较小的一种预测方法作为较好的预报方法。利用公式（1），我们在Matlab中应用编程求解，程序代码见附录2.2。求解之后得到91个观测点41天164个时段的预测值，整理成91164矩阵，然后把预测矩阵和实测矩阵对应元素值相减取平方作残差平方和，再作平均误差，最后结果如下表5-6：表5-6 雨量预测计算结果数据表结果残差平方和平均误差预测方法一174290.000.4553预测方法二195580.000.4753由此可知，雨量预

16、报方法一的准确性更高一些。（2）针对问题2由问题1我们可以整理得到91个观测点41天的预测值和测量值对应的两个91164矩阵，再根据问题2中气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵，继续利用残差平方和以及平均误差来作为评价的标准。残差平方和以及平均误差数值越小，表明预报越准确度越高，相应公众感受就越好。表5-7 雨量等级计算结果数据表结果残差平方和平均误差预测方法一27740.1730预测方法二28060.1745由结果可知，预测方法一的准确度更高。七、模型优化针对问题2，由于残差平方和和平均误差难以反

17、映真实汇报的准确度，我们将模型改进优化。由问题1我们可以整理得到91个观测点41天的预测值和测量值对应的两个91164矩阵，再根据问题2中气象部门将降雨的等级分为6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，得到两个预报等级矩阵，将两个预报等级矩阵和实测等级矩阵做差值运算，得到两个等级差矩阵，对等级差作绝对值处理后，我们就可以从中进行等级差统计。我们利用预测准确度检验法对两种预报进行评价。预测准确度（）等于预报正确次数（）（即运算之差为0的情况）和预测次数（）之比，即 4准确度越高，表明预报越准确度越高，相应公众感受就越好。雨量分级和统计程序见附录4，数据统计结果见下表7-1：表

18、7-1 两种预报等级和实测等级的级差和比例表雨量等级差预报1等级差个数预报1等级差个数比例预报2等级差个数预报2等级差个数比例0124250.8326124040.8311121220.142221280.14262550.0037520.00353100.0007100.0007410.000120.00015000060000由表7-1可见预报1的预报准确度为83.26%高于预报2的准确度83.11%。预报1和实测雨量的符合度高，那么两种预报方法中，公众更易接受第一种预报方法。八、模型的评价和推广优点：本文思路自然流畅，较多地使用了原始数据，对91个站点的预测值比较准确。缺点：认为降雨量在

19、不同的地区之间的变化有连续性，过于理想化。推广：本文对降雨量的预测方法也可以运用在工农业生产中其他方面，比如一段时期某地区的气温变化、全国粮食产量预测等。本文应用的散乱节点插值和残差分析法可以应用在很多方面。九、附录附录1：数据初始化 File_FORECAST1=dir(预测方法一)File_FORECAST1 =166x1 struct array with fields: name date bytes isdir File_FORECAST2=dir(预测方法二)File_FORECAST2 = 166x1 struct array with fields: name date byt

20、esisdir File_MEASURING=dir(MEASURING)File_MEASURING = 43x1 struct array with fields: name date bytes isdir【mypro1.m】rain1=zeros(53,47,164);for n=1:164 filename=File_FORECAST1(n+2).name; rain1(:,:,n)=importdata(filename);end【mypro2.m】rain2=zeros(53,47,164);for n=1:164 filename=File_FORECAST2(n+2).nam

21、e; rain2(:,:,n)=importdata(filename);end【mypro3.m】temp=zeros(91,7,41);for n=1:41 filename=File_MEASURING(n+2).name; temp(:,:,n)=importdata(filename);endrain3=zeros(91,164);i=1;for j=1:41 for k=4:7 rain3(:,i)=temp(:,k,j); i=i+1; endend附录2：问题1模型的建立和求解 calc1pjwc = 0.4553ccpfh = 1.7429e+005 calc2pjwc =

22、0.4753ccpfh = 1.9558e+005【calc1.m】lat;lon;lat1;lon1;ccpfh=0;pjwc=0;for i=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);endpjwc=pjwc/(91*164);pjwcccpfh【calc2.m】

23、lat;lon;lat1;lon1;ccpfh=0;pjwc=0;for i=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);endpjwc=pjwc/(91*164);pjwcccpfh附录3：问题2模型的建立和求解 calc3pjwc1 = 0.1730ccpfh1 =

24、 2774pjwc2 = 0.1745ccpfh2 = 2806【calc3.m】rank1=zeros(91,164);rank2=zeros(91,164);rank=zeros(91,164);ccpfh1=0;pjwc1=0;ccpfh2=0;pjwc2=0;for i=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

25、);endfor i=1:164 for j=1:91 if weap1(j,i)0.1 rank1(j,i)=0; elseif weap1(j,i)=2.5 rank1(j,i)=1; elseif weap1(j,i)=6.0 rank1(j,i)=2; elseif weap1(j,i)=12.0 rank1(j,i)=3; elseif weap1(j,i)=25.0 rank1(j,i)=4; elseif weap1(j,i)=60.0 rank1(j,i)=5; else rank1(j,i)=6; end endendfor i=1:164 for j=1:91 if wea

26、p2(j,i)0.1 rank2(j,i)=0; elseif weap2(j,i)=2.5 rank2(j,i)=1; elseif weap2(j,i)=6.0 rank2(j,i)=2; elseif weap2(j,i)=12.0 rank2(j,i)=3; elseif weap2(j,i)=25.0 rank2(j,i)=4; elseif weap2(j,i)=60.0 rank2(j,i)=5; else rank2(j,i)=6; end endendfor i=1:164 for j=1:91 if rain3(j,i)0.1 rank(j,i)=0; elseif rai

27、n3(j,i)=2.5 rank(j,i)=1; elseif rain3(j,i)=6.0 rank(j,i)=2; elseif rain3(j,i)=12.0 rank(j,i)=3; elseif rain3(j,i)=25.0 rank(j,i)=4; elseif rain3(j,i)=60.0 rank(j,i)=5; else rank(j,i)=6; end endendchazhi1=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);pjwc1ccpfh1pjwc2ccpfh2附录4：问题2的模型优化求解