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

雨量预报分析的评价模型数学建模 Pleasure Group Office【T985AB-B866SYT-B182C-BS682T-STT18】雨量预报分析的评价模型 一、摘要 我们将 FORECAST 文件夹中的数据按日期先后顺序导入 Matlab，建立5347164 的三维矩阵 rain1 和 rain2；把 MEASURING 文件夹中的数据以同样方法导入 91741的三维矩阵 temp中，然后建立循环将 temp矩阵中每一层的后 4 列提取，另存入一个 91164的 rain3 矩阵；在命令窗中直接导入预测点的经度和纬度存入矩阵 lon 和 lat中，导入实测点的经度和纬度存入矩阵 lon1 和lat1中，并对其作图，得到实测点和预测点的经纬度图。整理得到 91个观测点 41天的预测值和测量值对应的两个 91164 矩阵，根据气象部门将降雨的等级分为 6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵。针对问题（1），利用插值基点为散乱节点的插值函数griddata 1在Matlab 中进行三次样条插值处理，将 91 个观测站点 41 天 164 个时段的雨量情况进行预测。利用残差平方和21()nijiiweapwear以及平均误差11nijiiavgweapwearn来作为评价的标准。残差平方和与平均误差avg值较小的一种预测方法作为较好的预报方法。残差平方和以及平均误差数值越小，表明预报越准确度越高。预测方法一的残差平方和为，平均误差为。预测方法二的残差平方和为，平均误差为。雨量预报方法一的准确性更高一些。针对问题（2），两个预报等级矩阵，继续利用残差平方和以及平均误差来作为评价的标准。残差平方和以及平均误差数值越小，表明预报越准确度越高，相应公众感受就越好。预测方法一的残差平方和为 2774，平均误差为。预测方法二的残差平方和为 2806，平均误差为。雨量预报方法一的准确性更高一些。由于残差平方和与平均误差难以反映真实汇报的准确度，我们将模型改进优化。把矩阵中相应的降雨量值转化为其所对应等级值，得到两个预报等级矩阵，将两个预报等级矩阵与实测等级矩阵做差值运算，得到两个等级差矩阵，对等级差作绝对值处理，进行等级差统计。我们利用预测准确度检验法对两种预报进行评价。预测准确度（H）等于预报正确次数（R）（即运算之差为 0的情况）和预测次数（T）之比，即100%RHT。准确度越高，表明预报准确度越高，相应公众感受就越好。预报 1的预报准确度为%高于预报 2 的准确度%，公众更易接受第一种预报方法。关键字：散乱节点插值 残差平方和 平均误差 预报等级矩阵 预测准确度 二、问题重述 雨量预报对农业生产和城市工作和生活有重要作用，但准确、及时地对雨量作出预报是一个十分困难的问题，广受世界各国关注。我国某地气象台和气象研究所正在研究 6小时雨量预报方法，即每天晚上 20点预报从 21点开始的 4个时段（21点至次日 3 点，次日 3 点至 9 点，9 点至 15点，15点至 21点）在某些位置的雨量，这些位置位于东经 120度、北纬 32度附近的 5347 的等距网格点上。同时设立 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毫米视为无雨。(1)请建立数学模型来评价两种 6小时雨量预报方法的准确性；(2)气象部门将 6 小时降雨量分为 6 等：毫米为小雨，6 毫米为中雨，12毫米为大雨，25毫米为暴雨，60 毫米为大暴雨，大于 60.1 毫米为特大暴雨。若按此分级向公众预报，如何在评价方法中考虑公众的感受 三、名词和符号说明 表示预测点的经度值 表示预测点的纬度值 表示预测点已知的预测值 表示观察站点的经度值 表示观察站点的纬度值 表示三次样条插值的参数选项 观察站点的预测值 表示两种预测值矩阵中第 i 个元素的测量值 表示实测矩阵中第 i 个元素的测量值 表示残差平方和 表示均值误差 表示预测准确度 预报正确次数 预测次数 四、模型假设 1L：假设题目中全部数据真实可靠，忽略误差；2L：假设观测站所在位置的经纬度准确无误；3L：假设天气预报针对的位置在所给网格点附近；4L：假设雨量在各网点之间的变动是连续的；五、问题分析 针对问题 1，我们将两种预测方法的所有预测值构造成两个以有序时间段对应的预测值为列,以网格点的个数为行的 2491164矩阵，对于 91 观测站点 41天的实测值做同样的处理，构造成 91164的矩阵。这样，繁琐的数据经过预处理后就整理成了三个矩阵。由于观测站点相应位置没有两种预测方法对应的预测值，无法直接进行评价，我们采用了三次样条插值的方法进行插值预处理，到了 91个观测站点两种预测方法的相应时刻的预测值，然后将两种预测方法雨量预测值与雨量实测值进行比较，从而判断出两种预测方法的准确性。针对问题 2，我们根据要求的雨量分级方式来考虑观众的感受。我们将问题 1 中 91个观测站点预测处理后雨量预测值构成的两个 91164矩阵和实际雨 量观测值构成的 91164 这个三个矩阵分别采用雨量等级记法构造出三个新的矩阵，然后分别把两个预测值构成的降雨量等级矩阵和观测值构成的等级矩阵对应元素相减并取绝对值，并进行等级统计，再利用预测准确度检验法进行判断，准确度越高说明我们预报的误差越小，表明预测方法更准确。六、模型建立 1、数据预处理（1）针对问题 1 根据上面的分析，我们先对数据进行预处理。处理方法为：把 FORECAST文件夹中的第一种和第二种预测方式得到的数据分开两个文件夹，分别以记事本格式按照日期的先后顺序有序的导入 Matlab 的 workspace 工作空间中，然后建立 m 文件编辑公式将两部分的数据导入 5347164的三维矩阵 rain1 和 rain2中；把 MEASURING 文件夹中的数据以同样方法导入 91741 的三维矩阵temp中，然后建立循环将 temp矩阵中每一层的后 4 列提取，另存入一个91164的 rain3 矩阵；在命令窗中直接导入预测点的经度和纬度存入矩阵 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 日第一个时间段的预测值构造成的矩阵 预测值 B1 B2 B3.B47 A1 .A2 .A3 .A52 .A53 .表 5-3 91 个站点在 41 天共 164 个时段的雨量实测值 预测值 D1 D2 D3.D164 P1 0 0 0.0 P2 0 0 0.0 P3 0 0 0.0.P90 0 0 .0 P91 0 0 .0（2）针对问题 2 由问题 1 我们可以整理得到 91 个观测点 41 天的预测值和测量值对应的两个 91164 矩阵，再根据问题 2 中气象部门将降雨的等级分为 6 个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为 0，得到两个预报等级矩阵，如下表 5-4 和 5-5：表 5-4 预测方法一在 6 月 18 日第一个时间段的预报等级构造成的矩阵 预报等级 B1 B2 B3.B47 A1 0 0 0.0 A2 0 0 0.0 A3 0 0 0.0.A52 0 0 0.0 A53 0 0 0.0 表 5-5 预测方法二在 6 月 18 日第一个时间段的预报等级构造成的矩阵 预报等级 B1 B2 B3.B47 A1 0 0 0.0 A2 0 0 0.0 A3 0 0 0.0.A52 0 0 0.0 A53 0 0 0.0 2、模型建立与求解（1）针对问题 1 由于 91 个观测站点没有相应的预测值，因此不能够直接对实测值进行评价，属于离散的散乱节点，我们利用插值基点为散乱节点的插值函数griddata1在Matlab 中进行三次样条插值处理，插值函数griddata为：(,1,1,)weapgriddata lat lon wea lat loncubic 1 其中lon表示预测点的经度值，lat表示预测点的纬度值，wea表示预测点的已知的预测值，1lat表示观察站点的纬度值，1lon表示观察站点的经度值，cubic表示三次样条插值的参数选项，weap观察站点的预测值。将 91 个观测站点 41 天 164 个时段的雨量情况进行预测之后，我们可以建立模型来评价这两种预测方法。这里我们利用残差平方和2以及平均误差avg3来作为评价的标准：21(),nijiiweapwear 1,2j 2 11,nijiiavgweapwearn 1,2j 3 最后，我们根据和avg的值进行评价，取值越大，表明预报的值准确性越低。因此，残差平方和与平均误差avg值较小的一种预测方法作为较好的预报方法。利用公式（1），我们在 Matlab 中应用编程求解，程序代码见附录。求解之后得到 91 个观测点 41 天 164 个时段的预测值，整理成 91164 矩阵，然后把预测矩阵和实测矩阵对应元素值相减取平方作残差平方和，再作平均误差，最后结果如下表 5-6：表 5-6 雨量预测计算结果数据表 结果 残差平方和 平均误差avg 预测方法一 预测方法二 由此可知，雨量预报方法一的准确性更高一些。（2）针对问题 2 由问题 1 我们可以整理得到 91 个观测点 41 天的预测值和测量值对应的两个 91164 矩阵，再根据问题 2 中气象部门将降雨的等级分为 6 个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为0，得到两个预报等级矩阵，继续利用残差平方和以及平均误差来作为评价的标准。残差平方和以及平均误差数值越小，表明预报越准确度越高，相应公众感受就越好。表 5-7 雨量等级计算结果数据表 结果 残差平方和 平均误差avg 预测方法一 2774 预测方法二 2806 由结果可知，预测方法一的准确度更高。七、模型优化 针对问题 2，由于残差平方和与平均误差难以反映真实汇报的准确度，我们将模型改进优化。由问题 1我们可以整理得到 91个观测点 41天的预测值和 测量值对应的两个 91164矩阵，再根据问题 2中气象部门将降雨的等级分为 6个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，得到两个预报等级矩阵，将两个预报等级矩阵与实测等级矩阵做差值运算，得到两个等级差矩阵，对等级差作绝对值处理后，我们就可以从中进行等级差统计。我们利用预测准确度检验法对两种预报进行评价。预测准确度（H）等于预报正确次数（R）（即运算之差为 0 的情况）和预测次数（T）之比，即 100%RHT 4 准确度越高，表明预报越准确度越高，相应公众感受就越好。雨量分级与统计程序见附录 4，数据统计结果见下表 7-1：表 7-1 两种预报等级与实测等级的级差与比例表 雨量等级差 预报 1 等级差个数 预报 1 等级差个数比例 预报 2 等级差个数 预报 2 等级差个数比例 0 12425 12404 1 2122 2128 2 55 52 3 10 10 4 1 2 5 0 0 0 0 6 0 0 0 0 由表 7-1 可见预报 1 的预报准确度为%高于预报 2的准确度%。预报 1与实测雨量的符合度高，那么两种预报方法中，公众更易接受第一种预报方法。八、模型的评价与推广 优点：本文思路自然流畅，较多地使用了原始数据，对91个站点的预测值比较准确。缺点：认为降雨量在不同的地区之间的变化有连续性，过于理想化。推广：本文对降雨量的预测方法也可以运用在工农业生产中其他方面，比如一段时期某地区的气温变化、全国粮食产量预测等。本文应用的散乱节点插值和残差分析法可以应用在很多方面。九、附录 附录 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 bytes isdir File_MEASURING=dir(MEASURING)File_MEASURING=43x1 struct array with fields:name date bytes isdir【】rain1=zeros(53,47,164);for n=1:164 filename=File_FORECAST1(n+2).name;rain1(:,:,n)=importdata(filename);end【】rain2=zeros(53,47,164);for n=1:164 filename=File_FORECAST2(n+2).name;rain2(:,:,n)=importdata(filename);end【】temp=zeros(91,7,41);for n=1:41 filename=File_MEASURING(n+2).name;temp(:,:,n)=importdata(filename);end rain3=zeros(91,164);i=1;for j=1:41 for k=4:7 rain3(:,i)=temp(:,k,j);i=i+1;end end 附录 2：问题 1 模型的建立与求解 calc1 pjwc=ccpfh=+005 calc2 pjwc=ccpfh=+005【】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);end pjwc=pjwc/(91*164);pjwc ccpfh【】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);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;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);end for i=1:164 for j=1:91 if weap1(j,i)rank1(j,i)=0;elseif weap1(j,i)=rank1(j,i)=1;elseif weap1(j,i)=rank1(j,i)=2;elseif weap1(j,i)=rank1(j,i)=3;elseif weap1(j,i)=rank1(j,i)=4;elseif weap1(j,i)=rank1(j,i)=5;else rank1(j,i)=6;end end end for i=1:164 for j=1:91 if weap2(j,i)rank2(j,i)=0;elseif weap2(j,i)=rank2(j,i)=1;elseif weap2(j,i)=rank2(j,i)=2;elseif weap2(j,i)=rank2(j,i)=3;elseif weap2(j,i)=rank2(j,i)=4;elseif weap2(j,i)=rank2(j,i)=5;else rank2(j,i)=6;end end end for i=1:164 for j=1:91 if rain3(j,i)rank(j,i)=0;elseif rain3(j,i)=rank(j,i)=1;elseif rain3(j,i)=rank(j,i)=2;elseif rain3(j,i)=rank(j,i)=3;elseif rain3(j,i)=rank(j,i)=4;elseif rain3(j,i)calc4 count=12425 12404 2428 2451 60 56 10 11 1 2 0 0 0 0 rate=0 0 0 0【】w1=abs(rank1-rank);w2=abs(rank2-rank);count=zeros(7,2);rate=zeros(7,2);for i=1:164 for j=1:91 if w1(j,i)=0 count(1,1)=count(1,1)+1;elseif w1(j,i)=1 count(2,1)=count(2,1)+1;elseif w1(j,i)=2 count(3,1)=count(3,1)+1;elseif w1(j,i)=3 count(4,1)=count(4,1)+1;elseif w1(j,i)=4 count(5,1)=count(5,1)+1;elseif w1(j,i)=5 count(6,1)=count(6,1)+1;elseif w1(j,i)=6 count(7,1)=count(7,1)+1;end end end for i=1:164 for j=1:91 if w2(j,i)=0 count(1,2)=count(1,2)+1;elseif w2(j,i)=1 count(2,2)=count(2,2)+1;elseif w2(j,i)=2 count(3,2)=count(3,2)+1;elseif w2(j,i)=3 count(4,2)=count(4,2)+1;elseif w2(j,i)=4 count(5,2)=count(5,2)+1;elseif w2(j,i)=5 count(6,2)=count(6,2)+1;elseif w2(j,i)=6 count(7,2)=count(7,2)+1;end end end rate=count./(91*164);count rate