中型垃圾焚烧厂对于周边居民的潜在风险进行分析建模论文-本科论文.doc

摘要
在国内垃圾焚烧厂出现各种环境污染问题的背景下，本组针对深圳市计划建立的中型垃圾焚烧厂对于周边居民的潜在风险进行分析，确定一种可行的环境动态监控体系，并设计了合理的周围居民风险承担经济补偿方案。对于本问题的研究，我们建立高斯大气扩散模型和概率模型进行总体分析。
对于问题1，为了解决垃圾焚烧厂附近污染物浓度随空间位置变化的规律，我们建立了污染物扩散问题中最为经典的高斯大气扩散模型。首先，利用谷歌地球中的经纬度划分，把当地周围划分为多个网格。其次，根据附件4所给数据，用Excel计算11年4月至12年3月的平均风速，风向等指标。随后，使用MATLAB软件编程，根据深圳市气象局所提供的数据，包括风力，降水量，光照时间，得到了08年~10年每天白天和夜间的大气稳定度等级。最后，根据大气稳定度，确定污染物浓度函数中的各个参数，进而给出最终的每季度各污染物随空间位置变化的函数，并用MATLAB三维绘图做出函数的图像。建立模型后，为确立合理的经济补偿方案，本组人员查找相关资料，通过计算当地的空气污染指数，按照对居民环境侵害的强度，时间和地点，根据当地的经济赔偿相关标准，设计了合理的经济补偿方案。
对于问题2，通过观察附件2中焚烧炉的在线排放监测记录，很容易发现，颗粒物的浓度明显超标，而二氧化硫和氮氧化合物的浓度并不超标。对于颗粒物的监测方法，我们提出了优化措施，建立有关污染物扩散的微分方程模型，即根据透过受污染地区的透过的光的光强来近似判断颗粒物的浓度。为解决故障发生概率下的污染物浓度问题，我们建立概率模型，将排放监测记录中污染物超标的频率近似为除尘装置发生故障的概率。对3月内颗粒物的排放浓度取平均值，在考虑故障发生概率的情况下，再次使用高斯大气扩散模型，得到此时的颗粒物浓度随空间位置的变化规律。最后在发生故障的概率条件下，对问题1中的经济补偿方案进行优化，针对颗粒物所造成的空气污染指数变化实行更为合理的补偿方案。
关键词 高斯大气扩散模型，微分方程模型，概率模型

目录

问题重述 1

问题分析 1

模型假设 2

定义与符号说明 2

模型的建立与求解 3

问题1的模型：高斯大气扩散模型 3

问题2的模型：微分方程模型 11

模型的优化与评价 15

参考文献 16

附录 17

