1、雨量预报分析的评价模型一、摘要我们将 FORECAST 文件夹中的数据按日期先后顺序导入 Matlab，建立 5347164 的三维矩阵 rain1 和 rain2；把 MEASURING 文件夹中的数据以同样方法导入 91741 的三维矩阵 temp 中，然后建立循环将 temp 矩阵中每一层的后 4 列提取，另存入一个 91164 的rain3 矩阵；在命令窗中直接导入预测点的经度和纬度存入矩阵 lon 和 lat 中，导入实测点的经度和纬度存入矩阵 lon1 和 lat1 中，并对其作图，得到实测点和预测点的经纬度图。整理得到 91 个观测点 41 天的预测值和测量值对应的两个 9116
2、4 矩阵，根据气象部门将降雨的等级分为 6 个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为 0，得到两个预报等级矩阵。针对问题（1） ，利用插值基点为散乱节点的插值函数 1在 Matlab 中进行三gridat次样条插值处理，将 91 个观测站点 41 天 164 个时段的雨量情况进行预测。利用残差平方和 以及平均误差 来作为评价的标准。21()nijiiweapr1nijiiavwep残差平方和 与平均误差 值较小的一种预测方法作为较好的预报方法。残差平方和以vg及平均误差数值越小，表明预报越准确度越高。预测方法一的残差平方和为 174290.00，平均误
3、差为 0.4553。预测方法二的残差平方和为 195580.00，平均误差为 0.4753。雨量预报方法一的准确性更高一些。针对问题（2） ，两个预报等级矩阵，继续利用残差平方和以及平均误差来作为评价的标准。残差平方和以及平均误差数值越小，表明预报越准确度越高，相应公众感受就越好。预测方法一的残差平方和为 2774，平均误差为 0.1730。预测方法二的残差平方和为 2806，平均误差为 0.1745。雨量预报方法一的准确性更高一些。由于残差平方和与平均误差难以反映真实汇报的准确度，我们将模型改进优化。把矩阵中相应的降雨量值转化为其所对应等级值，得到两个预报等级矩阵，将两个预报等级矩阵与实测等
4、级矩阵做差值运算，得到两个等级差矩阵，对等级差作绝对值处理，进行等级差统计。我们利用预测准确度检验法对两种预报进行评价。预测准确度（ ）等于预报正H确次数（ ） （即运算之差为 0 的情况）和预测次数（ ）之比，即 。准RT10%RT确度越高，表明预报准确度越高，相应公众感受就越好。预报 1 的预报准确度为 83.26%高于预报 2 的准确度 83.11%，公众更易接受第一种预报方法。关键字：散乱节点插值 残差平方和 平均误差 预报等级矩阵 预测准确度二、问题重述雨量预报对农业生产和城市工作和生活有重要作用，但准确、及时地对雨量作出预报是一个十分困难的问题，广受世界各国关注。我国某地气象台和气
5、象研究所正在研究 6 小时雨量预报方法，即每天晚上 20 点预报从 21 点开始的 4 个时段（21 点至次日 3 点，次日3 点至 9 点，9 点至 15 点，15 点至 21 点）在某些位置的雨量，这些位置位于东经 120 度、北纬 32 度附近的 5347 的等距网格点上。同时设立 91 个观测站点实测这些时段的实际雨量，由于各种条件的限制，站点的设置是不均匀的。气象部门希望建立一种科学评价预报方法好坏的数学模型与方法。气象部门提供了 41 天的用两种不同方法的预报数据和相应的实测数据。预报数据在文件夹 FORECAST 中，实测数据在文件夹 MEASURING 中，其中的文件都可以用
6、Windows 系统的 “写字板 ”程序打开阅读。FORECAST 中的文件 lon.dat 和 lat.dat 分别包含网格点的经纬度，其余文件名为 _dis1 和_dis2，例如 f6181_dis1 中包含 2002 年 6 月 18 日晚上 20 点采用第一种方法预报的第一时段数据（其 2491 个数据为该时段各网格点的雨量） ，而 f6183_dis2 中包含 2002 年 6 月 18 日晚上 20 点采用第二种方法预报的第三时段数据。MEASURING 中包含了 41 个名为.SIX 的文件，如 020618.SIX 表示 2002 年 6 月 18日晚上 21 点开始的连续 4
7、 个时段各站点的实测数据（雨量） ，这些文件的数据格式是：站号 纬度 经度 第 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.0000 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
8、1.5000 1.9000雨量用毫米做单位，小于 0.1 毫米视为无雨。(1) 请建立数学模型来评价两种 6 小时雨量预报方法的准确性；(2) 气象部门将 6 小时降雨量分为 6 等：0.1 2.5 毫米为小雨，2.66 毫米为中雨，6.112毫米为大雨，12.1 25 毫米为暴雨，25.160 毫米为大暴雨，大于 60.1 毫米为特大暴雨。若按此分级向公众预报，如何在评价方法中考虑公众的感受？三、名词和符号说明 lon表示预测点的经度值at表示预测点的纬度值we表示预测点已知的预测值1l表示观察站点的经度值t表示观察站点的纬度值cubi表示三次样条插值的参数选项eap观察站点的预测值ijw表
9、示两种预测值矩阵中第 i 个元素的测量值r表示实测矩阵中第 i 个元素的测量值表示残差平方和avg表示均值误差H表示预测准确度R预报正确次数T预测次数四、模型假设 1L：假设题目中全部数据真实可靠，忽略误差；2：假设观测站所在位置的经纬度准确无误；3：假设天气预报针对的位置在所给网格点附近；4：假设雨量在各网点之间的变动是连续的；五、问题分析针对问题 1，我们将两种预测方法的所有预测值构造成两个以有序时间段对应的预测值为列, 以网格点的个数为行的 2491164 矩阵，对于 91 观测站点 41 天的实测值做同样的处理，构造成 91164 的矩阵。这样，繁琐的数据经过预处理后就整理成了三个矩阵
10、。由于观测站点相应位置没有两种预测方法对应的预测值，无法直接进行评价，我们采用了三次样条插值的方法进行插值预处理，到了 91 个观测站点两种预测方法的相应时刻的预测值，然后将两种预测方法雨量预测值与雨量实测值进行比较，从而判断出两种预测方法的准确性。针对问题 2，我们根据要求的雨量分级方式来考虑观众的感受。我们将问题 1 中 91 个观测站点预测处理后雨量预测值构成的两个 91164 矩阵和实际雨量观测值构成的 91164这个三个矩阵分别采用雨量等级记法构造出三个新的矩阵，然后分别把两个预测值构成的降雨量等级矩阵和观测值构成的等级矩阵对应元素相减并取绝对值，并进行等级统计，再利用预测准确度检验
11、法进行判断，准确度越高说明我们预报的误差越小，表明预测方法更准确。六、模型建立1、数据预处理（1 ）针对问题 1根据上面的分析，我们先对数据进行预处理。处理方法为：把 FORECAST 文件夹中的第一种和第二种预测方式得到的数据分开两个文件夹，分别以记事本格式按照日期的先后顺序有序的导入 Matlab 的 workspace 工作空间中，然后建立 m 文件编辑公式将两部分的数据导入 5347164 的三维矩阵 rain1 和 rain2 中；把 MEASURING 文件夹中的数据以同样方法导入 91741 的三维矩阵 temp 中，然后建立循环将 temp 矩阵中每一层的后 4 列提取，另存入
12、一个 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 日第一个时间段的预测值构造成的矩阵预测值 B1 B2 B3 . B47A1 0.03
13、10 0.0321 0.0313 . 0.0081A2 0.0302 0.0310 0.0315 . 0.0078A3 0.0340 0.0341 0.0337 . 0.0079. . . . . .A52 0.0088 0.0089 0.0084 . 0.0061A53 0.0089 0.0086 0.0086 . 0.0062表 5-2 预测方法二在 6 月 18 日第一个时间段的预测值构造成的矩阵预测值 B1 B2 B3 . B47A1 0.0319 0.0293 0.0322 . 0.0077A2 0.0288 0.0311 0.0303 . 0.0085A3 0.0325 0.032
14、3 0.0313 . 0.0082. . . . . .A52 0.0095 0.0089 0.0079 . 0.0056A53 0.0085 0.0094 0.0093 . 0.0057表 5-3 91 个站点在 41 天共 164 个时段的雨量实测值预测值 D1 D2 D3 . D164P1 0 0 0 . 0P2 0 0 0 . 0P3 0 0 0 . 0. . . . . .P90 0 0 4.1 . 0P91 0 0 0.2 . 0（2）针对问题 2由问题 1 我们可以整理得到 91 个观测点 41 天的预测值和测量值对应的两个 91164矩阵，再根据问题 2 中气象部门将降雨的等级
15、分为 6 个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为 0，得到两个预报等级矩阵，如下表 5-4 和 5-5：表 5-4 预测方法一在 6 月 18 日第一个时间段的预报等级构造成的矩阵预报等级 B1 B2 B3 . B47A1 0 0 0 . 0A2 0 0 0 . 0A3 0 0 0 . 0. . . . . .A52 0 0 0 . 0A53 0 0 0 . 0表 5-5 预测方法二在 6 月 18 日第一个时间段的预报等级构造成的矩阵预报等级 B1 B2 B3 . B47A1 0 0 0 . 0A2 0 0 0 . 0A3 0 0 0 . 0. .
16、 . . . .A52 0 0 0 . 0A53 0 0 0 . 02、模型建立与求解（1）针对问题 1由于 91 个观测站点没有相应的预测值，因此不能够直接对实测值进行评价，属于离散的散乱节点，我们利用插值基点为散乱节点的插值函数 1在 Matlab 中进行三次样条gridat插值处理，插值函数 为：gridat1(,1,)wepltonweltoncubi其中 表示预测点的经度值， 表示预测点的纬度值， 表示预测点的已知的lon e预测值， 表示观察站点的纬度值， 表示观察站点的经度值， 表示三次样条1at插值的参数选项， 观察站点的预测值。将 91 个观测站点 41 天 164 个时段的
17、雨量情况进行预测之后，我们可以建立模型来评价这两种预测方法。这里我们利用残差平方和2以及平均误差 3来作为评价的标准：vg221(),nijiiweapr1,j3,ijiiv,最后，我们根据 和 的值进行评价，取值越大，表明预报的值准确性越低。因此，残ag差平方和 与平均误差 值较小的一种预测方法作为较好的预报方法。利用公式（1） ，我们在 Matlab 中应用编程求解，程序代码见附录 2.2。求解之后得到 91 个观测点 41 天 164 个时段的预测值，整理成 91164 矩阵，然后把预测矩阵和实测矩阵对应元素值相减取平方作残差平方和，再作平均误差，最后结果如下表 5-6：表 5-6 雨量
18、预测计算结果数据表结果 残差平方和 平均误差 avg预测方法一 174290.00 0.4553预测方法二 195580.00 0.4753由此可知，雨量预报方法一的准确性更高一些。（2 ）针对问题 2由问题 1 我们可以整理得到 91 个观测点 41 天的预测值和测量值对应的两个 91164矩阵，再根据问题 2 中气象部门将降雨的等级分为 6 个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，其中，预测中的零全部记为 0，得到两个预报等级矩阵，继续利用残差平方和以及平均误差来作为评价的标准。残差平方和以及平均误差数值越小，表明预报越准确度越高，相应公众感受就越好。表 5-7 雨量等级
19、计算结果数据表结果 残差平方和 平均误差 avg预测方法一 2774 0.1730预测方法二 2806 0.1745由结果可知，预测方法一的准确度更高。七、模型优化针对问题 2，由于残差平方和与平均误差难以反映真实汇报的准确度，我们将模型改进优化。由问题 1 我们可以整理得到 91 个观测点 41 天的预测值和测量值对应的两个91164 矩阵，再根据问题 2 中气象部门将降雨的等级分为 6 个等级的分法，把矩阵中相应的降雨量值转化为其所对应等级值，得到两个预报等级矩阵，将两个预报等级矩阵与实测等级矩阵做差值运算，得到两个等级差矩阵，对等级差作绝对值处理后，我们就可以从中进行等级差统计。我们利用
20、预测准确度检验法对两种预报进行评价。预测准确度（ ）等于预报正确H次数（ ） （即运算之差为 0 的情况）和预测次数（ ）之比，即RT10%RH4准确度越高，表明预报越准确度越高，相应公众感受就越好。雨量分级与统计程序见附录 4，数据统计结果见下表 7-1：表 7-1 两种预报等级与实测等级的级差与比例表雨量等级差 预报 1等级差个数预报 1等级差个数比例预报 2等级差个数预报 2等级差个数比例0 12425 0.8326 12404 0.83111 2122 0.1422 2128 0.14262 55 0.0037 52 0.00353 10 0.0007 10 0.00074 1 0.0
21、001 2 0.00015 0 0 0 06 0 0 0 0由表 7-1 可见预报 1 的预报准确度为 83.26%高于预报 2 的准确度 83.11%。预报 1 与实测雨量的符合度高，那么两种预报方法中，公众更易接受第一种预报方法。八、模型的评价与推广优点：本文思路自然流畅，较多地使用了原始数据，对 91 个站点的预测值比较准确。缺点：认为降雨量在不同的地区之间的变化有连续性，过于理想化。推广：本文对降雨量的预测方法也可以运用在工农业生产中其他方面，比如一段时期某地区的气温变化、全国粮食产量预测等。本文应用的散乱节点插值和残差分析法可以应用在很多方面。九、附录附录 1：数据初始化 File_
22、FORECAST1=dir(预测方法一 )File_FORECAST1 =166x1 struct array with fields:namedatebytesisdir File_FORECAST2=dir(预测方法二 )File_FORECAST2 = 166x1 struct array with fields:namedatebytesisdir File_MEASURING=dir(MEASURING)File_MEASURING = 43x1 struct array with fields:namedatebytesisdir【mypro1.m】rain1=zeros(53,4
23、7,164);for n=1:164filename=File_FORECAST1(n+2).name;rain1(:,:,n)=importdata(filename);end【mypro2.m】rain2=zeros(53,47,164);for n=1:164filename=File_FORECAST2(n+2).name;rain2(:,:,n)=importdata(filename);end【mypro3.m】temp=zeros(91,7,41);for n=1:41filename=File_MEASURING(n+2).name;temp(:,:,n)=importdata
24、(filename);endrain3=zeros(91,164);i=1;for j=1:41for k=4:7rain3(:,i)=temp(:,k,j);i=i+1;endend附录 2：问题 1 模型的建立与求解 calc1pjwc =0.4553ccpfh =1.7429e+005 calc2pjwc =0.4753ccpfh =1.9558e+005【calc1.m】lat;lon;lat1;lon1;ccpfh=0;pjwc=0;for i=1:164wea=rain1(:,:,i);wear=rain3(:,i);weap=griddata(lat,lon,wea,lat1,l
25、on1);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】lat;lon;lat1;lon1;ccpfh=0;pjwc=0;for i=1:164wea=rain2(:,:,i);wear=rain3(:,i);weap=griddata(lat,lon,wea,lat1,lon1);chazhi=wear-weap;pingfang=chazhi.2;ju
26、edui=abs(chazhi);ccpfh=ccpfh+sum(pingfang,1);pjwc=pjwc+sum(juedui);endpjwc=pjwc/(91*164);pjwcccpfh附录 3：问题 2 模型的建立与求解 calc3pjwc1 =0.1730ccpfh1 =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:164wea1=rain1(:,:,
27、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);endfor i=1:164for j=1:91if weap1(j,i) calc4count =12425 124042428 245160 5610 111 20 00 0rate =0.8326 0.83110.1627 0.16420.0040 0.00380.0007 0.00070.0001 0.00010 0
28、0 0【calc4.m】w1=abs(rank1-rank);w2=abs(rank2-rank);count=zeros(7,2);rate=zeros(7,2);for i=1:164for j=1:91if w1(j,i)=0count(1,1)=count(1,1)+1;elseif w1(j,i)=1count(2,1)=count(2,1)+1;elseif w1(j,i)=2count(3,1)=count(3,1)+1;elseif w1(j,i)=3count(4,1)=count(4,1)+1;elseif w1(j,i)=4count(5,1)=count(5,1)+1;
29、elseif w1(j,i)=5count(6,1)=count(6,1)+1;elseif w1(j,i)=6count(7,1)=count(7,1)+1;endendendfor i=1:164for j=1:91if w2(j,i)=0count(1,2)=count(1,2)+1;elseif w2(j,i)=1count(2,2)=count(2,2)+1;elseif w2(j,i)=2count(3,2)=count(3,2)+1;elseif w2(j,i)=3count(4,2)=count(4,2)+1;elseif w2(j,i)=4count(5,2)=count(5,2)+1;elseif w2(j,i)=5count(6,2)=count(6,2)+1;elseif w2(j,i)=6count(7,2)=count(7,2)+1;endendendrate=count./(91*164);countrate
