数学建模竞赛2005年C题论文.doc

雨量预报方法评价的数学模型摘要：降雨是日常生活中最常见的天气现象，及时准确地对降雨量作出预报是一个十分困难的问题。题目给出了用两种降雨量预测方法预测得出的41天各时段网格结点的预报数据，我们对预报方法进行评价。先使用二元三次样条插值法，对91个站点位置进行插值，从而求出各站点预报数据，然后用方差分析计算两种预测方法每一天的拟合优度，从而评价两种预报方法的准确性。第一种方法拟合优度的平均值为0.9087,第二种方法为0.9054,可见这两种方法的预测准确性都比较高，但第一种比第二种的精确度更高。在评价方法中考虑公众的感受，即考虑公众对预报结果的满意率和满意度。本文建立了两个模型，即满意率模型和满意度模型。满意率指满意数与总数的比率；满意度指达到公众满意的程度。两个模型分别运用算术平均数和加权平均数的方法以及概率的知识进行求解，第一种方法的，第二种方法的，可见公众对这两种方法的预测值都比较满意，但第一种方法给公众的感受优于第二种。关键词：二元三次样条插值法 拟合优度 满意率 满意度1、问题的重述本题是一道由气象部门提出，期望建立一种科学评价雨量预报方法好坏的数学模型与方法。在本题中，气象部门所采用了两种方案研究6小时雨量预报方法。所谓6小时雨量预报方法是指每天晚上20点预报从21点开始的4个时段（21点至次日3点，次日3点至9点，9点至15点，15点至21点）在某些位置的雨量，这些位置位于东经120度、北纬32度附近的5347的等距网格点上。同时设立91个观测站点实测这些时段的实际雨量，由于各种条件的限制，站点的设置是不均匀的。气象部门提供了41天的用两种不同方法的预报数据和相应的实测数据。预报数据在文件夹FORECAST中，实测数据在文件夹MEASURING中。FORECAST中的文件lon.dat和lat.dat分别包含网格点的经纬度，其余文件名为_dis1和_dis2，例如f6181_dis1中包含2002年6月18日晚上20点采用第一种方法预报的第一时段数据（其2491个数据为该时段各网格点的雨量），而f6183_dis2中包含2002年6月18日晚上20点采用第二种方法预报的第三时段数据。MEASURING中包含了41个名为.SIX的文件，如020618.SIX表示2002年6月18日晚上21点开始的连续4个时段各站点的实测数据（雨量）。题目中有如下假设：第一，雨量用毫米做单位，小于0.1毫米视为无雨。第二，气象部门将6小时降雨量分为6等：0.12.5毫米为小雨，2.66毫米为中雨，6.112毫米为大雨，12.125毫米为暴雨，25.160毫米为大暴雨，大于60.1毫米为特大暴雨。该题目的问题如下：问题一：请建立数学模型来评价两种6小时雨量预报方法的准确性；问题二：若上述分级向公众预报，如何在评价方法中考虑公众的感受？2、问题的分析2.1 问题一的分析要评价两种预报方法的准确性，关键是需要把题目中给出的大量预报数据进行处理，使其能与各站点实测值进行比较。因此我们分为两个大的步骤处理该问题：第一，根据观测点的分布，将网格结点预报数据转化成观测站点预报数据；第二，根据网格结点预报数据插值成的各站点预报值与实测值比较，判断两种方法的优劣。在该问题中所使用的数学工具和方法有：第一，根据各网格结点的降雨量采用二元三次样条插值法求出各观测站点位置的降雨量；第二，通过拟合优度分析，判断出两种雨量预报方法的准确程度。2.2 问题二的分析在评价方法中考虑公众的感受，即公众对预报结果准确度的满意率和满意度。若要求总体的满意率和满意度，首先应该确定每一等级差下站点个数(n)与总站点个数的比率,通过概率的基本性质与平均数来建立比率与满意度的函数关系，最终得出总体的满意率和满意度。3、基本假设（1） 题目中所给的经纬度范围内的地球面是平面；（2） 题目中给出的每个降雨量数值是相应经纬度点数据，观测站点可以抽象为一个点，其面积大小不影响数据观测；（3） 降雨量小于0.1毫米的视为无雨；4、符号说明 样本的均值； 样本中第点的函数值； 回归直线上第点的函数值； 样本中第点的函数值与回归直线上第点的函数值的差值； 的总方差，即样本中各点的函数值与样本的均值差的平方和； 的剩余变差，即样本中各点的函数值与回归直线上各点函数值差的平方和； 的解释变差，即回归直线上各点函数值与样本均值差的平方和； 拟合优度，即对的回归方程所能解释的那部分变差占总变差的比例 拟合优度的平均值 拟合优度的方差 第一种方法每一天的拟合优度 第二种方法每一天的拟合优度 站点处的预报数据与其实测数据的等级差5、模型的建立与求解5.1 问题一将球面问题转化成平面问题解决。5.1.1 将91个观测站点及2491个网格结点的经纬度联系起来，用Matlab软件编程画图，如下图1所示：“.”表示2491个网格结点的位置，“*”表是91个站点位置；图1 原始网格结点及站点图通过观察发现，91个测量站点的位置都大约集中在经度118123度、纬度29.534度之间的范围内。问题一中要求对两种预报数据的方法进行评价，实际是91个站点位置的预报数据值与相应站点实测值的比较。由此，可选取上述范围作为工作区，将对结果无影响的网格结点排除出考虑范围，即只对在经度118123度、纬度29.534度的范围内进行数据的处理和评价，截取图1，如图2所示：图2 处理后的网格结点与站点图再用二元三次样条插值法求各观测站点的降雨量，具体如下：二维三次样条插值法可以直接利用网格结点的数值进行插值，但在网格结点处不能得到连续的二阶导数，其插值函数式可表示为：其中，为采样函数，为插值函数，（，）为插值点，为的在各点的函数值，为采样间隔，为采样函数的个数，、是插值内核。根据此插值法，我们运用Matlab软件中的二维三次样条函数的命令（interp2（x，y，z，xj，yj，spline）可以插值出91个站点41天每个时段的降雨量。我们把41天中91个站点的各时段预测数据(降雨量)保存在41个名为_1/2文件中，如f020703_1表示使用第一种预报方法的数据插值得到2002年7月03日晚上21点开始的连续4个时段各站点的预测数据。这些文件的数据格式如表1所示：表1文件_1/2的数据格式站号纬度经度第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中的完整数据参见附录表。将各个站点预测数据与实测数据进行比较，观察其插值效果。比如将7月03日晚上21点开始的第四个时段各个站点的预测降雨量与这个时段各个站点的实测降雨量进行比较，插值效果图如图3所示(其中红“*“代表的预测降雨量，绿“o“代表的实测降雨量)：图3 各个站点插值效果图由图3,可以看出预测降雨量与实测降雨量相差比较小，预测降雨量比较符合实际。5.1.2 所谓拟合优度原理，它是关于估计回归直线与数据之间拟合程度的一个很有用的度量，残差大说明拟合效果欠佳，反之，则拟合的好。残差的方差除以Y的偏差平方和可以用来衡量拟合优度（偏差平方和），,为了计算偏差平方和，可以将等式两边同时平方，且对所有观测值求和：用最小二乘残差的两个性质和（为残差），可得： Y的总变差 Y的剩余变差 Y的可解释变差 是的总变差中，对的回归方程所能解释的那一部分变差占总变差的比例，的值大则回归直线拟合的好，反之则不好。针对此题，我们通过回归变差，总变差求得,再求出每种方法下的平均值和方差，如表2所示：表2：拟合优度表日期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从上表数据中可以看出，的大于的的值, 的小于的的值,在均值相差不大的情况下，方差越小，数值的离散程度越小,所以第一种预报方法的准确度较高，与实测值之间的误差较小。注：（1）日期竖列中，例如，618表示6月18日； （2）当时,，表明此时拟合的最好。5.2 问题二依据问题二中的分段形式，我们将无雨，小雨，大雨，暴雨，大暴雨，特大暴雨转化为0，1，2，3，4，5，6的数值形式，同时把问题一中求得的两种预报数据与实测值依照分级标准转化为数值形式,预报等级值与实测值之间的等级范围为，每一等级差中所包含的站点个数（n）占总站点数（N）的比率可以运用Matlab求得，如表3所示：表3：等级差与概率-6-5-4-3-2-10123456000000.00550.95500.03850.00090000000000.00580.95360.03940.00090.0001000由于公众的感受体现在公众的满意率与满意度两方面，满意，所以用两种方法来加以解决5.2.1 模型一通过公众的满意率来体现公众的感受，由于满意率必须通过统计调查才能得出，因而我们可以假设当，当，当，当当由概率的性质和算术平均数的计算可知，与之间的关系：可得第一种方法的，第二种方法的，可知第一种方法比较好。5.2.2 模型二公众的满意度也可以体现公众的感受，公众根据满意程度用评分制方法给其打分，如下表4所示：表4：满意度与打分非常不满意比较不满意不满意一般满意比较满意非常满意0123456由概率的性质和加权平均数的计算公式可知，与有如下关系： 由此可得 第一种方法的，第二种方法的，可知第一种方法比较好。6、模型的改进对于问题一，利用二元三次样条插值法求出了各网格结点的降雨量预报数值，在91个站点经纬度位置的数值，计算结果精确度高，但是该方法不易理解，不能直观地反映数据的求解过程，我们可以采用“最近邻点法”较直观地计算上述的数据。最近邻点法的定义：在图二中，每个站点附近都分布着数个网格结点，由已知的网格结点的预报数据估计该站点经纬度位置的预报数据。公式：，其中表示附近第i个结点的降雨量值，表示从站点到该i点之间的距离，即纬度差、经度差的平方和。具体的做法：取一边长为0.4的正方形小方格，使小方格的绝对中心与站点的位置重合，在小方格内及边上的结点均为该站点的最近邻点，再由上述公式求得各站点经纬度的预报数据，然后继续问题一求解中的方法，求拟合度并确定预报方法的准确性。7、结果分析与评价7.1 问题一两种预报方法分析：第一种预报方法比第二种预报方法好，准确性高。将网格结点的降雨量利用二元三次样条插值法算出各观测站点位置的降雨量，再用拟合优度法与实测值进行拟合。该模型主要采用了二元三次样条插值法和拟合优度法，将拟合后的数据进行求平均值和相应的方差，反应出两种方法的准确性。7.2 问题二本模型通过满意度和满意率反映的信息不同，运用统计学中算术平均数和加权平均数建立它们与的函数关系，从样本推算总体，有一定的可行性。参考文献1孙家 ，样条函数与计算几何，北京：科学出版社，1982年。2丹尼尔L. 鲁宾费尔德 罗伯特S. 平狄克，钱小军译，计量经济模型与经济预测，北京：机械工业出版社。3D.Hanselman,B.Littlefield著，张航 黄攀译，精通MATLAB6，北京：清华大学出版社，2002年。4丑纪范 杜行远 郭秉荣，数值天气预报浅谈，北京：气象出版社，1980年。5韩志刚 孙隆和 李林森等，四次卷积插值及其在数字地图中的应用，兵工学报，第21卷第4期：311312页，2000年。附录：附录表1文件f020703_1中的数据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天各个站点的预测数据太大，在此不一一列举。程序1 二元三次样条插值求各站点预测降雨量：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