雨量預報數(shù)學模型及求解代碼

1、僅供個人學習參考賽區(qū)評閱編號（由賽區(qū)組委會填寫）：2015高教社杯全國大學生數(shù)學建模競賽承諾書我們仔細閱讀了全國大學生數(shù)學建模競賽章程和全國大學生數(shù)學建模競賽參賽規(guī)則（以下簡稱為“競賽章程和參賽規(guī)則”，可從全國大學生數(shù)學建模競賽網(wǎng)站下載）。我們完全明白，在競賽開始后參賽隊員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊外的任何人（包括指導教師）研究、討論與賽題有關的問題。我們知道，抄襲別人的成果是違反競賽章程和參賽規(guī)則的，如果引用別人的成果或其他公開的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻的表述方式在正文引用處和參考文獻中明確列出。我們鄭重承諾，嚴格遵守競賽章程和參賽規(guī)則，以

2、保證競賽的公正、公平性。如有違反競賽章程和參賽規(guī)則的行為，我們將受到嚴肅處理。我們授權全國大學生數(shù)學建模競賽組委會，可將我們的論文以任何形式進行公開展示（包括進行網(wǎng)上公示，在書籍、期刊和其他媒體進行正式或非正式發(fā)表等）。我們參賽選擇的題號（從A/B/C/D中選擇一項填寫）：我們的報名參賽隊號（12位數(shù)字全國統(tǒng)一編號）：參賽學校（完整的學校全稱，不含院系名）：參賽隊員（打印并簽名）：1張偉堯2熊霜3韓倩指導教師或指導教師組負責人（打印并簽名）：董銀麗日期：2015年9月6日賽區(qū)評閱編號（由賽區(qū)組委會填寫）：-2015高教社杯全國大學生數(shù)學建模競賽編號專用頁賽區(qū)評閱記錄（可供賽區(qū)評閱時使用）：全國

3、評閱隨機編號（由全國組委會填寫）：雨量預報方法的評價模型摘要本文建立了一個關于雨量預報方法的評價模型。首先，通過給定的大量數(shù)據(jù)（預測數(shù)據(jù)和實測數(shù)據(jù)）利用Matlab圖形處理功能的基本繪圖命令plot畫出2491個網(wǎng)格點和91個觀測點的位置。在可接受范圍內(nèi)，計算各觀測站點和等距網(wǎng)格點之間的距離，并按升序排序。取其前5個到觀測站點距離最小的等距網(wǎng)格點，再根據(jù)歐拉公式距離倒數(shù)加權的方法對它們賦權重，取其前5個網(wǎng)格點的雨量，分別乘以它們各自對應的權并求和，就是相對應觀測站點的雨量。分別對兩種方法預報的41天每天4個時段各網(wǎng)格點的雨量處理，最后對兩種方法得到站點的預測雨量進行分析。問題一：利用上述得到的

4、站點預測雨量，計算某個觀測點在某時段的預報偏差率，并對其預報偏差率求和再算術平方根得出標準率，作為評價準確性高低的指標，從而得到第一種雨量預報方法的準確率為S二104.5506，第二種雨量預報方法準確率為1S=104.6044。因此得出第一種方法比第二種方法的準確性高。2問題二：分別將兩種方法的站點預測雨量與實測雨量按照題目給定的等級劃分，得到等級矩陣，然后分別用兩種方法的預測等級矩陣和實測等級矩陣作差，統(tǒng)計同一級別、相差1級、相差2級相差6級數(shù)據(jù)的頻數(shù)和頻率。等級差越小，相對報錯率越低，公眾越滿意。因此得出第一種方法比第二種方法的公眾滿意度高。關鍵詞：歐拉公式權標準率等級差一、問題重述雨量預

5、報對農(nóng)業(yè)生產(chǎn)和城市工作和生活有重要作用，但如何準確、及時的對雨量做出預報是一個十分困難的問題。本題由氣象部門提出，希望建立一種科學評價預報方法好壞的數(shù)學模型與方法。我國某地氣象臺和氣象研究所正在研究6小時雨量預報方法，即每天晚上20點預報從21點開始的4個時段（21點至次日3點，次日3點至9點，9點至15點，15點至21點）在某些位置的雨量，這些位置位于東經(jīng)120度，北緯32度附近的53x47的等距網(wǎng)格點上。同時設立91個分布不均勻的觀測站點實測這些時段的實際雨量。氣象部門提供了41天用兩種不同方法的預報數(shù)據(jù)和相應的實測數(shù)據(jù)，預測數(shù)據(jù)在文件夾FORECAST中，實測數(shù)據(jù)在文件夾MEASURIN