问题重述
垃圾焚烧厂对环境的污染风险与建设投资规模、运行监管力度有直接关系。由于大小的垃圾焚烧炉各有利弊，因此对于不同城市来说，究竟该把大型焚烧厂的建设规模控制在什么水平，这是一个值得研究的课题。
深圳市某地点计划建立一个中型的垃圾焚烧厂，计划处理垃圾量1950吨/天（设置三台可处理垃圾650吨/天的焚烧炉，排烟口高度80米，每天24小时运转）。从构建环境动态监控体系、并根据潜在污染风险对周围居民进行合理经济补偿的需求出发，有关部门希望能综合考虑垃圾焚烧厂对周围带来环境污染以及其他危害的多种因素（例如，焚烧炉的污染物排放量、居住点离开垃圾焚烧厂的距离、风力和风向及降雨等气象条件、地形地貌以及建筑物的遮挡程度等等），在进行科学定量分析的基础上，确立一套可行的垃圾焚烧厂环境影响动态监控评估方法，并针对潜在环境风险制定出合理的经济补偿方案。
请你在收集相关资料的基础上考虑以下问题：
(1) 假定焚烧炉的排放符合国家新的污染物排放标准（参见附件1），根据垃圾焚烧厂周边环境设计一种环境指标监测方法，实现对垃圾焚烧厂烟气排放及相关环境影响状况的动态监控。以你设计的环境动态监控体系实际监控结果为依据，设计合理的周围居民风险承担经济补偿方案。
(2) 由于各种因素焚烧炉的除尘装置（如袋式除尘器）损坏或出现其他故障导致污染物的排放增加，致使相关各项指标将严重超标（如：烟尘浓度、二氧化硫、氮氧化物、一氧化碳、二恶英类及重金属等排放超标，附件2给出了一台可处理垃圾350吨/天的焚烧炉正常运作时的在线排放监测记录）。请在考虑故障发生概率的情况下修正你设计的监测方法和补偿方案。
问题分析
(一) 问题1的分析
本问题所研究的环境指标监测方法有着较强的实际意义。通过对周边环境的实时监控，可以让居民了解自己所处的生态环境，并在必要时采取适当的保护措施。本组通过查找相关资料，运用Excel和MATLAB软件，在焚烧炉的排放符合国家新的污染物排放标准的条件下得到了在受影响区域内各个方位的污染物浓度。我们运用高斯大气扩散模型，利用了附件4中的风向和风速，借鉴了深圳市气象局的降水量和日照时间等数据，先得到处理场周围在08年至10年各月的大气稳定度平均值，再用风速，降水，大气稳定度等数据，参考大气扩散模型中的相关计算公式，得到该区域各点的污染物浓度。最后，参考附件3中的生活垃圾焚烧污染控制标准，结合以上对于居民区污染物浓度的分析，研究污染物超标可能导致的对人体的健康威胁，据此设计合理的周围居民风险承担经济补偿方案。
(二) 问题2的分析
由于问题1中所做的假设“焚烧炉的排放符合国家新的污染物排放标准”过于理想，现实情况下难以做到。由于各种因素焚烧炉的除尘装置（如袋式除尘器）损坏或出现其他故障导致污染物的排放增加，致使相关各项指标将严重超标。附件2给出了3月中烟气的日均排放量，通过计算各污染物严重超标的频率，利用概率模型，可以得到故障发生的概率。在可能发生故障的情况下，再次利用高斯大气扩散模型，在第一问中监测方法可行的条件下，计算出各项指标严重超标时各地的污染物浓度，进而设计出更优化的补偿方案。
模型假设
1. 假设题目所给的数据真实可靠；
2. 在全部高度上，风速是均匀且稳定的。大气污染物与空气没有相对运动，污染物的运动是被动的，且在污染物扩散的过程中，污染物的总质量没有发生变化，符合质量守恒定律，而且污染物本身的化学组成不会与其周围大气中的物质发生化学反应。大气污染物的扩散将在地表被反射出去。
3. 大气污染物浓度在 Y 方向和 Z 方向上分布满足正态分布；
4. 污染源的污染物浓度分布是连续、均匀、稳定的；
定义与符号说明

rain_i_j
第i年第j个月的降雨量
s_i_j
第i年第j个月的光照时长
R_i
第i年的降雨量
s_i
第i年的光照时长
S_i
第i年的太阳辐射强度（0代表弱，1代表中，2代表强）
w_i
第i个月的风速
W
全年风速
day_i
第i年白天的大气稳定度
night_i
第i年晚上的大气稳定度
cloud_i
第i年云的情况
byi
大气稳定度为i时的值
bzij
大气稳定度为i时的值(j=1，2)
ave_day
按季度的白天大气稳定度平均值
ave_night
按季度的晚上大气稳定度平均值
ave_w
按季度的风力的平均值
q
各项污染源强度
c
运用高斯模型计算出的浓度表达式
IAQIP
污染物项目P的空气质量分指数
CP
污染物项目P的质量浓度值
BPHi
相应地区的空气质量分指数及对应的污染物项目浓度指数表中与CP相近的污染物浓度限值的高位值
BPLo
相应地区的空气质量分指数及对应的污染物项目浓度指数表中与CP相近的污染物浓度限值的低位值
IAQIHi
相应地区的空气质量分指数及对应的污染物项目浓度指数表中与BPHi对应的空气质量分指数
IAQILo
相应地区的空气质量分指数及对应的污染物项目浓度指数表中与BPLo对应的空气质量分指数

