雨量預報論文

1、雨量預報方法的評價摘要本文研究了雨量預報的準確性問題,運用軟件編程，通過插值和對比計數(shù)的方法，比較方法一和方法二預測的準確性；問題二研究的是公眾對雨量預報的滿意程度，運用模糊數(shù)學的知識，將模糊的評價等級數(shù)字化，用偏小型柯西隸屬函數(shù)來刻畫公眾對兩種方法的滿意度。問題一 分別對每一天每一時段兩種方法預測的雨量進行插值，得到兩種方法在91個站點上的預測雨量，并分別與實際雨量對應(yīng)進行比較，哪種方法預測的誤差較小的總的次數(shù)最多，即為預測的比較準確的方法，由編程求解得知方法一預報的比方法二更準確。問題二 為了刻畫公眾對兩種預報方法的滿意度，我們引進了偏小型柯西隸屬函數(shù)，同樣地對每一天每一時段兩種方法預測的

2、雨量進行插值，得到兩種方法在91個站點上的預測雨量并按照給定的范圍等級化，求出兩種方法與實測數(shù)據(jù)的等級差，并代入確定的隸屬函數(shù)，最終將各天各時段各站點的滿意度加起來，即求得方法一、方法二的總滿意度、。由編程解得=14394 ，=14923 ，即方法二的滿意度高于方法一。關(guān)鍵詞:插值偏小型柯西隸屬函數(shù)對比計數(shù)1、 問題重述雨量預報對農(nóng)業(yè)生產(chǎn)和城市工作和生活有重要作用，但準確、及時地對雨量作出預報是一個十分困難的問題，廣受世界各國關(guān)注。我國某地氣象臺和氣象研究所正在研究6小時雨量預報方法，即每天晚上20點預報從21點開始的4個時段（21點至次日3點，次日3點至9點，9點至15點，15點至21點）在

3、某些位置的雨量，這些位置位于東經(jīng)120度、北緯32度附近的53×47（緯度）的等距網(wǎng)格點上。同時設(shè)立91個觀測站點實測這些時段的實際雨量，由于各種條件的限制，站點的設(shè)置是不均勻的。氣象部門希望建立一種科學評價預報方法好壞的數(shù)學模型與方法。氣象部門提供了41天的用兩種不同方法的預報數(shù)據(jù)和相應(yīng)的實測數(shù)據(jù)。預報數(shù)據(jù)在文件夾FORECAST中，實測數(shù)據(jù)在文件夾MEASURING中，其中的文件都可以用Windows系統(tǒng)的“寫字板”程序打開閱讀。FORECAST中的文件lon.dat和lat.dat分別包含網(wǎng)格點的經(jīng)緯度，其余文件名為<f日期i>_dis1和<f日期i>_

4、dis2，例如f6181_dis1中包含2002年6月18日晚上20點采用第一種方法預報的第一時段數(shù)據(jù)（其2491個數(shù)據(jù)為該時段各網(wǎng)格點的雨量），而f6183_dis2中包含2002年6月18日晚上20點采用第二種方法預報的第三時段數(shù)據(jù)。MEASURING中包含了41個名為<日期>.SIX的文件，如020618.SIX表示2002年6月18日晚上21點開始的連續(xù)4個時段各站點的實測數(shù)據(jù)（雨量），這些文件的數(shù)據(jù)格式是：站號 緯度 經(jīng)度 第1段 第2段 第3段 第4段 58138 32.9833 118.5167 0.0000 0.2000 10.1000 3.100058139 33

5、.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 1.5000 1.9000雨量用毫米做單位，小于0.1毫米視為無雨。(1) 請建立數(shù)學模型來評價兩種6小時雨量預報方法的準確性；(2) 氣象部門將6小時降雨量分為6等：0.12.5毫米為小雨，2.66毫米為中雨，6.112毫米為大雨，12.1

6、25毫米為暴雨，25.160毫米為大暴雨，大于60.1毫米為特大暴雨。若按此分級向公眾預報，如何在評價方法中考慮公眾感受.2、 問題分析針對問題一，在53×47的經(jīng)緯度等距網(wǎng)格點上有方法一、二的預測雨量，在91個觀測站點上有41天每天4個時段的實測雨量，而且站點的設(shè)置不與網(wǎng)格格點重合。把經(jīng)緯度分別設(shè)為x，y軸，雨量設(shè)為z軸，對方法一預測的數(shù)據(jù)進行插值，求出由方法一得到的91個站點的預測雨量，同理求出91個站點由方法二得到的預測雨量，再與91個站點的實測雨量進行作差得到某站點的誤差，誤差較小者則相應(yīng)地讓自增1,91個站點比完之后再對4個時段進行比較，中較大者相應(yīng)地讓自增1，然后再對41

