垃圾焚烧厂周边环境的建模论文

1、 承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们的题目是： 垃圾焚烧厂的经济补偿问题 我们参赛年级是（一年级，二年级以上）： 二年级 所属学院（请填写完整的全名，可填多个）：

2、 数学科学学院、英才实验学院、通信与信息工程学院 参赛队员 1. 王吉 2. 张博文 3. 樊聪丽 指导教师或指导教师组负责人 ： 张勇 日期： 2014 年 7 月 16日垃圾焚烧厂污染物高斯烟羽扩散动态监测方案摘要本文综合分析了计划建立的垃圾焚烧厂的污染物对周边环境的影响情况，建立了垃圾焚烧厂周边环境污染监测指标，并为垃圾焚烧厂周边的居民设计了风险承担经济补偿方案。假设该厂焚烧炉排放的污染物符合国家污染物排放新标准的情况下，基于高斯烟羽扩散理论科学定量地分析了周边环境中的污染物的浓度情况，全方位考虑垃圾焚烧厂的垃圾处理量，排烟筒高度，厂址周边的地形地貌，周边的居民点以及包括风速风向在内的气

3、象条件等因素，建立了垃圾焚烧厂周边的环境监测指标，实现了对垃圾焚烧厂烟气排放及相关环境影响状况的动态监控。以实际的监测结果为依据，计算出了空气质量指标AQI的等级，评估出了垃圾焚烧厂对其周边环境所带来的潜在环境污染风险，结合当地经济发展状况以及垃圾焚烧厂的支付意愿等因素为厂址周围承受环境污染风险的居民制订了合理的经济补偿方案。考虑故障发生概率情况下，利用泊松分布理论对故障发生的概率进行估计，计算故障发生情况下污染物产生量的变化，并对以上建立的模型进行了修正，建立更合理的环境污染动态监测指标和居民风险承担经济补偿方案。（还没有建立具体的监测方法，需要在文章所有内容完善后在进行完善摘要的工作）关键

4、词 高斯烟羽扩散，环境动态监测，AQI，经济补偿方案，泊松分布，容斥原理一、问题重述“垃圾围城”问题是世界性难题，在我国这个问题尤为突出，目前我国三分之二的城市面临“垃圾围城”的困境，垃圾焚烧成为我国解决垃圾问题的主要手段之一。但是解决垃圾处理问题的同时垃圾焚烧厂产生的污染物又带来了垃圾焚烧厂周边的环境污染问题。深圳市某地点计划建立一个日处理垃圾量1950吨，拥有3台可处理垃圾650吨/天的焚烧炉的中型垃圾焚烧厂，垃圾焚烧产生的污染物会对周边的环境产生污染，现需要我们解决以下几个问题：1. 基于目前主要是在焚烧厂内进行测量监控，缺少从周边环境出发的外围监控体系，需要综合分析垃圾焚烧厂的污染物对

5、周边环境的影响情况，全方位考虑焚烧厂垃圾处理量，排烟筒高度，厂址周边的地形地貌，包括风速风向在内的气象条件，距离居民点的距离等因素，在假设焚烧炉排放的污染物符合国家标准的情况下，利用科学定量分析的方法建立垃圾焚烧厂周边的环境监测指标模型，实现对环境污染的全方位的动态监测，2. 评估出垃圾焚烧厂对其周边环境所带来的潜在环境污染风险，并为厂址周围的居民设计出合理的经济补偿方案。垃圾焚烧厂在实际运转中，3. 由于各种因素会导致如焚烧炉除尘装置损坏等故障的情况致使污染物的排放量增加，在考虑到故障发生概率的情况下对以上模型进行合理修正。二、问题分析 问题一:假设在焚烧炉排放的污染物符合国家污染物排放新标