模型的建立与求解
问题1的模型：高斯大气扩散模型
高斯大气污染扩散模型定义如下参数和条件：
1、 高斯坐标系
原点设定为污染源排放中心点在地表的投射点，X 轴的正方向为该地区的平均风
向，Y 轴的正方向为该地区平均风向的切向，在地表的水平面上垂直于 X 轴，Z 轴为垂直于水平地表面 XOY，正方向向上。
高斯大气污染物扩散模型的坐标系如图2所示：

图1
2、 理想状态下，高斯污染扩散模型建立的假设
a) 大气污染物浓度在 Y 方向和 Z 方向上分布满足正态分布；
b) 在全部高度上，风速是均匀且稳定的；
c) 污染源的污染物浓度分布是连续、均匀、稳定的；
d) 扩散中，污染物的质量是守恒的（不考虑转化等其他因素）；
在本题中，需要考虑地形及建筑物等因素对污染物扩散的影响。运用谷歌地球，得到当地周围的地形等高线，近似表示为下图：

图2
在地形图中可以看出，顺风方向地势较为平坦，海拔较低，居民区较为集中。在本问题中，地势因素对污染物浓度的影响可以近似忽略。
定义：
我们假定大气污染物的浓度c，而且假定污染物的浓度与其所在的地理位置有关；
大气污染源的强度设为 Q；
沿 x 轴正方向上的风速设为 u。
首先，我们需要利用谷歌地球的经纬度，将污染扩散的预测地区进行网格划分，将该区域划分成a*b的矩阵，矩阵中的每一个点对将对应着一个实际的地理坐标。我们只需要计算出每一个坐标点上的大气污染物浓度，就完成了对该地区大气污染程度的预测。
然后，我们还需要创建大气扩散的参数矩阵，以便保存高斯大气污染扩散模型的计算中所需要的参数。
第三，我们将读取预测区域当前的气象条件，通过用户输入污染源强度、当地平均风速、所要计算的大气的高度和当时大气稳定程度，来读入当前气象条件。然后根据读入的实际的大气状态来确定高斯系数矩阵中的相应系数计算网格中每一测量点上的污染状况。
具体流程如下图所示：

图3：高斯大气扩散模型的设立
由正态分布假定，可以推定下风向中任意一点的污染物浓度分布，因此大气污染物的浓度 c 将满足下式：

（1）
其中，方差形式的表达式为：
（2）
（3）
由于单位时间内污染物的质量守恒，由假定可以推出d污染源强的积分公式：
（4）
4 式中的未知数为：c 代表污染物的浓度，函数 A(x)待定，系数 a，b 待定，通过查表或者将式 4 按照级数展开可得：

再代入（2）式，可得：
（5）
将 1 式、5 式、6 式代入 4 式中，可得：
（6）
将 5 式、6 式代入 1 式中，可得高斯大气扩散模型的最终表达式：

（7）
由 7 式中的高斯大气污染物浓度表达式可以知道，污染物的浓度除了受到空间位置的影响之外，还会受到大气污染源的强度 q、沿 X 轴正方向的平均风速 u 以及污染扩散系数a,c,f等因素的影响。
本系统中，高斯大气污染扩散模型的实现采用的是 matlab 技术，由于 matlab具有高效率的数值计算功能和符号运算功能，可以把编程人员从复杂的数学运算分析过程中解放出来，同时 matlab 还具有比较完备的图形图像处理功能，可以很方便地实现计算结果和编程的可视化，由于 matlab 的语法和思维方式比较接近人类的自然语言,使得初学者容易学习和应用，而且 matlab 提供了功能强大的一些实用的工具，为编程人员提供了大量便于使用的处理工具，所以本项目使用 matlab软件进行编程实现。
3、 高斯污染物扩散模型的参数的确定
高斯污染物扩散模型的参数是描述大气污染物沿不同方向的扩散程度，其计算公
式为：