7、天進行比較，中較大者相應(yīng)地讓自增1，最后可得到兩種方法在41天內(nèi)更準確的天數(shù)；針對問題二，循環(huán)分別讀取出某天某時段的91個站點的實測雨量，和網(wǎng)格點預測雨量，用散亂節(jié)點插值函數(shù)插值求得某天某時段91個站點的預測雨量，由此便得到某天某時段91個站點的實測雨量、方法一預測雨量、方法二預測雨量，接著讓、和按照題目給出的范圍進行雨量的等級劃分，然后對方法一、二分別求出等級誤差、和雨量誤差、；若等級誤差為0，則滿意度為1，否則，按照分類討論中的表格進行等級誤差的確定，并代入函數(shù)求出某站點的滿意度；最后把41天的4個時段的91個站點的滿意度分別加起來得到兩種方法的總滿意度和. 3、模型假設(shè)（1）實際測量的雨

8、量沒有誤差，不受風速等影響；（2）方法一、二預測的雨量不會因為人為因素而錯誤。（3）忽略地球的海拔高低，把經(jīng)緯度分別作為橫縱坐標建立平面直角坐標系；4、符號說明符號含義方法一預測的觀測站點的雨量方法二預測的觀測站點的雨量方法一預測的觀測站點的雨量與實際雨量的誤差方法二預測的觀測站點的雨量與實際雨量的誤差 某天某時段91個站點方法一誤差比較小的次數(shù) 某天某時段91個站點方法二誤差比較小的次數(shù)某天4個時段方法一更準確的次數(shù)某天4個時段方法二更準確的次數(shù)41天中方法一更準確的天數(shù) 41天中方法二更準確的天數(shù)方法一預測網(wǎng)格點的雨量方法二預測網(wǎng)格點的雨量41天j91個站點k4個時段41天實測數(shù)據(jù)的存儲地

9、址字符串某天的實測數(shù)據(jù)實測經(jīng)度實測緯度某天4個時段的實測站點降雨量網(wǎng)格點的經(jīng)度網(wǎng)格點的緯度某天某個時段方法一預測數(shù)據(jù)的字符串地址某天某個時段方法二預測數(shù)據(jù)的字符串地址某天某個時段方法一預測的網(wǎng)格點降雨量某天某個時段方法二預測的網(wǎng)格點降雨量方法一插值得到的91個站點預測雨量方法二插值得到的91個站點預測雨量方法一的等級差方法二的等級差方法一的雨量差方法二的雨量差方法一、二對應(yīng)的等級差方法一某天某時段某站點的滿意度方法二某天某時段某站點的滿意度方法一、二的總滿意度方法一某天某時段的滿意度方法二某天某時段的滿意度方法一某天的滿意度方法二某天的滿意度5、 分析與建立模型5.1問題一5.1.1模型建立（

10、一）目標確定比較方法一、方法二預測的準確性，總體來看哪種方法預測的誤差較小。（二）模型的建立分別以經(jīng)、緯度為x、y坐標建立平面直角坐標系，得到網(wǎng)格點與實測點的相對位置，如下圖：通過循環(huán)分別讀取出某天某時段的91個站點的實測雨量，和網(wǎng)格點預測雨量，用散亂節(jié)點插值函數(shù)插值求得某天某時段91個站點的預測雨量，由此便得到某天某時段91個站點的實測雨量、方法一預測雨量、方法二預測雨量，接著讓、分別和進行比較，91個站點循環(huán)完后，比較與，即可得到某天某時段中更準確的方法，再進行4個時段的循環(huán)，比較和即可得到某天內(nèi)更準確的方法，最后再對41天循環(huán)，比較和， 即可得出41天中更為準確的方法；具體實施步驟參考模

11、型求解。5.1.2 模型求解（一）求解:利用函數(shù)從6月18日到7月30日循環(huán)讀取實測站點數(shù)據(jù)，并賦給矩陣，從矩陣提取出第2、3、4到7列分別賦給矩陣、 、；：讀取網(wǎng)格點的經(jīng)、緯度賦給矩陣、；：讓時段從1到4循環(huán)，分別讀取方法一、二預測的網(wǎng)格點的雨量賦給矩陣、，并通過函數(shù)插值得到方法一、二在91個站點上預測雨量、；：站點從1到91循環(huán)，對應(yīng)站點的預測雨量和實測雨量相減取絕對值，分別得到方法一、二的誤差、；：與中小者則對應(yīng)讓或自增1，知直到91個站點循環(huán)完；然后再讓與比較，較大者則對應(yīng)地讓或自增1，直到4個時段循環(huán)完；接著讓與比較，兩者中較大者對應(yīng)地讓或自增1，直到41天循環(huán)完（二）求解結(jié)果與分析