6、准的情况下，建立垃圾焚烧厂周边的环境监测指标的问题。 由于垃圾焚烧厂产生的污染物主要为废气，建立环境中污染物的检测指标，需要根据污染物的扩散方式，采用高斯烟羽扩散理论分析污染物在周边环境中的扩散情况，定量计算出在具体的地点空气中污染物的浓度大小，由污染物的浓度大小可以得到AQI（空气质量指数），以AQI作为评估该地点的空气污染情况的具体指标。问题二:设计合理的居民风险承担经济补偿方案 以第一个模型中计算出的AQI级别为参考指标，另外需要考虑要进行经济补偿需要综合考虑接受补偿的个体类别是家庭还是个人，对常住人口还是包括流动人口在内的在一定时期内所有遭受环境污染风险的所有人员进行补偿。补偿是一次性

7、的还是每一段时间进行一次的，以及补偿方式等因素。另外还需要综合考虑当地经济发展水平，以及垃圾焚烧厂的支付意愿。问题三: 考虑到故障发生概率的情况下对模型进行修正只考虑影响污染物产生量的故障，故障的具体类型可能有很多种类，根据题目中所给出的附件二的在线监测数据，将故障抽象为一种类型，以焚烧炉故障作为代表。要建立的焚烧厂有三台垃圾焚烧炉，所以发生主要故障的数目也要考虑在内。对于机器故障概率的分析可以用泊松分布理论。计算故障发生概率下污染物产生量的变化，并且利用所得结果对环境污染动态监测模型和居民风险承担经济补偿模型进行合理修正。三、模型假设1. 在不同年份中同一月的气候状况基本相同。2. 焚烧炉污

8、染物排放在不考虑地形地貌因素、降水因素的条件下，气体在空气中的扩散规律符合高斯烟羽扩散模型。3. 排放的污染物中的氮氧化合物可以全部转化为,且有转化比例【1】。4. 不考虑建筑物对污染物扩散的影响。5. 机器故障发生概率符合泊松分布。6. 以8点的天气计算出来的污染物浓度代替全天整体的污染物浓度。7. 日处理垃圾量为650吨/天的焚烧炉的排烟筒内径为5米，日处理垃圾量为350吨/天的焚烧炉的排烟筒直径为3米。并且焚烧炉的垃圾处理量与污染物排放量成正比，各台焚烧炉之间工作相互独立。四、变量说明变量符号含义单位有效源高度为H，坐标为的点的污染物浓度源强，单位时间内的排放污染物质量平均风速水平、垂直

9、方向的标准差，即y, x方向的扩散参数空间中的点到极点的距离极轴逆时针旋转到达该点的角度radIAQI空气质量分指数无AQI空气质量指数无五、模型建立与求解5.1建模前的准备工作根据高斯烟羽扩散理论计算垃圾焚烧厂周边空间中各点的污染物浓度5.11：空间直角坐标系中某点污染物浓度的计算表达式 计算空间中坐标为的点的污染物浓度:建立以风向为轴正方向的右手螺旋空间直角坐标系，根据高斯烟羽模扩散理论则可以得到该点的污染物浓度为： (1)其中表示有效源高度为H的点的污染物浓度，单位：，分别表示颗粒物，三种污染物；表示第种污染物的源强，即是单位时间内的排放污染物质量，单位：，表示平均风速，单位：，分别为水

10、平、垂直方向的标准差，即y,x方向的扩散参数，单位：。5.1.2：处的点的污染物浓度在本文中为了分析的方便只考虑污染物在地表处的浓度大小，将(1)式中的，得到： (2)5.1.3：将直角坐标系中的污染物浓度表达式转化为在柱坐标系中的污染物浓度表达式由于风向的不断变化，在柱坐标系中计算某点污染物浓度更加简单明了。以风向为极轴正方向，以水平面为极坐标平面，污染物源点在的极坐标平面下在的投影点为坐标原点，建立柱坐标系，则直角坐标系空间中的点就可以表示为，将其代人(2)可得： (3)r表示点到极点的距离，单位：；表示在极轴逆时针旋转到达该点的角度，单位：rad 。5.1.4:污染物浓度表达式中的参数计

