数学建模论文基于贵州省部分区域雷达监测降雨量的分析

第一届“九章数学杯”数学建模挑战赛承 诺 书 我们仔细阅读了第一届“九章数学杯”数学建模挑战赛的竞赛规则。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们接受相应处理结果。我们允许海南大学数学建模协会公布论文，以供同学之间学习交流，海南大学数学建模协会网站以非商业目的的论文交流不需要提前取得我们的同意。 参赛队员 (签名) ：队员1：赵荥仁 队员2：李锐 队员3：秦昊 基于贵州省部分区域雷达监测降雨量的分析摘要针对雷达检测系统对贵州省2015年5月7日0时至11时降雨天气的记录，本文将采用插值拟合、积分求和等方式求解总降雨量。模型修正改进分析中利用查阅资料所得结果总结出最佳有效模型。对于问题一，我们首先使用相关绘图软件对雷达回波图中的贵阳市进行剪辑分割，对贵阳市区域面积上各时间段的雷达回波图的降水量进行样点采集和预处理。其次，在不考虑雷达系统出现故障以及每半小时内出现降雨量不均等等问题下使用MATLAB将预处理之后的雷达回波图转化为灰度图，将之转化成矩阵以知道每个像素点的数值，再通过数值统计计算出每种降雨量程度的分布面积。最终根据面积与基本反射率乘积在时间上进行累加，从而计算出贵阳市逐小时降雨量及总降雨量。对于问题二，根据附录表格及雷达监测范围125km以内确定包括贵阳市在内的清镇市、修文县、平坝县、花溪区、长顺县、开阳县、惠水县、龙里县、贵定县、都匀市、麻江县、息烽县等13个区域进行降雨量总量的强度的比较。在问题一模型的基础上运用插值拟合分析，并在整个降雨时间内积分求和以分别计算出各个地区的总降雨量，最终比较13个区域的总降雨量，得出最大总降雨量的区域。对于判断看图识降雨量的方式精确问题，我们先分析引起产生误差的因素，再次，根据题意我们先使用卡尔曼滤波模型计算出的总降雨量为基准参量，再利用雷达回波图计算出的总降雨量与的比值为估算的准确程度。关键词：总降雨量 基本反射率 灰度图 插值拟合 卡尔曼滤波一、 问题重述2015年5月7日，贵州省迎来的一场大范围降雨天气，此次降雨天气过程从00:00开始至11:00结束，整个过程被雷达记录了下来（每半小时一次），基本上反映了这次降雨天气的整个过程。再此基础上计算出贵阳市逐小时降雨量及总降雨量。并计算此次降雨过程中降雨量最大的地方可能会出现在哪个区域；最终对这种看图识降雨量的方式准确程度做出判断估算。二、 问题分析对于问题一，针对贵阳市逐小时降雨量及总降雨量的求解，在雷达回波图能正确反映降雨量的条件下以面雨量大小在时间上的累积进行计算。所以必须求出贵阳市面积以及降雨量情况分布面积。运用MATLAB软件将预处理之后的雷达回波图转化为灰度图，将之转化成矩阵以知道每个像素点的数值，再通过数值统计计算出每种降雨量程度的分布面积，再进行叠加求和。而对于求解贵阳市逐小时降雨量由于给出的雷达回波图是每间隔半小时监测出的数据，因此则还需进一步的进行预测处理。而对问题二而言，贵阳市也处于监测范围内，它也可能是总降雨量最大的区域，因此需在问题一的基础上进行总降雨量最大的区域估计。其次运用插值拟合和积分方法对其余12个区域进行总降雨量求解，最终比较13个区域的总降雨量来估计总降雨量最大的地方。针对问题三，根据题意求解对看图识降雨量的方式准确程度判断问题，应先采用最佳有效模型测算出贴合实际降雨情况的总降雨量，再应用问题一二的模型求解出的总降雨量与进行比值推算，从而作为看图识降雨量的准确程度估算值。三、 符号说明面雨量特定区域面积有限元上雨量各基本反射率在雷达回波图中的像素点个数基本反射率，单位6mm/m3表示不同基本反射率的下标各区域的总降雨量，单位m3表示区域名称表示逐小时间隔段准确程度准确程度的比较基准参量四、 模型假设4.1对问题一二的基本假设（1）雷达监测系统运行正常。（2）在雷达监测每间隔半小时之内的大气对流稳定。（3）在125数据范围内各区域地形地貌相似。（4）在层状云降雨时，不同的降雨空间分辨率相同。（5）地面物质对雷达监测不产生反射影响。4.2对问题三的基本假设（1）雷达监测系统运行正常。（2）在125数据范围内各区域地形地貌相似。五、 模型建立与模型求解5.1贵阳市逐小时降雨量及总降雨量模型求解5.1.1基本反射率下的面积整合使用相关绘图软件对雷达回波图中的贵阳市进行剪辑分割，对贵阳市区域面积上各时间段的雷达回波图的降水量进行样点采集和预处理（部分图片见下图1所示）。 00:24:45 03:03:13 05:22:50 09:02:30图1贵阳部分回波图预处理然后运用MATLAB中H=imread(文件名)语句将预处理之后的雷达回波图转化为灰度图，将之转化成矩阵以知道每个像素点的数值，具体执行程序见附录。再通过数值统计计算出每种降雨量程度的分布面积大小，再进行叠加求和，即像素点个数总和。5.1.2面雨量的定义在分析预报水情变化时面雨量运用相当广泛，一般将其定义为相关流域（或区域）各个雨量站测得的降水量（点雨量）计算出流域（或区域）平均降水量，这种降水量称为面雨量。在水文学中，将面雨量表示为： （1）式中为面雨量，为特定区域面积，为有限元上雨量1。区域面雨量能客观的反映该区域降雨情况。5.1.3逐小时降雨量与总降雨量求解根据预处理之后的雷达回波图和MATLAB的数值转换分别统计出每种降雨量程度的像素点个数总和，设基本反射率为，降雨量为立方米，假设像素点近似等于单位面积，则有关系式：（2）在计算贵阳市逐小时降雨量时，题目给出的附件中雷达回波图是以每半个小时为间隔进行监测测量，而待解决的问题是要求以逐小时为单位时间间隔求出每个小时内的降雨量，因此需要进行时间段修正，将每半小时之内检测一次的降雨量延长至每一个小时再进行降雨量累积。公式（2）变为： （3）此时对应各时间段进行的像素点数乘以基本反射率的积分结果就将为贵阳市逐小时降雨量。由以上式子运用MATLAB计算出贵阳市逐小时降雨量具体数值结果见表格（2）。表格（1）基本反射率对应降水量分布表dBZ颜色降雨量45-50红色20mm/h40-45土黄色15mm/h30-40亮黄色10mm/h25-30绿色5mm/h20-25浅绿色3mm/h10-20浅蓝色1mm/h0-10蓝色0mm/h表格（2） 贵阳市逐小时降雨量分布表记录时刻有效面积像素值红色像素值土黄色像素值亮黄色像素值绿色像素值浅绿色像素值浅蓝色像素值蓝色像素值降雨量（mm）5.6；23:45:228940742723780.05525.7；00:24:458940129558350.01465.7；01:09:47894017954717990.16135.7；01:37:548940721792660.06035.7；02:00:2489405831342720.21065.7；02:28:298940240511724599359917437223.94185.7；03:02:1389408112001368139518399723565.18585.7；03:30:18894097810351933115217378915897.41355.7；04:49:0389401841470781204911841233.31275.7；05:50:59894020291771912456852.37415.7；06:52:528940103112961861260972.50255.7；07:26:408940631180711236823761.66625.7；08:23:03894078260110170.18675.7；09:02:30894046105117500.09985.7；10:27:01894023210.00265.7；11:00:5089400总降雨量：21.63。注：降雨量=（相应颜色像素值分别/有效面积像素值）*相应颜色的降雨量 每个时刻（逐时刻）降雨量已用红体字标出，总降雨量即逐时刻降雨量之和。对表格（2）贵阳市逐小时降雨量分布表进行绘图整理如图2所示：图2 贵阳市逐小时降雨量分布图5.2最大降雨量区域评估依据题目附录表格中给出的数据，可以假定雷达监测的125km范围内的清镇市、修文县、平坝县、花溪区、长顺县、开阳县、惠水县、龙里县、贵定县、都匀市、麻江县、息烽县12个区域做为总降雨量最大的区域监测评估。由于贵阳市完全处在雷达监测的范围之类，而且在求解问题一的基础上已经求出贵阳市的总降雨量强度较大，由此可推断贵阳市也有可能成为总降雨量最大的区域。因此为减小模型的求解误差，我们将贵阳市也列如总降雨量最大的区域监测评估之中。至此将有13个地区参与评估。在求解贵阳市逐小时降雨量及总降雨量的过程中，已经假设雷达回波图能准确反映降雨实际降雨量，因此可以在第一问贵阳市逐小时降雨量及总降雨量的解答基础上采用插值拟合方法拟合计算出各个区域降雨量的拟值函数方程，再根据各个区域的预处理之后的雷达回波图和MATLAB的数值转换分别统计出每种降雨量程度的像素点个数总和、基本反射率和降雨量三者的关系，再在整个监测过程中进行时间上的积分求和，列出下列关系式： （4）由公式（4）可计算出各个地区的总降雨量。根据公式（4）和附录表格中的数据用MATLAB插值拟合分别得出各个区域的函数表达式、拟合图形以及进行积分求和之后的总降雨量：清镇市降雨量：；当t=10，得y=-0.2963当y=0，得t=9.3420s所以降雨量s=0.5*（9.3420*4.2067-0.6580*0.2963）=19.5520mm修文县降雨量：当t=10，y=0.024故s=0.5*（2.0600+0.024）*10=10.4200mm平坝县降雨量：当t=10s，y=-0.6713当y=0，t=8.91s故s=0.5*(5.4867*8.91-1.09*0.6713)=24.0774mm花溪区降雨量：当t=10，y=0.7217，故s=0.5*（0.7217+4.0067）*10=23.6420mm长顺县降雨量： 当t=10，y=2.5747，故s=0.5*（7.1867+2.5747）*10=48.807mm开阳县降雨量：当t=10，y=0.0417故s=0.5*（0.3267+0.0417）*10=1.8420mm惠水县降雨量： 当t=10，y=1.6837所以s=0.5*（1.6837+8.1867）*10=49.352mm龙里县降雨量：当t=10，y=1.4所以s=0.5*（1.4+2.8）*10=21mm贵定县降雨量：当t=10，y=0.2483所以s=0.5*（0.2483+0.7333）*10=4.908mm都均市降雨量：当t=10，y=-3.3193当y=0，t=0.7950所以s=-0.5*（0.2867*0.7950-9,205*3.3193）=15.1631mm麻江县降雨量：当t=10，y=-3.6697当y=0，t=0.09所以s=-0.5（0.09*0.0333-9.91*3.6697）=18.1819mm息烽县降雨量：当t=10，y=-0.02当y=0，t=9所以s=0.5*（9*0.18-1*0.02）=0.8mm贵阳市降雨量：当t=10，y=1.2714所以s=0.5*（2.5794+1.2714）*10=19.254mm综上所述，进行最终总降雨量比较得降雨量最大的地区是惠水县。在雷达监测图中的位置用黑灰色标示出，如图（3）所示：图（3） 总降雨量最大区域图5.3看图识降雨量的方式准确程度判断根据Pessoa等2用不同空间分辨率的雷达降雨研究结论，不同的降雨空间分辨率对于模拟的精度有很大的影响，降雨积累的时间越短，结果的误差越大，相对长时间的降雨累计可以消除由于短时间降雨累计的误差，然而降雨的空间分辨率并不是越小越好。采用多大的时间和空间分辨率更合理，与所选的模型、流域特征、降雨的类型、降雨的强度以及流域的尺度等因素有关，并且大气对流程度以及水分蒸发对雷达监测的基本反射率都存在一定影响，其次，看图识降雨量主要是依据雷达监测的基本反射率来判断，方式简单直观，但基本反射率对降水粒子直径有所要求，且影响雷达监测的因素较多，因此看图识降雨量的方式必定存在一定误差。查阅资料可知采用卡尔曼滤波算法得出的监测结果更加准确，更加契合实际降雨量情况，因此我们近似将卡尔曼滤波算法得出的总降雨量作为准确精度的比较基准参量，设看图识降雨量的方式计算出的总降雨量为，准确程度为，则： （5）由公式（5）可算出看图识降雨量的方式的准确程度。六、 模型结果分析与模型优缺点在本次模型建立中，对于问题一贵阳市逐小时降雨量和总降雨量的求解，模型相对简单，计算比较容易也容易理解，利用降雨量在时间和空间上的分段累积求出的结果大体接近真实值，但仍然存在一定误差，不能真正准确的计算出实际降雨量值。而对于问题二的求解过程中，本模型是基于问题一的模型上进行进一步的优化改进，即插值拟合与积分累加的联合求解，使总降水量的计算过程更加方便简洁。不足之处就是在于模型的建立必须依赖假设条件的多样性，模型成立的要求条件较多，例如要求地形相似，忽略大气对流及水汽蒸发等，然而这些因素在实际情况中是仍然存在的，所以模型求解出来的误差相对较大。而针对问题三，由于雷达监测站网不够多，数据量小，采用卡尔曼滤波算法不能完整求出精确的总降雨量，因此看图识降雨量的方式的准确程度只能是一个估量值。七、 改进方向以上模型结果分析与模型优缺点可清楚的知道，此次模型的建立条件假设约束太多与实际上的总降雨量存在较大的出入，所以要想减小模型结果误差需要对模型进行优化改进。在模型求解面雨量及总降雨量过程中，不仅仅只用MATLAB将图片转化为灰度图，再根据像素点总数乘以基本反射率，也应结合附录表格中给出的对贵州省12个区域的监测数据运用Voronoi图3进行区域降水面积求解，再根据相应的权重比比较两种方法哪一种测出的总降雨量与实际降雨相近，误差更小，更合理，由此确定最佳有效模型。而看图识降雨量方式的准确程度判断应在最佳有效模型求解的基础上进行影响检测因素的权重分析，再按影响检测结果的因素所占权重比对看图识降雨量方式及最佳有效模型进行评估对比，进一步估算出准确程度。八、 参考文献1葛徽衍，张永红.流域面雨量的一种估算方法-网络插值法J.陕西气象，2006，（1）：34-35.2Pessoa M L，Bras R L，Williams E RUse of weather radar for flood forecasting in the sieve river basin：A sensitivity analysisJJApplMeteor，1993,32：462-4753司守奎，孙玺菁，张徳存，周刚，韩庆龙，数学建模算法与应用北京：国防工业出版社，2013.1九、 附录1.1雷达回波图转化为灰度图H=imread(文件名)1.2贵阳市逐小时降雨量x=1:1:12;y=0.05520.01460.06033.94187.41353.31272.37412.50251.66620.18670.09980.0026;plot(t,y)title(贵阳市逐小时降雨量)xlabel(时间)ylabel(降雨量)2.1第二问代码 x1=1:1:10;y1=0.1 6.1 5.4 1.1 2.5 2.1 0 0 0 0;plot(x1,y1,*) xlabel(时间) ylabel(降雨量) title(清镇市逐小时降雨量) a=polyfit(x1,y1,1)a = -0.4503 4.2067 x2=1:1:10; y2=0 0.9 5.6 0.5 0.1 2.3 0 0 0 0; plot(x2,y2,*) xlabel(时间) ylabel(降雨量) title(修文县) a=polyfit(x2,y2,1)a = -0.2036 2.0600 x3=1:1:10; y3=0 6.2 10.8 0.4 3.3 0.3 0 0 0 0; plot(x3,y3,*) xlabel(时间) ylabel(降雨量) title(平坝县) a=polyfit(x3,y3,1)a = -0.6158 5.4867 x4=1:1:10;y4=0 0 4.4 12.3 0.8 4 0.5 0 0 0; plot(x4,y4,*)xlabel(时间)ylabel(降雨量)title(花溪区)a=polyfit(x4,y4,1)a = -0.3285 4.0067 x5=1:1:10;y5=0 0 1 28.8 0.7 16 0 0 0 0; plot(x5,y5,*)xlabel(时间)ylabel(降雨量)title(长顺县)a=polyfit(x5,y5,1)a = -0.4612 7.1867 x6=1:1:10;y6=0 0 0.9 0.4 0.1 0 0.2 0.1 0 0; plot(x6,y6,*)xlabel(时间)ylabel(降雨量)title(开阳县)a=polyfit(x6,y6,1)a = -0.0285 0.3267 x7=1:1:10;y7=0 0 0 34.8 7.8 2.8 0.7 0 0 0; plot(x7,y7,*)xlabel(时间)ylabel(降雨量)title(惠水县)a=polyfit(x7,y7,1)a = -0.6503 8.1867 x8=1:1:10;y8=0 0 0 8.1 8.2 1.4 2.5 0.1 0 0; plot(x8,y8,*)xlabel(时间)ylabel(降雨量)title(龙里县)a=polyfit(x8,y8,1)a = -0.1400 2.8000 x9=1:1:10;y9=0 0 0 0.5 4 2.8 1.5 1.1 0.1 0; plot(x9,y9,*)xlabel(时间)ylabel(降雨量)title(贵定县)a=polyfit(x9,y9,1)a = 0.0485 0.7333 x10=1:1:10;y10=0 0 0 0 0.9 10.3 6.5 2.5 2.2 0.3; plot(x10,y10,*)xlabel(时间)ylabel(降雨量)title(都均市)a=polyfit(x10,y10,1)a = 0.3606 0.2867 x11=1:1:10;y11=0 0 0 0 1.5 3.8 9.7 5.1 0.6 0; plot(x11,y11,*)xlabel(时间)ylabel(降雨量)title(麻江县)a=polyfit(x11,y11,1)a = 0.3703 0.0333 x12=1:1:10;y12=0 0 0.6 0.1 0 0 0 0 0 0; plot(x12,y12,*)xlabel(时间)ylabel(降雨量)title(息烽县)a=polyfit(x12,y12,1)a = -0.0200 0.1800 x13=1:1:13;y13=0.0552 0.0146 0.0603 3.9418 7.4135 3.3127 2.3741 2.5025 1.6662 0.1867 0.0998 0.0026 0;plot(x13,y13,*)xlabel(时间)ylabel(降雨量)title(贵阳市)a=polyfit(x13,y13,1)a = -0.1308 2.5794