搜索
2020/1/10黄建华制作0（五）MATLAB应用举例例5.12005高教社杯全国大学生数学建模竞赛题目C雨量预报方法的评价雨量预报对农业生产和城市工作和生活有重要作用，但准确、及时地对雨量作出预报是一个十分困难的问题，广受世界各国关注。我国某地气象台和气象研究所正在研究6小时雨量预报方法，即每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度、北纬32度附近的53×47的等距网格点上。同时设立91个观测站点实测这些时段的实际雨量，由于各种条件的限制，站点的设置是不均匀的。2020/1/10黄建华制作1例5.1雨量预报方法的评价气象部门希望建立一种科学评价预报方法好坏的数学模型与方法。气象部门提供了41天的用两种不同方法的预报数据和相应的实测数据。预报数据在文件夹FORECAST中，实测数据在文件夹MEASURING中，其中的文件都可以用Windows系统的“写字板”程序打开阅读。FORECAST中的文件lon.dat和lat.dat分别包含网格点的经纬度，其余文件名为f日期i_dis1和f日期i_dis2，例如f6181_dis1中包含2002年6月18日晚上20点采用第一种方法预报的第一时段数据（其2491个数据为该时段各网格点的雨量），而f6183_dis2中包含2002年6月18日晚上20点采用第二种方法预报的第三时段数据。2020/1/10黄建华制作2例5.1雨量预报方法的评价MEASURING中包含了41个名为日期.SIX的文件，如020618.SIX表示2002年6月18日晚上21点开始的连续4个时段各站点的实测数据（雨量），这些文件的数据格式是：站号纬度经度第1段第2段第3段第4段5813832.9833118.51670.00000.200010.10003.10005813933.3000118.85000.00000.00004.60007.40005814133.6667119.26670.00000.00001.10001.40005814333.8000119.80000.00000.00000.00001.80005814633.4833119.81670.00000.00001.50001.9000……雨量用毫米做单位，小于0.1毫米视为无雨。2020/1/10黄建华制作3例5.1雨量预报方法的评价（1）请建立数学模型来评价两种6小时雨量预报方法的准确性；（2）气象部门将6小时降雨量分为6等：0.1—2.5毫米为小雨，2.6—6毫米为中雨，6.1—12毫米为大雨，12.1—25毫米为暴雨，25.1—60毫米为大暴雨，大于60.1毫米为特大暴雨。若按此分级向公众预报，如何在评价方法中考虑公众的感受？（注：本题数据位于压缩文件C2005Data.rar中,可从下载）2020/1/10黄建华制作4例5.1雨量预报方法的评价模型的分析：本题的关键主要是采用Matlab软件对所提供的数据进行分析并得到主要结论。对于问题一，采用load命令和循环结构实现数据文件的载入，再从实测数据文件中提出实测点位置数据，并且依次从预测数据文件中通过曲面拟合命令griddata得到相应日期、时段、方法下的对应位置上的预测估计值。然后分别计算两种预报方法下的实测值和预测值的偏差并且求对应的总偏差平方和，根据两个总偏差平方和的大小来得出两种方法的优劣比较，通过Matlab程序的实现得到结果为：第一种方法比第二种方法好。2020/1/10黄建华制作5例5.1雨量预报方法的评价对于问题二，我们认为公众的感受主要体现在预测与实际之间偏差的大小程度，然而，实际测量值未知，因此我们考虑在分级预报中加入准确概率的方式来实现公众满意度的提高。我们认为雨量的实测值应该在预测值点处服从正态分布，通过合理的假设和推导，我们得到：正态分布的均值可以取为雨量预测值，不同雨量分级区间上的方差可以近似取为对应区间上的总偏差平方和的平均值。运用Matlab软件编程，可以实现对每一个预报数据的预报内容的改变，并且通过几个不同的数据体现程序运行所得到的结果。最后，对模型的缺点进行了讨论和改进。2020/1/10黄建华制作6例5.1雨量预报方法的评价模型的假设假设测量雨量的工具正常不受任何因素的影响，且所得数据真实可靠；测量点所在位置的等距网格点均视为质点；符号说明x(i,j)表示第j种方法下的第i个实测点处的雨量实测值y(i,j)表示第j种方法下的第i个实测点处的雨量预测估计值i=1,2,3…….,91*164,j=1,2e(i,j)表示在第j种方法下的第i个实测点处的偏差Y表示实际值的随机变量y*表示预报值ST2(j)表示在第j种方法下的所有实测点处的总偏差平方和2020/1/10黄建华制作7例5.1雨量预报方法的评价问题一模型分析：考察下载的lon.dat、lat.dat、020日期.SIX及f日期i_disj等数据文件。我们可以看到雨量预报的网格点有53*47个。而雨量实测点只有91个，且分布不均匀，因此我们不能直接引用数据文件进行实测数据与预测数据之间偏差的计算。显然首先要对所提供的数据文件进行处理，由于所引用的文件比较多，而且具有一定的规律性，所以数据文件的载入可用MATLAB命令load和循环结构来实现。当有了具体的实测点位置对应的实测数据yi及预测数据yi*之后,我们可以求出每种方法下所有实测点位置上的总偏差平方和。然后比较两个总偏差平方和的大小，就可以判断两种方法的优劣性。2020/1/10黄建华制作8例5.1雨量预报方法的评价问题一模型建立：通过Matlab程序编程计算，具体过程如下：（1）使用MATLAB命令load和循环结构将lot.dat、lat.dat、f6181_dis1…、f7304_dis2、020618.six…、020730.six等文件输入MATLAB程序中备用，其中020618.six…、020730.six等文件系统自动在前面添上X。（2）从X020618文件中取第二、第三列作为实测点位置。其中第一列对应实测点的纬度，第二列对应实测点的经度。分别用lat、lon中的数据作为网格点的横坐标和纵坐标。依次取f6181_dis1……,f7304_dis1,f6181_dis2,……,f7304_dis2文件对应点数据作为竖坐标，张成一个预报值数据曲面，用曲面拟合命令griddata得出对应实测点处的预报估计值，得到一个2020/1/10黄建华制作9例5.1雨量预报方法的评价91*328维的对应实测点处的预报估计值矩阵yczjz。（3）依次从X020618，…….X020730文件中取第四到第七列组成所有日期，所有时段下的实测数据矩阵sczjz。（4）将sczjz分别与yczjz前164列、后164列对应元素相减，得出第一种方法下和第二种方法下所有实测点处的偏差，组成矩阵f1pcz和f2pcz。（5）分别求矩阵f1pcz和f2pcz所有元素的平方之和，赋予变量pcpfh(1)和pcpfh(2)。（6）比较变量pcpfh(1)、pcpfh(2)的大小，得出两种方法优劣的比较。各过程的MATLAB程序实现，可分别由M-文件sjzr.m、qsj.m、pcbj.m运行得到。2020/1/10黄建华制作10例5.1雨量预报方法的评价程序sjzr.m：rqjz=[618:628,701:730];%产生1*41维日期行向量loadlat.DAT;%载入53*47维纬度矩阵loadlon.DAT;%载入53*47维经度矩阵forrq=rqjz%用日期行向量作为循环向量load(['020',int2str(rq),'.SIX'])%载入53*47维实测值矩阵,系统自动在文件名前加大写字母Xforsd=1:4%指定时段循环向量forff=1:2%指定方法循环向量load(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)])%载入53*47维预测值矩阵，共328个endendend2020/1/10黄建华制作11例5.1雨量预报方法的评价将文件全部拷贝在默认的work文件加下，运行程序sjzr.m：2020/1/10黄建华制作12例5.1雨量预报方法的评价2020/1/10黄建华制作13例5.1雨量预报方法的评价程序huatu.m%画一个时点预测曲面图与实测散点图比较clf;%预先清除别的图象Z=eval(['f',int2str(618),int2str(1),'_dis',int2str(1)]);%取出某时段预报值作为竖坐标mesh(lat,lon,Z)%张成预测曲面axis([27,36,117,125,0,0.4])shadingflatxlabel('纬度');ylabel('经度');title('图象比较');holdonX1=X020618(:,2);Y1=X020618(:,3);Z1=X020618(:,4);%取出实测点数据plot3(X1,Y1,Z1,'r*')%画出实测点数据散点图2020/1/10黄建华制作14例5.1雨量预报方法的评价2020/1/10黄建华制作15例5.1雨量预报方法的评价程序qsj.mycr=1;%指定预测值矩阵列标yczjz=zeros(91,328);%预先指定预测值矩阵为全零矩阵forff=1:2%指定方法循环变量forrq=rqjz%指定日期循环变量forsd=1:4%指定时段循环变量z0=eval(['f',int2str(rq),int2str(sd),'_dis',int2str(ff)]);%依次取出雨量预测值数据yczjz(:,ycr)=griddata(lat,lon,z0,X020618(:,2),X020618(:,3));%在张成的预测数据曲面上拟合对应实测点的雨量预测值，依次放入预测值矩阵对应列ycr=ycr+1;%预测值矩阵列标向后一列end2020/1/10黄建华制作16例5.1雨量预报方法的评价endendscr=1;%指定实测值矩阵列标sczjz=zeros(91,164);%预先指定实测值矩阵为全零矩阵forrq=rqjz%指定日期循环变量sczjz(:,scr:scr+3)=eval(['X020',int2str(rq),'(:,4:7)']);%依次取出雨量实测值数据放入实测值矩阵scr=scr+4;%实测值矩阵列标向后4列end2020/1/10黄建华制作17例5.1雨量预报方法的评价程序pcbj.mff1=yczjz(:,1:164);%取出预测值矩阵的前164列ff2=yczjz(:,165:328);%取出预测值矩阵的后164列f1pcz=sczjz-ff1;%计算方法1下的偏差f2pcz=sczjz-ff2;%计算方法2下的偏差pcpfh=[0,0];%预先产生1*2维偏差平方和零矩阵forpch=1:91%指定偏差行循环变量forpcr=1:164%指定偏差列循环变量pcpfh(1)=pcpfh(1)+f1pcz(pch,pcr).^2;%计算方法1偏差平方和pcpfh(2)=pcpfh(2)+f2pcz(pch,pcr).^2;%计算方法2偏差平方和2020/1/10黄建华制作18例5.1雨量预报方法的评价endendpcpfhifpcpfh(1)==pcpfh(2)%比较两种方法下的偏差平方和display('一样好')elseifpcpfh(1)pcpfh(2)display('第一种方法好')elsedisplay('第二种方法好')end首先将数据文件减压缩在matlab程序的work文件夹下，然后依次运行M-文件sjzr.m、qsj.m、pcbj.m，得到如下结论：偏差平方和：226730243880第一种方法比第二种方法好。2020/