雨量预报评价模型.docx
《雨量预报评价模型.docx》由会员分享，可在线阅读，更多相关《雨量预报评价模型.docx（9页珍藏版）》请在冰豆网上搜索。
雨量预报评价模型
雨量预报方法的评价模型
李泳，易勋，贺望香
1.问题分析
雨量预报对农业生产和城市工作和生活有重要作用，但准确、及时地对雨量作出预报是一个十分困难的问题。
针对问题一，气象部门提供了41天实测数据并且希望建立数学模型来评价两种雨量预报方法的准确性。
通过分析数据可知，网格点数据代表的是53*47的网格点的预测雨量，而实测数据是以经纬度定位的观测点的实测雨量。
很显然无法对两个定位方式不同的数据进行插值。
所以我们选取与对应观测站点最近的5个网格点，通过计算网格点对对应观测站点的权重和网格点的预测雨量的乘积的和求出观测站点的预测雨量。
然后就可以通过对两种方法的预测雨量与实测雨量进行比较得出哪种方法更准确。
由于两种方法与实测雨量很相近，很难对比出两种方法的准确性，所以我们采用预报偏差率来比较两种6小时雨量预报方法的准确性。
2.模型假设
（1）观测站点之间距离的设置是不同的。
（2）雨量用毫米做单位，小于0.1毫米视为无雨
3.符号说明
：
表示第n个网格点到第i个观测站点的距离；（其中n=1,2,3,4,5;i=1-91）
：
表示与第i个观测站点的距离最小的前5个网格点的权重；
：
表示第i个网格点分别在某月某日某个时段的雨量值；（其中n=1,2,3,4,5）
：
表示第i个观测点第j个时段的预测雨量值；
：
表示第i个观测点第j个时段的实测雨量值。
4.模型的建立与求解
4.1.1首先，找出网格点与观测站点的散点图.从图中可以看出，观测站（红点）与网格点（蓝点）分布不均匀。
且观测站点分布在网格点中间部分。
4.1.2筛选出每一个站点周围的5个距离观测站点最近的网格点。
5个网格点的选取通过先给定最大与最小值的范围，然后利用matlab的find（）函数找出符合筛选条件的全部点。
（见附录2）
4.1.3运用欧拉公式求出符合条件的5个网格点到观测站点的距离的集合，并对其进行从小到大排序。
欧拉距离公式：
4.1.4运用欧拉的倒数加权法计算出5个网格点到观测站的权重[1]
计算公式：
4.1.5计算观测点的雨量预测值，取得每一个5个网格某月某日某时段对应的预测雨量值乘以其对应的权重即等于各个观测点各个时段的预测雨量值
。
计算公式：
4.1.6已经得到了两种方法的观测点的预测雨量，接下来则与实测数据比较，分析哪种方法的准确性高。
①通过比较方法1与方法2的4个时段的预测雨量的平均值（见附录四）
预测雨量的平均值
第一时段
第二时段
第三时段
第四时段
方法一
1.2112
1.3588
3.9047
2.5623
方法二
2.4127
3.1605
2.7383
2.6722
实测雨量的平均值
第一时段
第二时段
第三时段
第四时段
实测
1.4827
2.1179
1.8007
1.6362
②通过对上面两个表格的分析可得平均值无法得到准确评估，故采用求预报偏差率的方法来判断两种方法的准确性。
计算公式：
③由于最终结果依然是一个164*91矩阵，所以对其求和之后再对他进行求列平均得到一个行向量，然后再对行向量平均得到一个均值
通过matlab程序，最终得到：
第一种方法预报偏离的平均值为124.3207
第二种方法预报偏离的平均值为124.4289
5.参考文献
[1]《雨量预报数学模型》论文
[2]欧拉,《基于欧式距离的最近邻改进算法》，
[3]H.P.Williams,《数学规划模型建立与计算机应用》，北京，国防工业出版社，1991年；
六附录
附录1
x=lat;
y=lon;
xi=A020618（:
2）;%%为了不让编译器报错，故将文件名第一个改为A开头
yi=A020618（:
3）;
plot（x,y,'.b',xi,yi,'.r'）
axis（[27.335.3116.8125]）;
xlabel（'经度'）,ylabel（'纬度'）,title（'观测站点与网格点分布图'）
附录2
x=lat;
y=lon;
xi=A020618（:
2）;%为了不让编译报错，故将文件改名为A字母开头
yi=A020618（:
3）;
max_x1=max（xi）;
min_x2=min（xi）;
max_y1=max（yi）;
min_y2=min（yi）;
b1=find（x>=min_x2&x<=max_x1）;%满足条件的点的x坐标位置
sWang_x=x（b1）;%通过筛选后的序号找到对应网格点的坐标
sWang_y=y（b1）;
b3=find（sWang_y>=min_y2&sWang_y<=max_y1）;%筛选后满足条件的点的y坐标位置
youXiao_x=sWang_x（b3）;%有效的点
youXiao_y=sWang_y（b3）;
plot（youXiao_x,youXiao_y,'.b',xi,yi,'.r'）;
xlabel（'经度'）,ylabel（'纬度'）,title（'筛选后有效点的分布图'）
附录3
x=lat;
y=lon;
xi=A020618（:
2）;%为了不让编译报错，故将文件改名为A字母开头
yi=A020618（:
3）;
max_x1=max（xi）;
min_x2=min（xi）;
max_y1=max（yi）;
min_y2=min（yi）;
b1=find（x>=min_x2&x<=max_x1）;%满足条件的点的x坐标位置
sWang_x=x（b1）;%通过筛选后的序号找到对应网格点的坐标
sWang_y=y（b1）;
b3=find（sWang_y>=min_y2&sWang_y<=max_y1）;%筛选后满足条件的点的y坐标位置
youXiao_x=sWang_x（b3）;%有效的点
youXiao_y=sWang_y（b3）;
xb4=xi;
yb4=yi;
l=length（youXiao_x）;
d=zeros（91,l）;
fori=1:
91
forj=1:
l
d（i,j）=sqrt（（xb4（i）-youXiao_x（j））^2+（yb4（i）-sWang_y（j））^2）;%用欧拉公式计算距离
end
end
附录4（实测平均值）
sum1=0;
sum2=0;
sum3=0;
sum4=0;
fori=1:
4:
164
forj=1:
91
sum1=sum1+B1（i,j）;
end
end
avg1=sum1/（41*91）;
fori=2:
4:
164
forj=1:
91
sum2=sum2+B1（i,j）;
end
end
avg2=sum2/（41*91）;
fori=3:
4:
164
forj=1:
91
sum3=sum3+B1（i,j）;
end
end
avg3=sum3/（41*91）;
fori=4:
4:
164
forj=1:
91
sum4=sum4+B1（i,j）;
end
end
avg4=sum4/（41*91）;
求预测平均值
fori=1:
4:
164
sum1=sum1+A1（i）;
end
avg
（1）=sum1/41;
fori=2:
4:
164
sum2=sum2+A1（i）;
end
avg
（2）=sum2/41;
fori=3:
4:
164
sum3=sum3+A1（i）;
end
avg（3）=sum3/41;
fori=4:
4:
164
sum4=sum4+A1（i）;
end
avg（4）=sum4/41;
fori=1:
4:
164
sum1=sum1+A2（i）;
end
avg（2,1）=sum1/41;
fori=2:
4:
164
sum2=sum2+A2（i）;
end
avg（2,2）=sum2/41;
fori=3:
4:
164
sum3=sum3+A2（i）;
end
avg（2,3）=sum3/41;
fori=4:
4:
164
sum4=sum4+A2（i）;
end
avg（2,4）=sum4/41;
end
附录5（参考：
《雨量预报数学模型论文》）
计算偏差率
w=reshape（a',455,1）';%将91*5个点转换成455*1
A11=A1（:
w）;%第一类164*455个降雨量
A22=A2（:
w）;%第二类164*455个降雨量
gg=repmat（HH,164,1）%为使HH与A11，A22大小相同，扩展平铺HH，gg大小为164*455
A1_1=gg.*A11%第一类方法权与降雨量乘积
A2_2=gg.*A22%第二类方法权与降雨量乘积大小为
%%%%%%%%分别对A1_1，A2_2每隔5列求和并返回sum1,sum2%%%%%%%%%
sum1=[];
fori=1:
5:
454
sum1=[sum1A1_1（:
i）+A1_1（:
i+1）+A1_1（:
i+2）+A1_1（:
i+3）+A1_1（:
i+4）];
end
sum2=[];
fori=1:
5:
454
sum2=[sum2A2_2（:
i）+A2_2（:
i+1）+A2_2（:
i+2）+A2_2（:
i+3）+A2_2（:
i+4）];
end
B1=[];
formonth1=6:
7
ifmonth1==6;
fordata1=18:
28
%str3=['020'int2str（month1）int2str（data1）];
str3=['020'int2str（month1）int2str（data1）'.SIX'];
tmp3=load（str3）;
B1=[B1;tmp3（:
4:
7）'];
end
else
fordata1=1:
30
ifdata1<10strtmp=['0'int2str（data1）];
else
strtmp=int2str（data1）;
end
str3=['020'int2str（month1）strtmp'.SIX'];
tmp3=load（str3）;
B1=[B1;tmp3（:
4:
7）'];
end
end
end
附录6
B_1=abs（sum1-B1）./sum1;%第一种方法的偏差率
B_2=abs（sum2-B1）./sum2;%第二种方法的偏差率
B_1（find（isnan（B_1）==1））=0;%将B_1中NaN替换成为0
B_2（find（isnan（B_2）==1））=0;%将B_2中NaN替换成为0
sum_1=sum（B_1）;
sum_2=sum（B_2）;
D1=mean（sum_1）%列平均，得到行向量
D_1=mean（D1）%再对行向量平均，得到一个值
D2=mean（sum_2）
D_2=mean（D2）
升级会员 