雨量預報分析的評價模型-數(shù)學建模(共16頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上雨量預報分析的評價模型一、摘要我們將FORECAST文件夾中的數(shù)據(jù)按日期先后順序?qū)隡atlab，建立53×47×164的三維矩陣rain1和rain2；把MEASURING文件夾中的數(shù)據(jù)以同樣方法導入91×7×41的三維矩陣temp中，然后建立循環(huán)將temp矩陣中每一層的后4列提取，另存入一個91×164的rain3矩陣；在命令窗中直接導入預測點的經(jīng)度和緯度存入矩陣lon和lat中，導入實測點的經(jīng)度和緯度存入矩陣lon1和lat1中，并對其作圖，得到實測點和預測點的經(jīng)緯度圖。整理得到91個觀測點41天的預測值和測量值

2、對應的兩個91×164矩陣，根據(jù)氣象部門將降雨的等級分為6個等級的分法，把矩陣中相應的降雨量值轉(zhuǎn)化為其所對應等級值，其中，預測中的零全部記為0，得到兩個預報等級矩陣。針對問題（1），利用插值基點為散亂節(jié)點的插值函數(shù) 1在Matlab中進行三次樣條插值處理，將91個觀測站點41天164個時段的雨量情況進行預測。利用殘差平方和以及平均誤差來作為評價的標準。殘差平方和與平均誤差值較小的一種預測方法作為較好的預報方法。殘差平方和以及平均誤差數(shù)值越小，表明預報越準確度越高。預測方法一的殘差平方和為.00，平均誤差為0.4553。預測方法二的殘差平方和為.00，平均誤差為0.4753。雨量預報方

3、法一的準確性更高一些。針對問題（2），兩個預報等級矩陣，繼續(xù)利用殘差平方和以及平均誤差來作為評價的標準。殘差平方和以及平均誤差數(shù)值越小，表明預報越準確度越高，相應公眾感受就越好。預測方法一的殘差平方和為2774，平均誤差為0.1730。預測方法二的殘差平方和為2806，平均誤差為0.1745。雨量預報方法一的準確性更高一些。由于殘差平方和與平均誤差難以反映真實匯報的準確度，我們將模型改進優(yōu)化。把矩陣中相應的降雨量值轉(zhuǎn)化為其所對應等級值，得到兩個預報等級矩陣，將兩個預報等級矩陣與實測等級矩陣做差值運算，得到兩個等級差矩陣，對等級差作絕對值處理，進行等級差統(tǒng)計。我們利用預測準確度檢驗法對兩種預報進

4、行評價。預測準確度（）等于預報正確次數(shù)（）（即運算之差為0的情況）和預測次數(shù)（）之比，即。準確度越高，表明預報準確度越高，相應公眾感受就越好。預報1的預報準確度為83.26%高于預報2的準確度83.11%，公眾更易接受第一種預報方法。關鍵字：散亂節(jié)點插值 殘差平方和 平均誤差 預報等級矩陣 預測準確度二、問題重述雨量預報對農(nóng)業(yè)生產(chǎn)和城市工作和生活有重要作用，但準確、及時地對雨量作出預報是一個十分困難的問題，廣受世界各國關注。我國某地氣象臺和氣象研究所正在研究6小時雨量預報方法，即每天晚上20點預報從21點開始的4個時段（21點至次日3點，次日3點至9點，9點至15點，15點至21點）在某些位置