其中参数 a、c、f 是根据不同地区的大气稳定程度而进行确定取值，不同的地区、
季节，大气稳定度系数不同。为计算大气稳定度，本组人员对08年至10年的降水量，光照时间，风速进行综合分析，确定了大气稳定度等级的判别方法，具体方法如下：

地面风u (m/s)
白天太阳辐射
阴天的白天或夜间
有云的夜间
强
中
弱
薄云遮天或低云大于等于5/10
云量小于5/10
<2
A
A—B
B
D

2—3
A—B
B
C
D
E
F
3—5
B
B—C
C
D
D
E
5—6
C
C—D
D
D
D
D
>6
C
D
D
D
D
D

判定了三年内每天的白天和夜间的大气稳定度之后，按照季度去平均值，得到每个季度的大气稳定度平均水平。本项目中对大气稳定程度的分类采用国际上比较常用的帕斯奎尔分类法，这种分类法将空气稳定度分为六个等级：分别为极不稳定、不稳定、较不稳定、中性、较稳定和稳定六个级别，从不稳定到稳定的等级分别为 A、B、C、D、E、F。我们根据不同的大气稳定度，来最终确定计算大气污染状况模型中所采用的参数，详见表1、表2和表3:

表 1：σy扩散系数表

A
B
C
D
E
F
a
213
156
104
68
50.5
34
f
0
0
0
0
0
0
c
0.894
0.894
0.894
0.894
0.894
0.894

表2：σz（x>1000）扩散系数表

A
B
C
D
E
F
a
440.8
106.6
61
33.2
22.8
14.35
f
-9.27
-3.3
0
1.7
1.3
0.35
c
1.941
1.149
0.911
0.725
0.678
0.74

表3: σz（x<=1000）扩散系数表

A
B
C
D
E
F
a
459.7
108.2
61
44.5
55.4
62.6
f
9.6
-2
0
13
34
48.6
c
2.094
1.098
0.911
0.516
0.305
0.18

运用以上数据，可以得到参数σy、σz。
对于污染源强q，根据附件1中各污染物的标准浓度和附件2中3月份的平均排放量，得到各种污染物每秒的排放量，即为源强。
最后，将各参数代入高斯大气扩散模型的最终表达式，即可得到周边各位置的污染物浓度c(x,y,z)。
对于颗粒物，本组成员运用MATLAB软件编辑三维作图程序，将其浓度的表达式化简，得到：
C=(1432975.*exp(-1./(21632.*X.^(447./250))))./(32684288.*pi.*X.^(361./200));
则其浓度随空间变化的图为：

图4
对于第一季度的氮氧化合物浓度，我们运用同样的方法，得到其浓度的变化函数为：
C=(615458076098559969.*exp(-1./(21632.*X.^(447./250))))./(1123023584424361984.*pi.*X.^(361./200))
浓度随空间变化的图像：

图5
对于第一季度的二噁英浓度，其函数为：
C=(57319.*exp(-1./(21632.*X.^(447./250))))./(261474304.*pi.*X.^(361./200))
图像如下：