12、最后求得=17，=10,；從結(jié)果可知，在41天中對兩種方法進行4個時段的比較，其中有17天方法一更準確，有10天方法二更準確，還有14天，方法一、二在4個時段中分別都有兩個時段更準確，即兩種方法在這14天的4個時段均一樣準確。綜上所述，方法一比方法二更為準確；5.2問題二5.2.1模型建立（一）目標確定通過引入偏小型柯西隸屬分布函數(shù)，求出某天某時段某個站點的滿意度，然后對41天的4個時段的91個站點的滿意度進行相加，求出兩種方法的總滿意度和；（二）模型準備1、引入偏小型柯西隸屬分布函數(shù)為了描述滿意度，我們引入了偏小型柯西隸屬分布函數(shù)，其中橫坐標表示降雨量等級差：當=1時，=0.8,；當=6時，

13、=0.01；令=0把以上兩點代入函數(shù)中科解得=0.25，=3.34，即柯西隸屬分布函數(shù)為函數(shù)圖像如下圖：2、分類討論 為了防止誤差太大，我們細化了幾種情況：方法一、二的站點預測降雨量 1 2 1 2 3 1 2 4 4 5 6（三）模型的建立通過循環(huán)分別讀取出某天某時段的91個站點的實測雨量，和網(wǎng)格點預測雨量，用散亂節(jié)點插值函數(shù)插值求得某天某時段91個站點的預測雨量，由此便得到某天某時段91個站點的實測雨量、方法一預測雨量、方法二預測雨量，接著讓、和按照題目給出的范圍進行雨量的等級劃分，然后對方法一、二分別求出等級誤差、和雨量誤差、；若等級誤差為0，則滿意度為1，否則，按照分類討論中的表格進行

14、等級誤差的確定，并代入函數(shù)求出某站點的滿意度；最后把41天的4個時段的91個站點的滿意度分別加起來得到兩種方法的總滿意度和5.2.2 模型求解（一）求解:利用函數(shù)從6月18日到7月30日循環(huán)讀取實測站點數(shù)據(jù)，并賦給矩陣，從矩陣提取出第2、3、4到7列分別賦給矩陣、 、,讓矩陣=，并讓按照題目給出的雨量范圍進行降雨量等級劃分；：讀取網(wǎng)格點的經(jīng)、緯度賦給矩陣、；：讓時段從1到4循環(huán)，分別讀取方法一、二預測的網(wǎng)格點的雨量賦給矩陣、，并通過函數(shù)插值得到方法一、二在91個站點上預測雨量、，分別讓=、=，并讓、按照題目給出的雨量范圍進行降雨量等級劃分；：站點從1到91循環(huán)，對應(yīng)站點的余量等級對應(yīng)相減取絕對

15、值，分別得到方法一、二的等級誤差、，對應(yīng)站點的預測雨量和實測雨量相減取絕對值，分別得到方法一、二的雨量誤差、；：對兩種方法，若等級差為0，則滿意度s=0；否則，則按分類討論中的情況編程即可，對應(yīng)的等級差代入隸屬函數(shù)；最終把41天的4個時段的91個站點的滿意度加起來得到方法一、二的總滿意度和 （二）求解結(jié)果與分析最終求得，方法一的總滿意度=14394，方法二的總滿意度=14923。由結(jié)果可知，方法二的滿意度略高于方法一。附錄問題一程序 ：clear allclcM1=0;M2=0;for i=618:628,701:730 fm='dataMEASURING','020&

16、#39;,num2str(i),'.SIX'%實測字符串41 tmp=dlmread(fm);%實測數(shù)據(jù)41 latm=tmp(:,2);%實測緯度 lonm=tmp(:,3);%實測經(jīng)度 rqm=tmp(:,4:7);%4個時段 lonf=dlmread('dataFORECASTlon.dat','');%預測網(wǎng)格點經(jīng)度 latf=dlmread('dataFORECASTlat.dat','');%預測網(wǎng)格點緯度 N1=0;N2=0; for k=1:4 ff1='dataFORECAST',

