關(guān)于雨量預(yù)報方法評價的模型.doc

關(guān)于預(yù)報雨量方法評價的數(shù)學(xué)模型摘要：降雨是日常生活中最常見的天氣現(xiàn)象，及時準(zhǔn)確地對降雨量作出預(yù)報是一個十分困難的問題。題目給出了用兩種降雨量預(yù)測方法預(yù)測得出的41天各時段網(wǎng)格結(jié)點的預(yù)報數(shù)據(jù)，我們對預(yù)報方法進(jìn)行評價。先使用二元三次樣條插值法，對91個站點位置進(jìn)行插值，從而求出各站點預(yù)報數(shù)據(jù)，然后用方差分析計算兩種預(yù)測方法每一天的擬合優(yōu)度，從而評價兩種預(yù)報方法的準(zhǔn)確性。第一種方法擬合優(yōu)度的平均值為0.9087,第二種方法為0.9054,可見這兩種方法的預(yù)測準(zhǔn)確性都比較高，但第一種比第二種的精確度更高。在評價方法中考慮公眾的感受，即考慮公眾對預(yù)報結(jié)果的滿意率和滿意度。本文建立了兩個模型，即滿意率模型和滿意度模型。滿意率指滿意數(shù)與總數(shù)的比率；滿意度指達(dá)到公眾滿意的程度。兩個模型分別運(yùn)用算術(shù)平均數(shù)和加權(quán)平均數(shù)的方法以及概率的知識進(jìn)行求解，第一種方法的，第二種方法的，可見公眾對這兩種方法的預(yù)測值都比較滿意，但第一種方法給公眾的感受優(yōu)于第二種。關(guān)鍵詞：二元三次樣條插值法 擬合優(yōu)度 滿意率 滿意度1、問題的重述本題是一道由氣象部門提出，期望建立一種科學(xué)評價雨量預(yù)報方法好壞的數(shù)學(xué)模型與方法。在本題中，氣象部門所采用了兩種方案研究6小時雨量預(yù)報方法。所謂6小時雨量預(yù)報方法是指每天晚上20點預(yù)報從21點開始的4個時段（21點至次日3點，次日3點至9點，9點至15點，15點至21點）在某些位置的雨量，這些位置位于東經(jīng)120度、北緯32度附近的5347的等距網(wǎng)格點上。同時設(shè)立91個觀測站點實測這些時段的實際雨量，由于各種條件的限制，站點的設(shè)置是不均勻的。氣象部門提供了41天的用兩種不同方法的預(yù)報數(shù)據(jù)和相應(yīng)的實測數(shù)據(jù)。預(yù)報數(shù)據(jù)在文件夾FORECAST中，實測數(shù)據(jù)在文件夾MEASURING中。FORECAST中的文件lon.dat和lat.dat分別包含網(wǎng)格點的經(jīng)緯度，其余文件名為_dis1和_dis2，例如f6181_dis1中包含2002年6月18日晚上20點采用第一種方法預(yù)報的第一時段數(shù)據(jù)（其2491個數(shù)據(jù)為該時段各網(wǎng)格點的雨量），而f6183_dis2中包含2002年6月18日晚上20點采用第二種方法預(yù)報的第三時段數(shù)據(jù)。MEASURING中包含了41個名為.SIX的文件，如020618.SIX表示2002年6月18日晚上21點開始的連續(xù)4個時段各站點的實測數(shù)據(jù)（雨量）。題目中有如下假設(shè)：第一，雨量用毫米做單位，小于0.1毫米視為無雨。第二，氣象部門將6小時降雨量分為6等：0.12.5毫米為小雨，2.66毫米為中雨，6.112毫米為大雨，12.125毫米為暴雨，25.160毫米為大暴雨，大于60.1毫米為特大暴雨。該題目的問題如下：問題一：請建立數(shù)學(xué)模型來評價兩種6小時雨量預(yù)報方法的準(zhǔn)確性；問題二：若上述分級向公眾預(yù)報，如何在評價方法中考慮公眾的感受？2、問題的分析2.1 問題一的分析要評價兩種預(yù)報方法的準(zhǔn)確性，關(guān)鍵是需要把題目中給出的大量預(yù)報數(shù)據(jù)進(jìn)行處理，使其能與各站點實測值進(jìn)行比較。因此我們分為兩個大的步驟處理該問題：第一，根據(jù)觀測點的分布，將網(wǎng)格結(jié)點預(yù)報數(shù)據(jù)轉(zhuǎn)化成觀測站點預(yù)報數(shù)據(jù)；第二，根據(jù)網(wǎng)格結(jié)點預(yù)報數(shù)據(jù)插值成的各站點預(yù)報值與實測值比較，判斷兩種方法的優(yōu)劣。在該問題中所使用的數(shù)學(xué)工具和方法有：第一，根據(jù)各網(wǎng)格結(jié)點的降雨量采用二元三次樣條插值法求出各觀測站點位置的降雨量；第二，通過擬合優(yōu)度分析，判斷出兩種雨量預(yù)報方法的準(zhǔn)確程度。2.2 問題二的分析在評價方法中考慮公眾的感受，即公眾對預(yù)報結(jié)果準(zhǔn)確度的滿意率和滿意度。若要求總體的滿意率和滿意度，首先應(yīng)該確定每一等級差下站點個數(shù)(n)與總站點個數(shù)的比率,通過概率的基本性質(zhì)與平均數(shù)來建立比率與滿意度的函數(shù)關(guān)系，最終得出總體的滿意率和滿意度。3、基本假設(shè)（1） 題目中所給的經(jīng)緯度范圍內(nèi)的地球面是平面；（2） 題目中給出的每個降雨量數(shù)值是相應(yīng)經(jīng)緯度點數(shù)據(jù)，觀測站點可以抽象為一個點，其面積大小不影響數(shù)據(jù)觀測；（3） 降雨量小于0.1毫米的視為無雨；4、符號說明 樣本的均值； 樣本中第點的函數(shù)值； 回歸直線上第點的函數(shù)值； 樣本中第點的函數(shù)值與回歸直線上第點的函數(shù)值的差值； 的總方差，即樣本中各點的函數(shù)值與樣本的均值差的平方和； 的剩余變差，即樣本中各點的函數(shù)值與回歸直線上各點函數(shù)值差的平方和； 的解釋變差，即回歸直線上各點函數(shù)值與樣本均值差的平方和； 擬合優(yōu)度，即對的回歸方程所能解釋的那部分變差占總變差的比例 擬合優(yōu)度的平均值 擬合優(yōu)度的方差 第一種方法每一天的擬合優(yōu)度 第二種方法每一天的擬合優(yōu)度 站點處的預(yù)報數(shù)據(jù)與其實測數(shù)據(jù)的等級差5、模型的建立與求解5.1 問題一將球面問題轉(zhuǎn)化成平面問題解決。5.1.1 將91個觀測站點及2491個網(wǎng)格結(jié)點的經(jīng)緯度聯(lián)系起來，用Matlab軟件編程畫圖，如下圖1所示：“.”表示2491個網(wǎng)格結(jié)點的位置，“*”表是91個站點位置；圖1 原始網(wǎng)格結(jié)點及站點圖通過觀察發(fā)現(xiàn)，91個測量站點的位置都大約集中在經(jīng)度118123度、緯度29.534度之間的范圍內(nèi)。問題一中要求對兩種預(yù)報數(shù)據(jù)的方法進(jìn)行評價，實際是91個站點位置的預(yù)報數(shù)據(jù)值與相應(yīng)站點實測值的比較。由此，可選取上述范圍作為工作區(qū)，將對結(jié)果無影響的網(wǎng)格結(jié)點排除出考慮范圍，即只對在經(jīng)度118123度、緯度29.534度的范圍內(nèi)進(jìn)行數(shù)據(jù)的處理和評價，截取圖1，如圖2所示：圖2 處理后的網(wǎng)格結(jié)點與站點圖再用二元三次樣條插值法求各觀測站點的降雨量，具體如下：二維三次樣條插值法可以直接利用網(wǎng)格結(jié)點的數(shù)值進(jìn)行插值，但在網(wǎng)格結(jié)點處不能得到連續(xù)的二階導(dǎo)數(shù)，其插值函數(shù)式可表示為：其中，為采樣函數(shù)，為插值函數(shù)，（，）為插值點，為的在各點的函數(shù)值，為采樣間隔，為采樣函數(shù)的個數(shù)，、是插值內(nèi)核。根據(jù)此插值法，我們運(yùn)用Matlab軟件中的二維三次樣條函數(shù)的命令（interp2（x，y，z，xj，yj，spline）可以插值出91個站點41天每個時段的降雨量。我們把41天中91個站點的各時段預(yù)測數(shù)據(jù)(降雨量)保存在41個名為_1/2文件中，如f020703_1表示使用第一種預(yù)報方法的數(shù)據(jù)插值得到2002年7月03日晚上21點開始的連續(xù)4個時段各站點的預(yù)測數(shù)據(jù)。這些文件的數(shù)據(jù)格式如表1所示：表1文件_1/2的數(shù)據(jù)格式站號緯度經(jīng)度第1段第2段第3段第4段5813832.9833118.51670.00000.000064.40001.62145813933.3000118.85000.00000.000031.10001.99885814133.6667119.26670.00000.00000.00000.64005814333.8000119.80000.00000.80006.20000.3382注：文件f020703_1中的完整數(shù)據(jù)參見附錄表。將各個站點預(yù)測數(shù)據(jù)與實測數(shù)據(jù)進(jìn)行比較，觀察其插值效果。比如將7月03日晚上21點開始的第四個時段各個站點的預(yù)測降雨量與這個時段各個站點的實測降雨量進(jìn)行比較，插值效果圖如圖3所示(其中紅“*“代表的預(yù)測降雨量，綠“o“代表的實測降雨量)：圖3 各個站點插值效果圖由圖3,可以看出預(yù)測降雨量與實測降雨量相差比較小，預(yù)測降雨量比較符合實際。5.1.2 所謂擬合優(yōu)度原理，它是關(guān)于估計回歸直線與數(shù)據(jù)之間擬合程度的一個很有用的度量，殘差大說明擬合效果欠佳，反之，則擬合的好。殘差的方差除以Y的偏差平方和可以用來衡量擬合優(yōu)度（偏差平方和），,為了計算偏差平方和，可以將等式兩邊同時平方，且對所有觀測值求和：用最小二乘殘差的兩個性質(zhì)和（為殘差），可得： Y的總變差 Y的剩余變差 Y的可解釋變差 是的總變差中，對的回歸方程所能解釋的那一部分變差占總變差的比例，的值大則回歸直線擬合的好，反之則不好。針對此題，我們通過回歸變差，總變差求得,再求出每種方法下的平均值和方差，如表2所示：表2：擬合優(yōu)度表日期6180.98080.98366190.96190.95626200.99990.99996210.73190.70186220.97980.98486230.97710.97556241.00011.00016250.99670.99676260.82250.76306270.99920.99906280.99870.99877011.00001.00007020.33760.36467030.99821.00357040.99460.99457050.90600.89997060.98580.98697071.00001.00007080.97850.98467091.00001.0000710NaNNaN711NaNNaN712NaNNaN7130.38360.36957140.99660.99897150.86600.76537160.71370.73107170.90530.92697181.00001.00007190.65610.55597200.87230.90177210.99950.99997220.99950.99567230.80930.82387240.89110.89957250.76940.82777260.89340.90897270.96100.98847280.94330.94207290.95460.97117300.99470.92050.90870.90540.02330.0249從上表數(shù)據(jù)中可以看出，的大于的的值, 的小于的的值,在均值相差不大的情況下，方差越小，數(shù)值的離散程度越小,所以第一種預(yù)報方法的準(zhǔn)確度較高，與實測值之間的誤差較小。注：（1）日期豎列中，例如，618表示6月18日； （2）當(dāng)時,，表明此時擬合的最好。5.2 問題二依據(jù)問題二中的分段形式，我們將無雨，小雨，大雨，暴雨，大暴雨，特大暴雨轉(zhuǎn)化為0，1，2，3，4，5，6的數(shù)值形式，同時把問題一中求得的兩種預(yù)報數(shù)據(jù)與實測值依照分級標(biāo)準(zhǔn)轉(zhuǎn)化為數(shù)值形式,預(yù)報等級值與實測值之間的等級范圍為，每一等級差中所包含的站點個數(shù)（n）占總站點數(shù)（N）的比率可以運(yùn)用Matlab求得，如表3所示：表3：等級差與概率-6-5-4-3-2-10123456000000.00550.95500.03850.00090000000000.00580.95360.03940.00090.0001000由于公眾的感受體現(xiàn)在公眾的滿意率與滿意度兩方面，滿意，所以用兩種方法來加以解決5.2.1 模型一通過公眾的滿意率來體現(xiàn)公眾的感受，由于滿意率必須通過統(tǒng)計調(diào)查才能得出，因而我們可以假設(shè)當(dāng)，當(dāng)，當(dāng)，當(dāng)當(dāng)由概率的性質(zhì)和算術(shù)平均數(shù)的計算可知，與之間的關(guān)系：可得第一種方法的，第二種方法的，可知第一種方法比較好。5.2.2 模型二公眾的滿意度也可以體現(xiàn)公眾的感受，公眾根據(jù)滿意程度用評分制方法給其打分，如下表4所示：表4：滿意度與打分非常不滿意比較不滿意不滿意一般滿意比較滿意非常滿意0123456由概率的性質(zhì)和加權(quán)平均數(shù)的計算公式可知，與有如下關(guān)系： 由此可得 第一種方法的，第二種方法的，可知第一種方法比較好。6、模型的改進(jìn)對于問題一，利用二元三次樣條插值法求出了各網(wǎng)格結(jié)點的降雨量預(yù)報數(shù)值，在91個站點經(jīng)緯度位置的數(shù)值，計算結(jié)果精確度高，但是該方法不易理解，不能直觀地反映數(shù)據(jù)的求解過程，我們可以采用“最近鄰點法”較直觀地計算上述的數(shù)據(jù)。最近鄰點法的定義：在圖二中，每個站點附近都分布著數(shù)個網(wǎng)格結(jié)點，由已知的網(wǎng)格結(jié)點的預(yù)報數(shù)據(jù)估計該站點經(jīng)緯度位置的預(yù)報數(shù)據(jù)。公式：，其中表示附近第i個結(jié)點的降雨量值，表示從站點到該i點之間的距離，即緯度差、經(jīng)度差的平方和。具體的做法：取一邊長為0.4的正方形小方格，使小方格的絕對中心與站點的位置重合，在小方格內(nèi)及邊上的結(jié)點均為該站點的最近鄰點，再由上述公式求得各站點經(jīng)緯度的預(yù)報數(shù)據(jù)，然后繼續(xù)問題一求解中的方法，求擬合度并確定預(yù)報方法的準(zhǔn)確性。7、結(jié)果分析與評價7.1 問題一兩種預(yù)報方法分析：第一種預(yù)報方法比第二種預(yù)報方法好，準(zhǔn)確性高。將網(wǎng)格結(jié)點的降雨量利用二元三次樣條插值法算出各觀測站點位置的降雨量，再用擬合優(yōu)度法與實測值進(jìn)行擬合。該模型主要采用了二元三次樣條插值法和擬合優(yōu)度法，將擬合后的數(shù)據(jù)進(jìn)行求平均值和相應(yīng)的方差，反應(yīng)出兩種方法的準(zhǔn)確性。7.2 問題二本模型通過滿意度和滿意率反映的信息不同，運(yùn)用統(tǒng)計學(xué)中算術(shù)平均數(shù)和加權(quán)平均數(shù)建立它們與的函數(shù)關(guān)系，從樣本推算總體，有一定的可行性。參考文獻(xiàn)1孫家 ，樣條函數(shù)與計算幾何，北京：科學(xué)出版社，1982年。2丹尼爾L. 魯賓費(fèi)爾德 羅伯特S. 平狄克，錢小軍譯，計量經(jīng)濟(jì)模型與經(jīng)濟(jì)預(yù)測，北京：機(jī)械工業(yè)出版社。3D.Hanselman,B.Littlefield著，張航 黃攀譯，精通MATLAB6，北京：清華大學(xué)出版社，2002年。4丑紀(jì)范 杜行遠(yuǎn) 郭秉榮，數(shù)值天氣預(yù)報淺談，北京：氣象出版社，1980年。5韓志剛 孫隆和 李林森等，四次卷積插值及其在數(shù)字地圖中的應(yīng)用，兵工學(xué)報，第21卷第4期：311312頁，2000年。附錄：附錄表1文件f020703_1中的數(shù)據(jù)5813832.9833118.51670.00000.000064.40001.62145813933.3000118.85000.00000.000031.10001.99885814133.6667119.26670.00000.00000.00000.64005814333.8000119.80000.00000.80006.20000.33825814633.4833119.81670.000012.20008.50000.39725814733.0333119.03330.00000.900043.70004.12525814833.2333119.30000.000020.800044.000018.43115815033.7667120.25000.00009.50002.90000.12865815433.3833120.15000.00002.60004.300002000120.48330.200054.20000.00000.25195823032.1000118.26670.00000.50000.00000.00005823632.3000118.30000.00000.00001.10000.00005823832.0000118.80000.00000.00000.00000.11305824032.6833119.01670.00008.80007.40000.00005824132.8000119.45000.000036.40001.00000.47275824332.9333119.83330.00006.60000.70000.46995824532.4167119.41670.00000.00000.10000.00005824632.3333119.93330.00000.00000.00000.10895824932.2000120.00000.00000.00000.00000.23715825132.8667120.31670.10006.00000.00000.64365825232.1833119.46670.00000.10000.00000.64735825432.5333120.45000.00000.00000.00000.22265825532.3833120.56670.00000.10000.10000.41215826432.3333121.18330.00001.20000.60000.30085826532.0667121.60000.00000.00000.00000.93735826931.8000121.66670.00000.300012.90003.36125833331.9500118.85000.00000.00000.80000.60005833431.3333118.38330.00000.50000.80000.00005833531.5667118.50000.00000.30000.00000.22585833631.7000118.51670.00000.00000.000000833118.18330.00000.00000.00000.00005834131.9833119.58330.00000.20000.60001.45375834231.7500119.55000.00000.00000.10002.40885834331.7667119.93330.00000.00000.00000.10725834431.9500119.16670.00000.00000.80000.42055834531.4333119.48330.00000.10000.00001.02555834631.3667119.81670.00006.60000.00000.28255834931.2667120.63330.00000.00000.90000.41965835131.8833120.26670.00002.70000.80000.63555835231.6500120.73330.00000.40002.60001.72195835431.5833120.31670.00004.20000.500004167120.95000.00000.00002.20002.12555835831.0667120.43330.00000.00001.50000.42235835931.1500120.63330.00000.00002.70000.79155836031.9000121.20000.00000.10001.30000.57235836131.1000121.36670.00000.000011.70001.45475836231.4000121.48330.00000.00009.90002.72555836531.3667121.25000.00001.700011.80002.37865836631.6167121.45000.30006.700016.00004.33265836731.2000121.43330.00001.20003.40001.85455836931.0500121.78330.00001.600010.70008.22295837031.2333121.53330.00001.40003.80002.58935837731.4667121.10000.00002.10006.50003.76035842630.3000118.13330.00000.00000.00006.01115843130.8500118.31670.00000.00000.00000.23145843230.6833118.40000.00000.00000.00000.26825843330.9333118.75000.00000.00000.10000.90135843530.3000118.53330.00000.40000.10000.41405843630.6167118.98330.00000.00001.70000.44655843830.0833118.58330.00000.00000.00000.32865844130.8833119.41670.00001.20000.70000.29255844231.1333119.18330.10002.60000.00000.31535844330.9833119.88330.00000.00000.10000.51775844630.9667119.68330.60002.60004.60000.51675844830.2333119.70000.50002.80002.20000.42685844930.0500119.95000.00001.200013.70000.26395845030.8500120.08331.70002.60001.10000.83865845130.8500120.90000.00000.500011.90000.48405845230.7833120.73330.10000.100010.50000.42095845330.0000120.63330.00000.40007.40000.65075845430.5333120.06670.10002.20000.50000.38465845530.5167120.68330.00000.20005.80000.11265845630.6333120.53330.00000.200012.00000.42135845730.2333120.16670.00002.900013.50000.26065845930.2000120.31670.00001.400017.20000.42925846030.8833121.16670.00001.100010.30001.37935846131.1333121.11670.00000.000011.20001.35325846231.0000121.25000.00000.30009.40001.05455846330.9333121.48330.00002.10006.40002.99795846430.6167121.08330.00000.30009.70001.52895846730.2667121.21670.40000.60001.40007.21325846830.0667121.15000.00004.70003.40008.01975847230.7333122.45000.00001.000011.000016.43225847730.0333122.10000.00005.800016.500067.97775848430.2500122.18330.00003.00009.100029.55875853029.8667118.43330.00000.00000.00000.29295853129.7167118.28330.00000.00000.00000.00005853429.7833118.18330.00000.00000.00000.00005854229.8167119.68330.00007.400016.60000.65625855029.7000120.25000.00003.90004.20000.38345856229.9667121.75000.00002.100010.100015.2070注：由于41天各個站點的預(yù)測數(shù)據(jù)太大，在此不一一列舉。程序1 二元三次樣條插值求各站點預(yù)測降雨量：format longload lon.dat;load lat.dat;lon1=;lat1=;for i=1:53lon1=lon1,lon(i,:);lat1=lat1,lat(i,:);endfor time=1:4 switch time case 1 load f7031_dis1; loadf=f7031_dis1; case 2 load f7032_dis1; loadf=f7032_dis1; case 3 load f7033_dis1; loadf=f7033_dis1; case 4 l