图6
4、 各种污染物的监测方法：
通过查找相关资料，本组提出了包括颗粒物、HCL、SO2、NOx、汞、铅和二噁英七种污染物的动态监控体系和实时监测方法，具体如下：
（1）颗粒物（包括PM2.5和PM10）：β射线法、微量震荡天平法；
(2) 二氧化硫：紫外荧光法、差分吸收光谱分析法 ；
（3）氮氧化合物：化学发光法、差分吸收光谱分析法
（4）氯化氢：滤膜采样氟离子选择电极法
（5）Hg：巯基棉富集-冷原子荧光分光光度法
（6）铅：火焰原子吸收分光光度法
（7）二噁英：色谱发、激光质谱法、免疫法和生物法
5、 赔偿标准的设计
本问题中的侵害属于环境侵害，它在本质上是一种侵权行为，但却是与传统的侵权行为有所不同的、具有独特性质的新型侵权行为。
在传统侵权行为中侵害的对象一般为特定的人或物, 而在环境侵害行为中, 环境侵害的对象常常是相当范围内不特定的多数人或物, 具有社会性。由于环境侵害与经济发展具有相伴相随的孪生关系, 造成环境侵害的原因事实其本身常常是各种创造社会财富、增进公众福利的活动在进行过程中的附带行为。因此环境侵害具有相当程度的价值正当性或社会有用性的特点,使得环境侵害必将长期存在、具有持续性特点。在环境侵害形成的过程中, 有时单一的排污行为或环境开发行为还不至于会造成环境侵害, 但众多的上述行为掺合在一起最终就形成了环境侵害。
对人身伤害应赔偿由此造成的损失历来是确定人身伤害赔偿的一般原则, 即由侵权行为造成他人人身伤害, 以因伤害人身而引起的财产损失作为标准, 损失多少财产就赔偿多少。污染环境所造成的人身损害多数具有潜在性和隐蔽性, 受害人往往在开始受害时显露不出明显的损害, 但随着时间的推移, 损害逐渐显露, 如早衰、人体功能减退等。对于这种潜在的危害虽未经医治、尚未损失医疗费等也应当予以赔偿, 以补偿对人体的潜在损害。
在环境侵权中, 有不少权利的损害很难归类于财产损害和人身伤害, 而只能归类于精神损害, 诸如公民环境权中的日照权、通风权、安宁权等。这些权利是环境权不可分割的部分,因而当这些权利受到侵害时, 应当予以赔偿。精神损害具有特殊性, 在考虑赔偿时应着重注意以下几方面:
第一, 侵害强度。精神损害往往因人而异,主观性较强, 难以具体度量, 而侵害强度则是有形的、可测量的, 在很大程度上能反映精神损害程度, 直接影响侵权人是否承担赔偿责任和赔偿的数额。侵害强度可参照环境标准确定。
第二, 侵害时间。通常持续性的或间歇短的环境侵害应承担责任, 一时性的、间歇长的可不承担。侵害时间还要区分是在夜间还是白天, 例如相同强度噪音对居民造成的侵害夜晚要大于白天。
第三, 侵害地点。侵害行为是在生活区还是工业区也影响损害赔偿标准。生活区对环境的要求相对较高, 工业区环境污染控制标准相对较松, 因而在生活区对环境权的侵害要承担比在工业区更重的责任。
环境侵害的精神赔偿标准应综合考虑侵害强度、侵害时间、侵害地点。显然, 在生活区对公民进行持续性的高强度侵害要承担最重的责任。为了使精神损害赔偿更加具体化甚至定型化, 可以先将侵害强度、侵害时间、侵害地点分成不同的等级, 将侵害强度、侵害时间、侵害地点赋予一个不同的加权数(比例),然后算出结果, 最后赋予不同的结果以不同的赔偿额。
侵害程度(赔偿数额)=（空气污染指数*权数(1/2)+侵害时间*权数(3/10)+侵害地点*权数(1/5)）*单位侵害程度赔偿数额（λ）。
即：P=（AQI*50%*+T*30%+D*20%）*λ
注：λ的取值依据不同地区的相关政策与发展水平确定
其中空气污染指数AQI计算方法如下：
空气质量分指数
对照各项污染物的分级浓度限值，以细颗粒物（PM2.5）、可吸入颗粒物（PM10）、二氧化硫（SO2）、氮氧化物（NOx）、铅（Pb）、二恶英、汞（Hg）、氯化氢（HCL）等各项污染物的实测浓度值（其中PM2.5、PM10为24小时平均浓度）分别计算得出空气质量分指数