6、G中，其中的文件都可以用Windows系統(tǒng)的“寫字板”程序打開閱讀。經(jīng)（lon.dat）緯（lat.dat）度也分別包括在文件夾FORECAST中，其余文件名為f日期i_dis1和f日期i_dis2表示該日期用的第一或第二種方法。題目中有如下假設：1、雨量用毫米作單位，小于0.1毫米的視為無雨。2、氣象部門將6小時降雨量分為6等：0.12.5毫米為小雨，2.66毫米為中雨，6.112毫米為大雨，12.125毫米為暴雨，25.160毫米為大暴雨，大于60.1毫米為特大暴雨。題目的問題如下：問題一：請建立數(shù)學模型來評價兩種6小時雨量預報方法的準確性；問題二：若按此分級向公眾預報，如何在評價方法中考

7、慮公眾的感受？二、問題分析我們從題目中了解分析得到：氣象臺每天晚上20點預報從21點開始的4個時段（21點至次日3點，次日3點至9點，9點至15點，15點至21點）在某些位置的雨量，這些位置位于東經(jīng)120度，北緯32度附近的5347的等距網(wǎng)格點上。同時設立91個分布不均勻的觀測站點實測這些時段的實際雨量。由于網(wǎng)格點比較多，且每個網(wǎng)格點的位置是以經(jīng)度和緯度表示在一定的區(qū)域，所以我們將緯度看作X軸，經(jīng)度看作y軸，采用Matlab圖形處理功能的基本繪圖命令plot畫出散點圖（圖1），（程序見附錄一）。圖中綠色表示網(wǎng)格結點的分布，紅色表示實測站點的分布。從圖中可以分析看出，氣象部門提供了2491個網(wǎng)格

8、點上41天4個時間段的大量預報數(shù)據(jù)（雨量），并且同樣給出了91個觀測站點的實測數(shù)據(jù)（雨量）。所以我們想通過網(wǎng)格上的預報數(shù)據(jù)來預測實測站點的數(shù)據(jù)。然而，觀測站點集中在所有網(wǎng)格的中央部分，而四周是大量的距離比較遠的網(wǎng)格點。因此，通過搜索出2491個網(wǎng)格點中對站點影響比較大的幾個網(wǎng)格點，再用搜索出來的幾個網(wǎng)格點的預測數(shù)據(jù)加權求出一個預測數(shù)據(jù)（雨量），進而和該站點實測數(shù)據(jù)進行比較，來評價兩種6小時雨量預報方法的準確性。在向公眾預報時，采用一種合理、準確的預測方法，增加雨量分等級預報的同級率，能對公眾的出行起到良好的指導作用，使人們對雨量預報有更深的理解，更多的關注。125124123122度121經(jīng)1

9、2011911811727282930313233343536緯度圖1、網(wǎng)格結點與實測站點的散點分布圖三、模型假設耀莓左；lymi拆1、觀測站點的設置是不均勻的；2、題中的網(wǎng)格式等間距的正方形網(wǎng)格(所謂“正方形網(wǎng)格”是指每個格子都是正方形的網(wǎng)格；網(wǎng)絡點是指縱線和橫線的交叉點)；3、一個x軸、y軸分別是緯度和經(jīng)度的坐標，通過把點的緯度和經(jīng)度分別看作橫坐標和縱坐標，用歐拉公式計算d=.;(a-a)2+(b-b)2來作為兩點之間的距離。ijij4、點到觀測站點的距離越短，則對觀測站點的雨量影響越大；5、單個網(wǎng)格點到觀測站點倒數(shù)與所取的5個網(wǎng)格點到觀測點倒數(shù)之和的比為它的權值；6、雨量用毫米做單位，小

10、于0.1毫米視為無雨；四、符號說明Q二(q)(i=1,2.,91,n二1,2,3,4,5):與第i個觀測站點的距離最小的前5個網(wǎng)格點ini.n的對應權矩陣；D二(d)(i二1,2,.91,n二1,2,3,4,5):與第i個觀測站點的距離最小的前5個網(wǎng)格點ini.n的距離矩陣；dx-x)2+(y-y)2():歐拉公式計算第i個觀測站點到第j個實測站點的距vijij離計算；p(x,y)(i二1,2,.91):第i個觀測站點的坐標點的緯度和經(jīng)度；iii8:可接受度數(shù)差；qin第i個實測站點的距離最小的前5個網(wǎng)格點-!l(i二l,2,.91,n二1,2,3,4,5)：f丄djTij到第i個觀測站點的權