17、'f',num2str(i),num2str(k),'_dis1'%第1方法預測雨量41x4 ff2='dataFORECAST','f',num2str(i),num2str(k),'_dis2' rqf1=dlmread(ff1);%第1方法預測雨量53x47 rqf2=dlmread(ff2); rqm_f1=griddata(lonf,latf,rqf1,lonm,latm,'V4');%插值得91個站點的預測值 rqm_f2=griddata(lonf,latf,rqf2,lonm,lat

18、m,'V4'); n1=0;n2=0; for j=1:91 a=abs(rqm(j,k)-rqm_f1(j); b=abs(rqm(j,k)-rqm_f2(j); if a<b n1=n1+1; else if a>b n2=n2+1; end end end if n1>n2 N1=N1+1; else if n1<n2 N2=N2+1; end end end if N1>N2 M1=M1+1; else if N2>N1 M2=M2+1; end end endM1M2問題二程序：clear allclcfor i=618:628,7

19、01:730 fm='F:5dataMEASURING','020',num2str(i),'.SIX' tmp=dlmread(fm); latm=tmp(:,2); lonm=tmp(:,3); rqm=tmp(:,4:7); r=rqm; r(r<0.1)=0; r(r>=0.1 &r<2.5)=1; r(r>=2.5 &r<6)=2; r(r>=6 &r<12)=3; r(r>=12 &r<25)=4; r(r>=25 &r<60)=

20、5; r(r>=60)=6; lonf=dlmread('F:5dataFORECASTlon.dat',''); latf=dlmread('F:5dataFORECASTlat.dat',''); for k=1:4 ff1='F:5dataFORECAST','f',num2str(i),num2str(k),'_dis1' ff2='F:5dataFORECAST','f',num2str(i),num2str(k),'_dis2

21、' rqf1=dlmread(ff1); rqf2=dlmread(ff2); rqm_f1=griddata(lonf,latf,rqf1,lonm,latm,'V4'); rqm_f2=griddata(lonf,latf,rqf2,lonm,latm,'V4'); r1=rqm_f1; r2=rqm_f2; r1(r1<0.1)=0; r1(r1>=0.1 & r1<2.5)=1; r1(r1>=2.5 & r1<6)=2; r1(r1>=6 & r1<12)=3; r1(r1>

22、;=12 & r1<25)=4; r1(r1>=25 & r1<60)=5; r1(r1>=60)=6; r2(r2<0.1)=0; r2(r2>=0.1 & r2<2.5)=1; r2(r2>=2.5 & r2<6)=2; r2(r2>=6 & r2<12)=3; r2(r2>=12 & r2<25)=4; r2(r2>=25 & r2<60)=5; r2(r2>=60)=6; for j=1:91 a1=abs( r1(j)-r(j,k)

23、); a2=abs( r2(j)-r(j,k) ); c1=abs(rqm_f1(j)-rqm(j,k); c2=abs( rqm_f2-rqm(j,k) ); if a1=0 x1(j)=0; s1(j)=1./(1+0.25*x1(j).3.34); else if rqm_f1(j)>=0 & rqm_f1(j)<6 if c1>=0.1 & c1<2.5 x1(j)=1; s1(j)=1./(1+0.25*x1(j).3.34); else if c1>=2.5 & c1<=5.9 x1(j)=2; s1(j)=1./(1+0

24、.25*x1(j).3.34); end end else if rqm_f1(j)>=6 & rqm_f1(j)<25 if c1>=0.1 &c1<6 x1(j)=1; s1(j)=1./(1+0.25*x1(j).3.34); else if c1>=6 & c1<12 x1(j)=2; s1(j)=1./(1+0.25*x1(j).3.34); else if c1>=12 & c1<=18.9 x1(j)=3; s1(j)=1./(1+0.25*x1(j).3.34); end end end else

25、if rqm_f1(j)>=25 if c1>=0.1 & c1<10 x1(j)=1; s1(j)=1./(1+0.25*x1(j).3.34); else if c1>=10 & c1<20 x1(j)=2; s1(j)=1./(1+0.25*x1(j).3.34); else if c1>=20 &c1<30 x1(j)=3; s1(j)=1./(1+0.25*x1(j).3.34); else if c1>=30 &c1<40 x1(j)=4; s1(j)=1./(1+0.25*x1(j).3.34)

26、; else if c1>=40 &c1<60 x1(j)=5; s1(j)=1./(1+0.25*x1(j).3.34); else if c1>=60 x1(j)=6; s1(j)=1./(1+0.25*x1(j).3.34); end end end end end end end end end end if a2=0 x2(j)=0; s2(j)=1./(1+0.25*x2(j).3.34); else if rqm_f2(j)>=0 & rqm_f2(j)<6 if c2>=0.1 & c2<2.5 x2(j)=1; s2(j)=1./(1+0.25*x2(j).3.34); else if c2>=2.5 & c2<=5.9 x2(j)=2; s2(j)=1./(1+0.25*