式中：
IAQIP——污染物项目P的空气质量分指数；
CP——污染物项目P的质量浓度值；
BPHi——相应地区的空气质量分指数及对应的污染物项目浓度指数表中与CP相近的污染物浓度限值的高位值；
BPLo——相应地区的空气质量分指数及对应的污染物项目浓度指数表中与CP相近的污染物浓度限值的低位值；
IAQIHi——相应地区的空气质量分指数及对应的污染物项目浓度指数表中与BPHi对应的空气质量分指数；
IAQILo——相应地区的空气质量分指数及对应的污染物项目浓度指数表中与BPLo对应的空气质量分指数。
空气质量指数:
AQI = max {IAQI1, IAQI2, IAQI3···IAQIn}
式中：IAQI——空气质量分指数；
n——污染物项目
简单来说，AQI就是在各IAQI中取最大值。
AQI大于50时，IAQI最大的污染物为首要污染物。若IAQI最大的污染物为两项或两项以上时，并列为首要污染物。IAQI大于100的污染物为超标污染物。

空气质量分指数（IAQI）
污染物项目限制
二氧化硫（SO2）24小时平均/（μg/m3）
Hg（最高允许排放速率，kg/h)
氮氧化物（Nox）24小时平均/（μg/m3）
Pb（最高允许排放速率，kg/h)
颗粒物（粒径小于等于10μm）24小时平均/（μg/m3）
HCL（最高允许排放速率，kg/h)
二恶英 (ng TEQ/Nm3)
0—100
0—150
0—0.03
0—80
0—0.1
0—225
0—6
0—0.1
100—200
150—800
0.03—0.05
80—280
0.1—0.2
225—500
6—12
0.1—0.5
200—300
800—1600
0.05—0.0
280—565
0.2—0.3
500—670
12—19
0.5—1
>300
>1600
>0.07
>565
>0.3
>670
>19
>1

问题2的模型：微分方程模型
通过对附件2中烟气排放连续监测日均月报表的分析，很容易发现，3月中大部分时间段颗粒物的排放量超过了国家标准，而二氧化硫和氮氧化合物的排放量一直稳定在国家标准以下，我们首先用Excel将全月的颗粒物，二氧化硫，氮氧化合物排放情况和国家标准分别进行作图，图像如下：

图7

图8

图9

通过查找相关资料，建立了烟雾扩散理论中经典的微分方程模型，通过分析透过污染地点的光线强度变化来分析颗粒物浓度的变化。
1、模型假设
（1）烟雾在无穷空间扩散，不受大地和风的影响；扩散服从热传导定律。
（2）光线穿过烟雾时光强的减少与烟雾浓度成正比；无烟雾的大气不影响光强。
（3）穿过烟雾进入仪器的光线只有明暗之分，明暗界限由仪器灵敏度决定。
2、模型建立
（1）穿过烟雾的光强的变化规律：
设I(x)：沿x方向的光强；C(x):沿x方向的颗粒物浓度。
根据经验，光强的减少大致与烟雾浓度成正相关，即：

设I（l0）=I0,，即未进入烟雾的初始光强为I0，
则求解简单的微分方程初值问题，得到上式的解：

（2）仪器灵敏度与烟雾明暗界限
由于烟雾浓度连续变化，烟雾中光强连续变化，导致不透光区域有扩大，缩小的过程。穿过烟雾进入仪器的光线只有明暗之分，明暗界限由仪器灵敏度决定。
设μ为仪器灵敏度，当I/I0<1-μ，观测结果为暗。
设光源在z=-¥, 仪器在z=+¥,则观测到的明暗界限为：

示意图如右图所示：
仪器
z
-¥
¥