11、重;y=fqf(iijininn=1=l,2,.,91,n=1,2,3,4,5):第i個觀測站點在第j個時段的預測雨量之和；M=(y)(i=1,2,.91,j=12.164):第i個觀測站點在第j個時段的預測雨量矩陣;ijijijI丄)(i=1,2,.91,j=1,2,.164):第i個觀測站點在第j個時段的實際雨量矩陣；Q=ijysijyijj(i=12.91,j=1,2,.164):第i個觀測站點在第j個時段的預報偏差率;164(q)(i=1,2,.,91,j=1,2,.,164):預報偏差率的算術平方根;i=1j=1J,J,.,J(c=1,2):第c種方法的預報數(shù)據(jù)與實測數(shù)據(jù)處在同一級別

12、，相差同一c0c1c6級別，相差1級相差6級的頻數(shù)五、模型的建立與求解5、1問題一模型建立及其求解由于問題涉及的數(shù)據(jù)量較大，為了模型計算的簡便，需要對數(shù)據(jù)進行處理。(一)數(shù)據(jù)處理：為了評價兩種6小時雨量預報方法的準確性，我們采用網(wǎng)格點上的預報數(shù)據(jù)來預測觀測站點的數(shù)據(jù)，再來和實際測得的數(shù)據(jù)相比，判斷其準確性，首先從預測數(shù)據(jù)選取可接受范圍內(nèi)的網(wǎng)格點，其次選取實測點附近5個距離最小的預測點，計算5個網(wǎng)格點分別到觀測站點的權重，并計算觀測站點的預測雨量。(1)可接受范圍內(nèi)網(wǎng)格點的選取以坐標p(a,b)(i二1,2,.,91)為基準點，給定一個可接受度數(shù)差(在求解中取iii=0.2，對任意的p(i二1,

13、2,.,91),搜索其任一個觀測站點在緯度和經(jīng)度都上下增加的i正方形內(nèi)所有等距網(wǎng)格點。若網(wǎng)格點x(m,n)緯度和經(jīng)度在同時滿足a-ma+和iib-nb+時，即認為該網(wǎng)格點時可接受范圍內(nèi)的網(wǎng)格點。ii(2)網(wǎng)格點到站點距離的計算找到可接受范圍內(nèi)的網(wǎng)格點后，我們計算網(wǎng)格點和這些觀測站點間的距離。再得到個觀測站點和等距網(wǎng)格點之間的距離后，將各觀測站點按距離從小到大排序后保存，我們用Matlab編程求得結果(程序見附錄二)。(3)權的計算各觀測站點和等距網(wǎng)格點之間的距離從小到大排序后，為了更好地用網(wǎng)格點的預報雨量來預測觀測站點的雨量，我們?nèi)∏?個到觀測站點距離最小的等距網(wǎng)格點。根據(jù)歐拉公式距離的倒數(shù)加

14、權的方法，先算出前5個網(wǎng)格點到網(wǎng)格點的距離，再分別對它們求倒，則5個網(wǎng)格點分別到觀測站點的權重為它們之間距離的倒數(shù)。1權的計算公式為：qindj=1ij-in,(i二1,2,.,91,n二1,2,3,4,5)4)觀測站點預測雨量的計算為了預測各觀測站點在某月某日某個時段的雨量值，我們采用距離的倒數(shù)加權的方法，取出5個等距網(wǎng)格點分別在某月某日某個時段的雨量值，然后分別乘以它們各自對觀測站點的權重，再求和就為預測降雨量。預測降雨量的計算公式：y=qf；ijininn=1每個時段中，每個觀測點對應有5個網(wǎng)格點預測雨量值，可計算出1個觀測站點預測雨量值，91個觀測站點，164個時段就可計算出一個第i個

15、觀測站點的第j個時段預測降雨量矩陣M=(y)(i=j=1,2,.,164)。iji.j這里我們用Matlab編程求得兩個方法對應的預測降雨量矩陣(程序見附錄三)。(二)模型建立及求解：將兩個方法對應的預測降雨量矩陣分別與實測降雨量矩陣進行比較，分析出哪一個的準確性高。這里我們用計算出預報偏差率的算術平方根作為一個準確性指數(shù)，來辨別準確性的高低。預測降雨量的矩陣M=(y)(i=j=1,2,.,164)和實際降雨量矩陣iji.jI=(s)(i=1,2,.,91,j=1,2,.,164)。則有:iji.j預報偏差率計算公式：b=九為(i=1,2,.,91,j=1,2,.,164)；ijyij取二sh