11、算1)计算有效源高度H，表示出烟口高度，单位m;表示烟流抬升高度，单位。根据霍兰德公式【2】 (4)式中：表示烟流抬升高度，单位；烟流出口速度，单位；烟流出口内径，单位，根据假设7可知；表示平均风速，单位：；表示烟气出口温度，单位；表示环境平均温度，单位；表示烟囱的热排效率，单位。2)计算扩散参数扩散参数是表示扩散范围及速率大小的特征量，也即正态分布函数的标准差。在计算中采取经验公式【3】得到： (5)称为扩散系数。计算太阳倾角：式中:表示太阳倾角，单位deg；表示一年中日期序数0，1,2,3365，从1月1日开始算起；,单位deg；计算太阳高度角: 式中表示太阳高度角，单位deg；表示当地纬

12、度，单位deg；表示北京时间，单位h;表示当地经度，单位deg。通过查询太阳幅度等级与云量（附录一）可以得到太阳辐射等级，查询太阳辐射等级与地面风速表（附录三）查出太阳稳定度等级，查询大气稳定度和下风距离（附录四、五）得到了扩散参数，将其带入(5)即可求得。3)源强Q的计算。根据量纲分析法知，其中表示第种污染物在出口时的浓度，v表示出烟口的速率，s出烟口内圆的面积。，可以推出 (6)表示出烟口内径，单位m,在计算中。将(3)(4)(5)(6)联立可得： (7)4） 坐标系的统一化由于有的地点存在风向的不同，建立的圆柱坐标系也不同，为了便于观察，将所有源强统一在同一个标准的坐标系中。确定一个标准

13、的柱坐标系，以西风确定的坐标系为标准，也就是极轴是一条向东的射线，表示西风、西南风、南风西北风所建立的坐标系下的极角，对于在不同坐标系下的点，转化到在标准坐标系下的点为:。5.2垃圾焚烧厂周边环境污染全方位动态监控体系模型以空气质量指数AQI作为环境污染监测指标5.2.1:计算垃圾焚烧厂周边空间中z=0处各坐标点的AQI1)IAQI的计算IAQI表示空气质量分指数，代表单项污染物的空气质量指数。污染物种类i的空气质量分指数： (8)式中表示在标准圆柱坐标下点处第i种(分别表示颗粒物，)污染物的空气质量分指数；表示在标准圆柱坐标系下点第i种污染物的浓度；,表示在表一中与相近的污染物浓度限值的高位

14、值；表示在表一中与相近的的污染物浓度限值的低位值；表示在表一中与对应的空气质量分指数；表示在表一中与对应的空气质量分指数。表1：空气质量分指数及其对应的污染物项目浓度限值空气质量分指数（IQAI）污染物项目浓度限值二氧化硫，24小时平均/二氧化氮，24小时平均/颗粒物（粒径小于10），24小时平均/0000505040501001508015015047518025020080028035030016005654204021007505005002620940600根据假设3可知，即二氧化氮的浓度=氮氧化合物浓度除以0.9，于是这里的表示的是转化后的浓度。根据（6）式可得： (9)式中表示空气

15、质量指数。表二中给出了空气质量指数的级别划分及其代表的空气质量的信息与人体健康的关系表2：空气质量指数及相关信息AOIAQI级别AQI类别对健康影响状况050一级优空气质量令人满意，基本无空气污染51100二级良空气质量可接受，但某些污染物开可能对极少数异常敏感人群健康有较弱影响101150三级轻度污染易感人群症状有轻度加剧，健康人群出现刺激症状151200四级中度污染进一步加剧易感人群症状，可能对健康人群心脏、呼吸系统有影响201300五级重度污染心脏病和肺病患者症状显著加剧，运动耐受力降低，健康人群普遍出现症状>300六级严重污染健康人群运动耐受力降低，有明显强烈症状，提前出现某些疾

16、病根据计算出的AQI，查表划分等级并在坐标系中绘制出垃圾焚烧厂周边环境中的AQI等级高于二级的AQI等高线，以等高线的动态变化作为环境污染的动态监控指标。附焚烧炉污染物的排放符合国家污染物排放新标准的情况下AQI的等高线图形。分析曲线中的点，通过观察AQI变化较为明显的地方确立检测点。5.3居民经济补偿方案 根据上面所建立的垃圾焚烧厂周边环境污染全方位动态监控体系模型中所计算出的AQI指数,将AQI级别在二级以上的区域规定为红线区域（即为遭受污染的区域），对于区域内的常住居民（非临时居住的流动人员）以家庭为单位，按照AQI等级的不同分别进行赔偿。赔偿与家庭中人员个数，家庭遭受不同污染程度的土地