图10
由附件2所提供的数据，经过Excel的数据分析发现，只有2号，3号和4号的颗粒物排放量达标，其余日期的排放量全部超标。由于超标的程度各不相同，我们对全月的排放量取平均值，其平均值为26.07 mg/m3，再算得3月的平均风速，将其代入问题1中的高斯大气扩散模型，得到此时的污染物浓度与空间位置的关系：

图 11
白天的浓度分布函数为：
c_day=(q/(2*pi*u*by4*bz41))*exp(-y^2/(2*by4^2))

图 12
晚上的浓度分布函数为：
c_night=(q/(2*pi*u*by3*bz31))*exp(-y^2/(2*by3^2))
由于颗粒物的浓度超标，将给周围居民的生活造成不良影响。大气中的颗粒物浓度对人体健康的危害是不完全一样的。不言而喻，浓度越高，危害越严重。此外，颗粒物的大小，对人体健康的危害程度也不一样。颗粒愈小，对人体的危害愈大。100微米以上的颗粒物，在大气中滞留时间很短，可以很快沉落地面。对人体健康影响不普遍。100微米以下的颗粒物，可以长期在空气中滞留，对人体健康和空气的能见度，有不良影响。10微米以下的颗粒物（PM10），可以进入呼吸道，引发呼吸道疾病。2.5微米以下的颗粒物（PM2.5），可进入肺泡，引发肺气肿等。1.0 微米以下的颗粒物（PM1），通过肺泡，进入血液系统，引发多种疾病。
基于颗粒物排放比较严重，厂家应当附加一定的居民对颗粒物的防护费用，例如购买口罩，除尘器，洗车等费用。
侵害程度(赔偿数额)=（空气污染指数*权数(1/2)+侵害时间*权数(3/10)+侵害地点*权数(1/5)）*单位侵害程度赔偿数额（λ）+居民防护费
模型的优化与评价
优点：
本文较为全面的考虑了各种因素对于污染物扩散的影响，包括降水，光照，风向，风速，大气稳定度等。对污染物浓度的计算较为准确。高斯大气扩散模型作为研究污染物扩散中最为经典，成熟的模型之一，考虑的因素较为全面，分析结果能够与实际结果较为吻合。由于本问题中的主要研究对象是平原地区，没有较大的地形落差，因此可以将地形因素忽略不计。此外，在现实的复杂情况下，地型数据和气象信息往往都是海量的，这样将会消耗很大的计算时间和空间，而本组利用高斯大气污染扩散模型公式进行计算，可以大大减少系统的计算量。在第二问中，我们所提出的用透过污染区域的光线的强度来估计颗粒物浓度的方法比较新颖，有较强的可行性。
不足：
由于时间紧迫以及数据量不足，我们在建立模型的过程中将一些影响较小的因素忽略，使得模型有点理想化，与实际情况可能有一定的偏差。而且由于条件有限及技术层面不够成熟，计算的结果可能不完全准确。

参考文献
[1] 维基百科：http://zh.wikipedia.org/wiki；
[2] 姜启源，谢金星，叶俊.数学模型[M].北京：高等教育出版社2011.1 第四版；
[3] 刘卫国.MATLAB程序设计教程[M].北京:中国水利水电出版社,2010.2 第二版；
[4] 张圣勤.MATLAB7.0实用教程[M].北京：机械工业出版社，2006.3 第一版； [5] 欧阳俊强.长春市环保局大气污染模拟系统的设计与实现[F].吉林大学，2013年；
[6] 中华人民共和国国家环境保护标准HJ633—2012；
[7] 中华人民共和国国家大气污染物综合排放标准;
[8] 赵胜军.环境侵害赔偿标准的具体化研究. 武汉：武汉大学环境法研究所，2000年；
[9] 曾书儿. 风速风向的矢量平均方法. 气象科学研究院大气探测研究所，1983年；
[10] 杨新兴，冯丽华，尉鹏. 大气颗粒物PM2.5及其危害. 北京：中国环境科学研究院，2012年。