5、的雨量，這些位置位于東經(jīng)120度、北緯32度附近的53×47的等距網(wǎng)格點上。同時設立91個觀測站點實測這些時段的實際雨量，由于各種條件的限制，站點的設置是不均勻的。氣象部門希望建立一種科學評價預報方法好壞的數(shù)學模型與方法。氣象部門提供了41天的用兩種不同方法的預報數(shù)據(jù)和相應的實測數(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>_dis2，例如f

6、6181_dis1中包含2002年6月18日晚上20點采用第一種方法預報的第一時段數(shù)據(jù)（其2491個數(shù)據(jù)為該時段各網(wǎng)格點的雨量），而f6183_dis2中包含2002年6月18日晚上20點采用第二種方法預報的第三時段數(shù)據(jù)。MEASURING中包含了41個名為<日期>.SIX的文件，如.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.3000 118.8500

7、 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.125毫米為暴雨，25.160

8、毫米為大暴雨，大于60.1毫米為特大暴雨。若按此分級向公眾預報，如何在評價方法中考慮公眾的感受？三、名詞和符號說明表示預測點的經(jīng)度值表示預測點的緯度值表示預測點已知的預測值表示觀察站點的經(jīng)度值表示觀察站點的緯度值表示三次樣條插值的參數(shù)選項觀察站點的預測值表示兩種預測值矩陣中第i個元素的測量值表示實測矩陣中第i個元素的測量值表示殘差平方和表示均值誤差 表示預測準確度預報正確次數(shù)預測次數(shù)四、模型假設：假設題目中全部數(shù)據(jù)真實可靠，忽略誤差；：假設觀測站所在位置的經(jīng)緯度準確無誤；：假設天氣預報針對的位置在所給網(wǎng)格點附近；：假設雨量在各網(wǎng)點之間的變動是連續(xù)的；五、問題分析針對問題1，我們將兩種預測方法的

9、所有預測值構造成兩個以有序時間段對應的預測值為列,以網(wǎng)格點的個數(shù)為行的2491×164矩陣，對于91觀測站點41天的實測值做同樣的處理，構造成91×164的矩陣。這樣，繁瑣的數(shù)據(jù)經(jīng)過預處理后就整理成了三個矩陣。由于觀測站點相應位置沒有兩種預測方法對應的預測值，無法直接進行評價，我們采用了三次樣條插值的方法進行插值預處理，到了91個觀測站點兩種預測方法的相應時刻的預測值，然后將兩種預測方法雨量預測值與雨量實測值進行比較，從而判斷出兩種預測方法的準確性。針對問題2，我們根據(jù)要求的雨量分級方式來考慮觀眾的感受。我們將問題1中91個觀測站點預測處理后雨量預測值構成的兩個91

10、5;164矩陣和實際雨量觀測值構成的91×164這個三個矩陣分別采用雨量等級記法構造出三個新的矩陣，然后分別把兩個預測值構成的降雨量等級矩陣和觀測值構成的等級矩陣對應元素相減并取絕對值，并進行等級統(tǒng)計，再利用預測準確度檢驗法進行判斷，準確度越高說明我們預報的誤差越小，表明預測方法更準確。六、模型建立1、數(shù)據(jù)預處理（1）針對問題1根據(jù)上面的分析，我們先對數(shù)據(jù)進行預處理。處理方法為：把FORECAST文件夾中的第一種和第二種預測方式得到的數(shù)據(jù)分開兩個文件夾，分別以記事本格式按照日期的先后順序有序的導入Matlab的workspace工作空間中，然后建立m文件編輯公式將兩部分的數(shù)據(jù)導入53

11、×47×164的三維矩陣rain1和rain2中；把MEASURING文件夾中的數(shù)據(jù)以同樣方法導入91×7×41的三維矩陣temp中，然后建立循環(huán)將temp矩陣中每一層的后4列提取，另存入一個91×164的rain3矩陣；在命令窗中直接導入預測點的經(jīng)度和緯度存入矩陣lon和lat中，導入實測點的經(jīng)度和緯度存入矩陣lon1和lat1中，并對其作圖，如圖5-1。實現(xiàn)的matlab語句已呈現(xiàn)在附錄2.1中。圖5-1 預測點（彩色實線）與實測點（藍色孤點）由于三維矩陣無法用表格的形式呈現(xiàn)，我們分別截取了rain1和rain2矩陣的第一層呈現(xiàn)在下表5-1

12、和5-2中，rain3是二維矩陣，將其數(shù)據(jù)呈現(xiàn)在表5-3中：表5-1 預測方法一在6月18日第一個時間段的預測值構造成的矩陣預測值B1B2B3.B47A10.03100.03210.0313.0.0081A20.03020.03100.0315.0.0078A30.03400.03410.0337.0.0079.A520.00880.00890.0084.0.0061A530.00890.00860.0086.0.0062表5-2 預測方法二在6月18日第一個時間段的預測值構造成的矩陣預測值B1B2B3.B47A10.03190.02930.0322.0.0077A20.02880.03110

13、.0303.0.0085A30.03250.03230.0313.0.0082.A520.00950.00890.0079.0.0056A530.00850.00940.0093.0.0057表5-3 91個站點在41天共164個時段的雨量實測值預測值D1D2D3.D164P1000.0P2000.0P3000.0.P90004.1.0P91000.2.0（2）針對問題2由問題1我們可以整理得到91個觀測點41天的預測值和測量值對應的兩個91×164矩陣，再根據(jù)問題2中氣象部門將降雨的等級分為6個等級的分法，把矩陣中相應的降雨量值轉(zhuǎn)化為其所對應等級值，其中，預測中的零全部記為0，得到

14、兩個預報等級矩陣，如下表5-4和5-5：表5-4 預測方法一在6月18日第一個時間段的預報等級構造成的矩陣預報等級B1B2B3.B47A1000.0A2000.0A3000.0.A52000.0A53000.0表5-5 預測方法二在6月18日第一個時間段的預報等級構造成的矩陣預報等級B1B2B3.B47A1000.0A2000.0A3000.0.A52000.0A53000.02、模型建立與求解（1）針對問題1由于91個觀測站點沒有相應的預測值，因此不能夠直接對實測值進行評價，屬于離散的散亂節(jié)點，我們利用插值基點為散亂節(jié)點的插值函數(shù)1在Matlab中進行三次樣條插值處理，插值函數(shù)為： 1 其中

15、表示預測點的經(jīng)度值，表示預測點的緯度值，表示預測點的已知的預測值，表示觀察站點的緯度值，表示觀察站點的經(jīng)度值，表示三次樣條插值的參數(shù)選項，觀察站點的預測值。將91個觀測站點41天164個時段的雨量情況進行預測之后，我們可以建立模型來評價這兩種預測方法。這里我們利用殘差平方和2以及平均誤差3來作為評價的標準： 2 3最后，我們根據(jù)和的值進行評價，取值越大，表明預報的值準確性越低。因此，殘差平方和與平均誤差值較小的一種預測方法作為較好的預報方法。利用公式（1），我們在Matlab中應用編程求解，程序代碼見附錄2.2。求解之后得到91個觀測點41天164個時段的預測值，整理成91×164矩

16、陣，然后把預測矩陣和實測矩陣對應元素值相減取平方作殘差平方和，再作平均誤差，最后結果如下表5-6：表5-6 雨量預測計算結果數(shù)據(jù)表結果殘差平方和平均誤差預測方法一.000.4553預測方法二.000.4753由此可知，雨量預報方法一的準確性更高一些。（2）針對問題2由問題1我們可以整理得到91個觀測點41天的預測值和測量值對應的兩個91×164矩陣，再根據(jù)問題2中氣象部門將降雨的等級分為6個等級的分法，把矩陣中相應的降雨量值轉(zhuǎn)化為其所對應等級值，其中，預測中的零全部記為0，得到兩個預報等級矩陣，繼續(xù)利用殘差平方和以及平均誤差來作為評價的標準。殘差平方和以及平均誤差數(shù)值越小，表明預報越

17、準確度越高，相應公眾感受就越好。表5-7 雨量等級計算結果數(shù)據(jù)表結果殘差平方和平均誤差預測方法一27740.1730預測方法二28060.1745由結果可知，預測方法一的準確度更高。七、模型優(yōu)化針對問題2，由于殘差平方和與平均誤差難以反映真實匯報的準確度，我們將模型改進優(yōu)化。由問題1我們可以整理得到91個觀測點41天的預測值和測量值對應的兩個91×164矩陣，再根據(jù)問題2中氣象部門將降雨的等級分為6個等級的分法，把矩陣中相應的降雨量值轉(zhuǎn)化為其所對應等級值，得到兩個預報等級矩陣，將兩個預報等級矩陣與實測等級矩陣做差值運算，得到兩個等級差矩陣，對等級差作絕對值處理后，我們就可以從中進行等

18、級差統(tǒng)計。我們利用預測準確度檢驗法對兩種預報進行評價。預測準確度（）等于預報正確次數(shù)（）（即運算之差為0的情況）和預測次數(shù)（）之比，即 4準確度越高，表明預報越準確度越高，相應公眾感受就越好。雨量分級與統(tǒng)計程序見附錄4，數(shù)據(jù)統(tǒng)計結果見下表7-1：表7-1 兩種預報等級與實測等級的級差與比例表雨量等級差預報1等級差個數(shù)預報1等級差個數(shù)比例預報2等級差個數(shù)預報2等級差個數(shù)比例0124250.8326124040.8311121220.142221280.14262550.0037520.00353100.0007100.0007410.000120.00015000060000由表7-1可見預報1

19、的預報準確度為83.26%高于預報2的準確度83.11%。預報1與實測雨量的符合度高，那么兩種預報方法中，公眾更易接受第一種預報方法。八、模型的評價與推廣優(yōu)點：本文思路自然流暢，較多地使用了原始數(shù)據(jù)，對91個站點的預測值比較準確。缺點：認為降雨量在不同的地區(qū)之間的變化有連續(xù)性，過于理想化。推廣：本文對降雨量的預測方法也可以運用在工農(nóng)業(yè)生產(chǎn)中其他方面，比如一段時期某地區(qū)的氣溫變化、全國糧食產(chǎn)量預測等。本文應用的散亂節(jié)點插值和殘差分析法可以應用在很多方面。九、附錄附錄1：數(shù)據(jù)初始化>> File_FORECAST1=dir('預測方法一')File_FORECAST1

20、=166x1 struct array with fields: name date bytes isdir>> File_FORECAST2=dir('預測方法二')File_FORECAST2 = 166x1 struct array with fields: name date bytesisdir>> File_MEASURING=dir('MEASURING')File_MEASURING = 43x1 struct array with fields: name date bytes isdir【mypro1.m】rain1=

21、zeros(53,47,164);for n=1:164 filename=File_FORECAST1(n+2).name; rain1(:,:,n)=importdata(filename);end【mypro2.m】rain2=zeros(53,47,164);for n=1:164 filename=File_FORECAST2(n+2).name; rain2(:,:,n)=importdata(filename);end【mypro3.m】temp=zeros(91,7,41);for n=1:41 filename=File_MEASURING(n+2).name; temp(:

22、,:,n)=importdata(filename);endrain3=zeros(91,164);i=1;for j=1:41 for k=4:7 rain3(:,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:164 wea=rain1(:,:,i); wear=rain

23、3(:,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);endpjwc=pjwc/(91*164);pjwcccpfh【calc2.m】lat;lon;lat1;lon1;ccpfh=0;pjwc=0;for i=1:164 wea=rain2(:,:,i); wear=rain3(:,i); weap=griddata(lat,lon,wea,

24、lat1,lon1); chazhi=wear-weap; pingfang=chazhi.2; juedui=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

25、);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');endfor i=1:164 for j=1:91 if weap1(j,i)<0.1 rank1(j,i)=0; elseif weap1(j,

26、i)<=2.5 rank1(j,i)=1; elseif weap1(j,i)<=6.0 rank1(j,i)=2; elseif weap1(j,i)<=12.0 rank1(j,i)=3; elseif weap1(j,i)<=25.0 rank1(j,i)=4; elseif weap1(j,i)<=60.0 rank1(j,i)=5; else rank1(j,i)=6; end endendfor i=1:164 for j=1:91 if weap2(j,i)<0.1 rank2(j,i)=0; elseif weap2(j,i)<=2.5