17、面积和不同污染级别的赔偿金额有关。赔偿模型的建立如下：对于有n个人的家庭，假设该家庭处于AQI为m级（）的土地面积为，则他应该得到的赔偿金额： 式中表示居住在AQI等级为m级的区域内的人均赔偿金额，单位元/人；表示家庭所拥有的土地位于AQI等级为m级区域内的土地面积，单位；表示对AQI等级为m级的土地的赔偿金额，单位：。以及不是一个恒定的常数，它是由当地的经济发展状况以及垃圾焚烧厂的支付意愿决定的。另外，可以进行除了直接金钱补偿外的基础设施方面的建设进行补偿。赔偿金额的发放每年一次。5.4考虑在故障发生率下的模型的修正5.4.1：附件二中的数据处理：分别计算出颗粒物、折算后的浓度超过国家污染物

18、排放新标准的天数所构成的集合为：。污染物的平均超标浓度分别为，单位；平均合格浓度为，单位。如果为空集，则。即全为平均合格浓度。5.4.2考虑故障概率的情况下的修正模型(1)利用附件二中的数据求解垃圾焚烧炉发生机器故障的概率：对于日垃圾处理量为350吨/天的焚烧炉，取污染物排放超标的天数，如果有一种污染物的排放超过标准排放就认为仪器发生了故障，依据容斥原理可知，泊松分布的期望，一台机器在一天内发生故障服从参数为的泊松分布。则一天内发生k次故障的概率，取得到一台焚烧炉在一天内发生故障的概率 （10）垃圾处理量为650吨/天的一天内的故障发生概率也近似为，依据假设7，三台垃圾焚烧炉是否出现故障相互独

19、立，则三台焚烧炉中发生故障的台数服从二项分布，记为： （11）(2)垃圾焚烧量为650吨/天的焚烧炉出烟口的污染物排放浓度与垃圾处理量为350吨/天的焚烧炉的污染物排放浓度近似相等，由于只考虑一天内的浓度变化，用一个月的每种污染物的平均排放浓度代替一天的排放浓度，即为：颗粒物，的平均超标浓度分别为，单位；平均合格浓度为，单位；污染物排放浓度的变化引起了源强Q的变化， （12）表示第种污染物有台仪器发生故障的前提下的源强，单位；D表示出烟口的内径，依据假设为；表示有台仪器发生故障，单位台。则在点的浓度为： (13)且依据离散型随机变量函数的分布率有：表3：第i种污染物在空间中点浓度的概率分布律用

20、数学期望来代替点的浓度: (14)依据公式（7）（10）（11）（12）（13）(14)可以得到：表示在点第i种污染物修正后的浓度，单位，表示出烟口内径。机器发生故障只会影响到污染物的浓度，所以环境动态监测模型和居民经济补偿模型只需改变污染物的源强。六、结果分析与模型检验七、模型优缺点分析与模型推广7.1模型优点 7.1.1计算污染物在各地点的扩散情况和浓度分布中采取了高斯羽化扩散模型，充分考虑了包括污染物排放量，周围环境等实际情况，计算的结果可信度较高。 7.1.2以国家空气质量指数AQI作为环境污染动态监测指标，具有可行性。7.2模型缺点 7.2.1高斯烟羽扩散理论是建立在地势平坦的基础上

21、的，而在实际情况中，地势有变化，会影响污染物的扩散情况，在本文的模型中未能准确考虑，所以模型会与实际情况有一定偏差。 7.2.2选取的动态环境污染监测指标中，以地面（z=0）处的污染物浓度作为建立指标的一个参数，这个参数的选取存在不合理之处。并且由于气象资料的限制对模型参数的选取并不准确。7.3模型推广研究垃圾焚烧厂的污染问题对解决“垃圾围城”这个世界难题有很重要的实际意义，这个模型不仅可以用于垃圾焚烧厂产生的污染物问题还可以进一步推广到化工厂，电厂等产生的有害废气问题中。如果对污染物的扩散方式进行修改还可以用于液体污染物的分析中。在环保呼声日高的今天，对这些问题的研究意义还是很大的。八、参考