附录
%大气稳定度的计算
%08年白天的大气稳定度
for i=1:31
for j=1:12
if s_08(i,j)<4
S_08(i,j)=0
end
if s_08(i,j)>=4&s_08(i,j)<8
S_08(i,j)=1
end
if s_08(i,j)>=8
S_08(i,j)=2
end
end
end
for i=1:31
for j=1:12
if R_08(i,j)>=1
day_08(i,j)=3
end
if R_08(i,j)<1
if S_08(i,j)==2
if W(i,j)<3
day_08(i,j)=6
end
if W(i,j)<5&W(i,j)>=3
day_08(i,j)=5
end
if W(i,j)>=5
day_08(i,j)=4
end
end
if S_08(i,j)==1
if W(i,j)<3
day_08(i,j)=5
end
if W(i,j)<5&W(i,j)>=3
day_08(i,j)=4
end
if W(i,j)>=5
day_08(i,j)=3
end
end
if S_08(i,j)==0
if W(i,j)<2
day_08(i,j)=5
end
if W(i,j)<5&W(i,j)>=2
day_08(i,j)=4
end
if W(i,j)>=5
day_08(i,j)=3
end
end
end
end
end
%08年晚上的大气稳定度
for i=1:30
for j=1:12
if S_08(i+1,j)==1||S_08(i+1,j)==2
cloud_08(i,j)=0
end
if S_08(i+1,j)==0
if R_08(i,j)<1
cloud_08(i,j)=1
end
if R_08(i,j)>=1
cloud_08(i,j)=2
end
end
if R_08(i+1,j)>1
if S_08(i,j)==1||S_08(i,j)==2
cloud_08(i,j)=1
end
else
cloud_08(i,j)=2
end
end
end
for i=31
for j=1:12
if R_08(i,j)>=1
cloud_08(i,j)=1
end
if R_08(i,j)<1
if S_08(i,j)==0
cloud_08(i,j)=2
else
cloud_08(i,j)=3
end
end
end
end
for i=1:31
for j=1:12
if W(i,j)>5
night_08(i,j)=3
end
if cloud_08(i,j)==0
night_08(i,j)=3
end
if cloud_08(i,j)==1
if W(i,j)<3
night_08(i,j)=1
end
if W(i,j)>=3&W(i,j)<5
night_08(i,j)=2
end
end
if cloud_08(i,j)==2
if W(i,j)<3
night_08(i,j)=2
end
if W(i,j)>=3
night_08(i,j)=3
end
end
end
end
%09年白天大气稳定度
for i=1:31
for j=1:12
if s_09(i,j)<4
S_09(i,j)=0
end
if s_09(i,j)>=4&s_09(i,j)<8
S_09(i,j)=1
end
if s_09(i,j)>=8
S_09(i,j)=2
end
end
end
for i=1:31
for j=1:12
if R_09(i,j)>=1
day_09(i,j)=3
end
if R_09(i,j)<1
if S_09(i,j)==2
if W(i,j)<3
day_09(i,j)=6
end
if W(i,j)<5&W(i,j)>=3
day_09(i,j)=5
end
if W(i,j)>=5
day_09(i,j)=4
end
end
if S_09(i,j)==1
if W(i,j)<3
day_09(i,j)=5
end
if W(i,j)<5&W(i,j)>=3
day_09(i,j)=4
end
if W(i,j)>=5
day_09(i,j)=3
end
end
if S_09(i,j)==0
if W(i,j)<2
day_09(i,j)=5
end
if W(i,j)<5&W(i,j)>=2
day_09(i,j)=4
end
if W(i,j)>=5
day_09(i,j)=3
end
end
end
end
end
%09年晚上大气稳定度
for i=1:30
for j=1:12
if S_09(i+1,j)==1||S_09(i+1,j)==2
cloud_09(i,j)=0
end
if S_09(i+1,j)==0
if R_09(i,j)<1
cloud_09(i,j)=1
end
if R_09(i,j)>=1
cloud_09(i,j)=2
end
end
if R_09(i+1,j