16、ice_x2&x二shice_x1);%滿足條件的點的x坐標位置xb2=x(b1);yb2=y(b1);plot(xb2,yb2,.g,xi,yi,.r);b3=find(yb2=shice_y2&yb2二shice_y1);%篩選后滿足條件的點的y坐標位置xb3=xb2(b3);%有效的點yb3=yb2(b3);plot(xb3,yb3,.g,xi,yi,.r);axis(29.334.3117.8122.8);xb4=xi;yb4=yi;l=length(xb3);d=zeros(91,l);fori=1:91forj=1:ld(i,j)二sqrt(xb4(i)-xb3(j)廠2+(yb4

17、(i)-yb3(j)廠2);%計算距離endend附錄三觀測站點預測降雨量的計算B,C=sort(d,2);%對d從小到大排序，并記錄在原矩陣中的坐標D=C(:,1,2,3,4,5);E=B(:,1,2,3,4,5);F=1./E;F(:,6)=sum(F,2);r,l=size(F);fori=l:rK(i,l:l)=F(i,:)/max(max(F(i,:);%每個最小值點權的計算endH=K(:,1,2,3,4,5);H(find(isnan(H)=l)=l%將H中NaN替換為1HH二reshape(H,455,l);%前五個最小值的權并排序a=bl(b3(D);%返回91*5個點在x中

18、的矩陣位置%預測數(shù)據(jù)調(diào)用和降雨量的篩選%將CurrentDirectory路徑改為數(shù)據(jù)所在位置%A1=;A2=;formonth1=6:7ifmonth1=6;fordata1=18:28fortime1=1:4str1二fint2str(month1)int2str(data1)int2str(time1)_dis1;str2=fint2str(month1)int2str(data1)int2str(time1)_dis2;tmp1=load(str1);tmp2=load(str2);A1=A1;tmp1(:);A2二A2;tmp2(:);endendelsefordata1=1:30i

19、fdata110strtmp=0int2str(data1);elsestrtmp=int2str(data1);endfortime1=1:4str1=fint2str(month1)strtmpint2str(time1)_dis1;str2=fint2str(month1)strtmpint2str(time1)_dis2;tmp1=load(str1);tmp2=load(str2);A1=A1;tmp1(:);A2=A2;tmp2(:);endendendend%A1為第一種方法統(tǒng)計的所有數(shù)據(jù)，A2為第二種方法統(tǒng)計的所有數(shù)據(jù)w二reshape(a,455,l);%將91*5個點轉換成

20、455*1All二Al(:,w);%第一類164*455個降雨量A22=A2(:,w);%第二類164*455個降雨量gg=repmat(HH,164,1)%為使HH與All，A22大小相同，擴展平鋪HH,gg大小為164*455A1_1=gg.*A11%第一類方法權與降雨量乘積A2_2=gg.*A22%第二類方法權與降雨量乘積大小為164*455%分別對A1_1，A2_2每隔5列求和并返回sum1,sum2%sum1=;fori=1:5:454sum1二sum1A1_1(:,i)+A1_1(:,i+1)+A1_1(:,i+2)+A1_1(:,i+3)+A1_1(:,i+4);endsum2=

21、;fori=1:5:454sum2二sum2A2_2(:,i)+A2_2(:,i+1)+A2_2(:,i+2)+A2_2(:,i+3)+A2_2(:,i+4);end附錄四模型的求解%91個站實測數(shù)據(jù)的導入%將CurrentDirectory路徑改為數(shù)據(jù)所在位置B1=;formonth1=6:7ifmonth1=6;fordata1=18:28%str3=020int2str(month1)int2str(data1);str3=020int2str(month1)int2str(data1).SIX;tmp3=load(str3);B1=B1;tmp3(:,4:7);endelsefordata1=1:30ifdata110strtmp=0int2str(data1);elsestrtmp=int2str(data1);endstr3=020int2str(month1)strtmp.SIX;tmp3=load(str3);B1=B1;tmp3(:,4:7);endendendB_1二abs(suml-Bl)./suml;%第一種方法的偏差率B_2=abs(sum2-Bl)./sum2;%第二種方法的偏差率B_l(find(isnan(B_l)=l)=O;%將B_1中NaN替換成為0B_2(find(isnan(B_2)=l)=0;%將B_2中NaN替換成為0B_11