22、文献【1】环境影响评价技术导则【HJ2.2-2008】【2】【3】制定地方大气污染物排放标准的技术方法【GB3840-91】【4】徐全智，杨晋浩 数学建模（第二版）M，北京：高等教育出版社，200806【5】魏楚，沈满洪 基于污染权角度的流域生态补偿模型及应用J，中国人口资源环境，2011,21（6）：137-14210、 附录附录清单：附录一：云量与太阳辐射等级附录二：2013年3月到2014年4月宝安区日云量附录三：大气稳定度的等级附录四：横向扩散参数附录五：垂直扩散参数附录六：模型求解的程序附录一：表格1：云量与太阳辐射等级总云量/低云量-1+1+2+30+1+2+300-1+1000+

23、10000附录二： 表格2：2013年3月到2014年4月宝安区日云量 单位：年月/日23456789101112131415161718192021222324252627282930312011年3月101374074811032922681111222811232011年4月415526351243548777253712397462011年5月5664251600565575218257561503002011年6月605504378765865656867678736572011年7月7446334846765756456782217845582011年8月8055025957806

24、835604575500506002011年9月9646214103797495647563476835782011年10月1083770516710641000000000313716322011年11月1728797253134107864430006100702011年12月21282803422002003040032138310322012年1月3346876402357372201078785315232012年2月48558665534884522356876874852012年3月5164253657588532212323761304072012年4月637374702367

25、46767665646554466附录三： 表格3:大气稳定度的等级地面风速太阳辐射等级+3+2+10-1-2<1.9AABBDEF22.9ABBCDEF34.9BBCCDDE55.9CCDDDDD>6DDDDDD附录四：表格4：横向扩散参数幂函数表达式系数值稳定度下风距离，mA0.9010740.8509340.4258090.60205201000>1000B0.9143700.8650140.2818460.39635301000>1000BC0.9193250.8750860.2295000.31423801000>1000C0.9242790.88515

26、70.1771540.23212311000>1000CD0.9268490.8869400.1439400gt;1000D0.9294180.8887230.1107260gt;1000DE0.9251180.8927940.09856310.12430811000>1000E0.9208180.8968640.08640010.10194711000>1000F0.9294180.8887230.05536340.073334801000>1000附录五：表格5：垂直扩散参数幂函数表达式系数值稳定度下风距离，mA1

27、.121511.513602.10881007999040.008517710.0002115450300300500>500B0.9614351.093560.1271900.0570250500>500BC0.9110151.007700.1146820.07571820500>500C0.9175950.106803>0CD0.8386280.7564100.8155750.1261520.2356670.136659O2000200010000>10000D0.8262120.6320230.555360.1046340.4001670.810763110

28、00100010000>10000DE0.7768640.5723470.4491490.1117710.52899221.0381002000200010000>10000E0.7883700.5651880.4147430.09275290.4333841.7324111000100010000>10000F0.7844000.5259690.3226590.06207650.3700152.4069111000100010000>10000附录五：模型求解的程序（1）% 第一问的高斯模型计算% zbw 2014.7.14编写% 2014.7.14 最后修改% 第一

29、问以报表中一个月的流量平均值作为典型流量值，通过国家标准计算排放量为QNQflu=832.05*(104); % 日流量QN1=20,80,250*Qflu; % 日排量 单位mg/dQN1=QN1/(24*3600); % 单位统一为mg/s% QN 第4列颗粒物，第7列so2，第10列noxQN2=1,26.97,24.87,0.16,33.22,30.63,0.19,53.86,49.65,0.31,883.83,8.35,95.23,0,0,02,19.61,29.94,0.09,42.38,64.72,0.19,45.93,70.13,0.2,648.28,13.36,82.11,0