27、 rank2(j,i)=1; elseif weap2(j,i)<=6.0 rank2(j,i)=2; elseif weap2(j,i)<=12.0 rank2(j,i)=3; elseif weap2(j,i)<=25.0 rank2(j,i)=4; elseif weap2(j,i)<=60.0 rank2(j,i)=5; else rank2(j,i)=6; end endendfor i=1:164 for j=1:91 if rain3(j,i)<0.1 rank(j,i)=0; elseif rain3(j,i)<=2.5 rank(j,i)=

28、1; elseif rain3(j,i)<=6.0 rank(j,i)=2; elseif rain3(j,i)<=12.0 rank(j,i)=3; elseif rain3(j,i)<=25.0 rank(j,i)=4; elseif rain3(j,i)<=60.0 rank(j,i)=5; else rank(j,i)=6; end endendchazhi1=rank1-rank;chazhi2=rank2-rank;pingfang1=chazhi1.2;pingfang2=chazhi2.2;juedui1=abs(chazhi1);juedui2=abs(chazhi2);ccpfh1=ccpfh1+sum(sum(pingfang1,1),2);ccpfh2=ccpfh2+sum(sum(pingfang2,1),2);pjwc1=pjwc1+sum(sum(juedui1),2);pjwc2=pjwc2+sum(sum(juedui2),2);pjwc1=pjwc1/(91*164);pjwc2=pjwc2/(91*16