30、,0,03,19.13,56.4,0.08,42.39,124.98,0.18,43.65,128.69,0.19,633.07,17.04,81.26,0,0,04,19.36,35.44,0.09,42.52,77.84,0.19,44.75,81.92,0.2,640.38,14.63,81.67,0,0,05,20.5,18.59,0.1,41.64,37.76,0.19,50.19,45.51,0.23,676.82,8.13,83.7,0,0,06,20.93,19.52,0.1,42.07,39.25,0.2,49.29,45.98,0.23,690.49,8.5,84.46,0

31、,0,07,20.99,17.93,0.1,41.18,35.18,0.2,44.65,38.15,0.21,692.6,7.34,84.58,0,0,08,21,18,0.1,40.43,34.66,0.19,44.7,38.32,0.21,693.01,7.39,84.6,0,0,09,20.93,16.94,0.1,38.43,31.1,0.18,44.41,35.95,0.21,690.71,6.59,84.47,0,0,010,20.9,18.53,0.1,37.5,33.24,0.18,44.27,39.25,0.21,689.58,7.84,84.41,0,0,011,21.6,

32、20.1,0.11,38.28,35.6,0.19,47.14,43.85,0.23,712.2,8.46,85.67,0,0,012,22.52,20.56,0.11,40.8,37.23,0.2,50.88,46.44,0.25,741.62,8.22,87.31,0,0,013,23.41,20.9,0.12,41.89,37.39,0.21,54.47,48.63,0.28,769.9,7.93,88.88,0,0,014,24.01,21.39,0.12,39.92,35.58,0.21,56.9,50.71,0.29,789,7.91,89.95,0,0,015,24.47,23.

33、16,0.13,39.01,36.93,0.2,58.77,55.63,0.31,803.75,8.68,90.77,0,0,016,25.13,24.75,0.13,39.74,39.14,0.21,61.46,60.54,0.33,824.91,9.16,91.95,0,0,017,26.77,25.8,0.15,43.3,41.72,0.24,57.43,55.34,0.32,877.41,8.89,94.87,0,0,018,28.6,27.54,0.17,49.33,47.5,0.29,58.14,55.98,0.35,936.04,8.88,98.14,0,0,019,29.74,

34、27.66,0.18,51.73,48.1,0.32,61.7,57.38,0.38,972.49,8.45,100.17,0,0,020,29.86,30.31,0.18,50.6,51.37,0.31,62.05,63,0.38,976.06,9.51,100.37,0,0,021,27.64,29.21,0.16,38.04,40.19,0.22,55.14,58.25,0.32,905.33,9.96,96.43,0,0,022,27.21,28.2,0.16,35.59,36.88,0.2,53.79,55.73,0.31,891.52,9.74,95.66,0,0,023,28.6

35、8,30.09,0.17,40.13,42.11,0.24,58.36,61.24,0.35,938.34,9.88,98.27,0,0,024,29.83,32.95,0.19,44.67,49.35,0.28,61.96,68.45,0.39,975.19,10.44,100.32,0,0,025,33.72,29.61,0.23,54.55,47.89,0.38,74.13,65.08,0.51,1099.65,7.71,107.25,0,0,026,34.3,31.94,0.24,54.01,50.3,0.39,75.93,70.71,0.54,1118.15,8.47,108.28,

36、0,0,027,32.7,30.89,0.22,44.32,41.86,0.3,70.94,67.01,0.49,1067.1,8.65,105.44,0,0,028,29.35,28.16,0.08,31.84,30.55,0.09,60.48,58.02,0.17,959.98,8.84,99.47,0,0,0;% 单位统一% QN(:,4,7,10)=QN(:,4,7,10).*(109)./(24*3600); % 排放量单位ug/s（2）% 计算第daynum天，位置为(r,theta)，出烟温度为Ts的地方的Q(颗粒，SO2，NOx)daynum=111;r=1000;theta=

37、1.3;Ts=80;% function C=calcu_C(daynum,r,theta,Ts)Vwind=windreadin; % 风速风向数据theta0=2*pi/365*Vwind(:,3); % theta0degsun=0.006918-0.339912*cos(theta0)+0.070257*sin(theta0)-0.006785*cos(2*theta0)+. 0.000907*sin(2*theta0)-0.002797*cos(3*theta0)+0.001480*sin(3*theta0); % 太阳倾角jingdu=114.097586; % 经度weidu=2

38、2.686033/180*pi; % 纬度h0=asin(sin(weidu)*sin(degsun)+cos(weidu)*cos(degsun)*cos(15*8+jingdu-300)/180*pi)/pi*180; % 太阳高度角(deg)lenh=length(h0);RayRate=zeros(lenh,1);StableRate=zeros(lenh,1);for i=1:lenh RayRate(i)=calcu_RayRate(h0(i),Cloudreadin(Vwind(i,3);endfor i=1:lenh StableRate(i)=calcu_StableRate

39、(RayRate(i),Vwind(i,1); % 稳定等级(1-6 A-F)end% yita1、aerfa1第一列是1000m以下的，第二列是1000m以上的yita1=0.425809,0.6020520.281846,0.3963530.177154,0.2321230.110726,0.1466690.0864001,0.1019470.053634,0.0733348;aerfa1=0.901074,0.8509340.914370,0.8650140.924279,0.8851570.929418,0.8887230.920818,0.8968640.929418,0.88872

40、3;% yita2、aerfa2第一列是300m以下的，第二列是300-500m,第三列是500-10000，第四列是10000以上aerfa2=1.12154,1.51360,2.10881,2.108810.961435,0.961435,1.09356,1.093560.917595,0.917595,0.917595,0.9175950.826212,0.826212,0.632023,0.555360.788370,0.788370,0.565188,0.4147430.784400,0.784400,0.525969,0.322659;yita2=0.0799904,0.00854

41、771,0.000211545,0.0002115450.127190,0.127190,0.057025,0.0570250.106803,0.106803,0.106803,0.1068030.104634,0.104634,0.400167,0.8107630.0927529,0.0927529,0.433384,1.732410.0620765,0.0620765,0.370015,2.40691;sigmay=yita1(StableRate(find(Vwind(:,3)=daynum),isRlargethan1000(r)*(r*cos(theta)aerfa1(StableR

42、ate(find(Vwind(:,3)=daynum),isRlargethan1000(r);sigmaz=yita2(StableRate(find(Vwind(:,3)=daynum),Rrange(r)*(r*cos(theta)aerfa2(StableRate(find(Vwind(:,3)=daynum);d=3; % 直径 mvs=12.7; % 烟气排出速度 m/sc0=20,80,250; % 排出烟的浓度 mg/m3Q=pi/4*d*d*vs*c0;h1=80;deltah=vs*d/Vwind(find(Vwind(:,3)=daynum),1)*(1.5+2.7*(T

43、s-25)/(Ts+273)*d);H=h1+deltah;C=(Q/(pi*Vwind(find(Vwind(:,3)=daynum),1)*sigmay*sigmaz)*exp(-0.5*(r*cos(theta)2)/(sigmay2)+(H2)/(sigmaz2);（3）% 根据年月求日期序号function daynum=calcu_date(month,day)switch month case 1 daynum=day-1; case 2 daynum=30+day; case 3 daynum=30+28+day; case 4 daynum=30+28+31+day; cas

44、e 5 daynum=30+28+31+30+day; case 6 daynum=30+28+31+30+31+day; case 7 daynum=30+28+31+30+31+30+day; case 8 daynum=30+28+31+30+31+30+31+day; case 9 daynum=30+28+31+30+31+30+31+31+day; case 10 daynum=30+28+31+30+31+30+31+31+30+day; case 11 daynum=30+28+31+30+31+30+31+31+30+31+day; case 12 daynum=30+28+31+30+31+30+31+31+30+31+30+day;end （4）% 本函数计算太阳辐射等级 输入h0和总云量function RayRate=calcu_RayRate(h0,Tcloud)if(h0<=15) if(Tcloud<=4) RayRate=1; else RayRate=0; endelse if(Tcloud<=7) RayRate=1; else RayRate=0; endend（5）% 本函数计算大气稳定度，根据辐射等级和风速f