本发明属于大气环境污染物检测技术领域，尤其涉及中尺度大气光化学污染的改进型预测模型。
背景技术：
空气质量模式，是在对污染物排入大气环境后传输、扩散、转化和清除等一系列物理和化学过程的认识基础上，利用气象、环境、物理、化学等学科的研究方法和计算机技术，实现模拟和预报不同空间尺度上空气污染物浓度分布状况及变化趋势的方法，在空气质量预报、大气污染控制、环境规划与管理、城市建设及公共卫生等方面均有重要的实际应用价值，具有广阔的发展前景。
CALGRID是由美国加州的ARB(Air Resources Board)开发的，是欧拉型的中尺度大气光化学模式，针对二次污染物如臭氧的模拟具有较好效果，主要适用于晴空条件下的光化学反应的模拟，包含了大气输送与扩散、气相化学反应、人为排放的点面线源、干沉降等过程。
随着城市化的进程加快和社会经济的快速发展，以酸雨、光化学污染、灰霾污染物相互耦合叠加为特征的大气复合污染成为了我国越来越突出的大气环境问题。近些年，对京津冀、长江三角洲、珠江三角洲等重点城市群的区域大气复合污染的形成机制和综合防治的研究是我国环境保护研究的热点问题。因此，要求空气质量模式能对多种大气污染物(痕量气体、气溶胶等)在不同尺度下(区域、城市等)的不同类型污染过程(气相、液相、非均相)进行模拟。为了更好的利用CALGRID模型认识和理解痕量气体的气相化学过程、气溶胶微物理化学过程和液相化学过程，为大气复合污染调控与治理提供了理论支持，需要CALGRID不仅要能模拟解决与气相化学反应有关的臭氧污染问题，还要能合理描述气溶胶的物理化学过程，云雨水中的液相化学过程，以及与降水有关的湿沉降过程。
此外，真实大气中各种化学成分以及各类化学反应过程是相互影响的。例如，臭氧和颗粒物之间，臭氧会影响无机盐气溶胶前体物SO2、NO2向硫酸盐、硝酸盐转化，会影响VOCs生成二次有机气溶胶(SOA)的过程；另一方面，颗粒物会直接影响大气辐射，作为云凝结核(CCN)形成云滴间接影响大气辐射，这些都会影响气相成分的光化学反应，并且颗粒物表面的非均相化学过程也会影响大气气相成分的浓度。可见，CALGRID模型只有包括了各个化学过程，才能合理准确的模拟预报臭氧、颗粒物等大气成分的浓度和沉降量。
有机气溶胶是大气气溶胶的重要成分，在污染严重的城市地区一般占PM2.5和PM10质量的20～60％，而在偏远地区大约占PM10的30～50％。无论在污染地区还是在偏远的背景地区，有机气溶胶都是由数百种有机化合物组成的混合物，其中许多具有致癌、致畸和致突变性，如多环芳烃、多氯联苯和其它含氯有机化合物。它们还能够影响大气能见度，是大气光化学烟雾、酸沉降的重要贡献者，可通过长距离传输对区域和全球环境产生重要影响。随着经济的高速发展，中国出现了城市和区域性大气颗粒物污染现象，有机气溶胶日益成为大气污染控制的关键污染物和控制的难点。因此，有必要在CALGRID中补充增加有机气溶胶的模拟能力。
技术实现要素：
针对现有技术存在的问题，本发明提供了一种增加有机气溶胶化学过程的大气污染模拟预测算法,有效的对大气颗粒物污染，尤其是对大气中挥发性有机物的气溶胶浓度变化进行了模拟预测，对大气污染的控制和防治提高了理论依据和指导。
为解决上述技术问题，本发明采用了以下技术方案：一种增加有机气溶胶化学过程的大气污染模拟预测算法，包括以下步骤：
步骤1：基于CALGRID的化学模式，考虑大气化学反应、大气输送与扩散、沉降、地面面源和高架排放源的影响，并引入对有机气溶胶浓度的影响变化项，得到化学物种浓度变化方程如式(1)，
式中，右侧多项式中，第一项为二阶湍流扩散项，第二项为扩散项，第三项(P-L)GAS为气相化学变化项，第四项CHEMorg为有机气溶胶浓度的变化项，第五项EANT为人为污染源排放项，第六项由干沉降引起的物种浓度变化；其中，C是化学物种平均浓度，V是三维风矢量的平均量，K是湍流扩散系数，式中，二阶湍流扩散项通过湍流扩散系数K理论闭合转化得到；
步骤2：有机气溶胶浓度的变化项CHEMorg的求解方法如下：
假定活泼性有机气体ROG通过光化学反应生成八类半挥发有机气体，并通过气粒转换成相应的有机气溶胶SOA，对所述八类活泼有机气体ROG进行种类划分，并对其转化的有机气溶胶SOA进行分类，各转换反应的参数列于表1：
上述活泼有机气体在大气中与OH自由基，NO3自由基和臭氧进行氧化反应，反应式如下：
TOL+OH→0.039CVARO1+0.108CVARO2(1)
XYL+OH→0.039CVARO1+0.108CVARO2(2)
CSL+OH→0.039CVARO1+0.108CVARO2(3a)
CSL+NO3→0.039CVARO1+0.108CVARO2(3b)
HC+OH→0.048CVALK1(4)
OLT+OH→0.008CVOLE1(5a)
OLT+NO3→0.008CVOLE1(5b)
OLT+O3→0.008CVOLE1(5c)
OLI+OH→0.008CVOLE1(6a)
OLI+NO3→0.008CVOLE1(6b)
OLI+O3→0.008CVOLE1(6c)
API+OH→fOH·(0.028CVAPI1+0.241CVAPI2)(7a)
API+NO3→fNO3·(0.028CVAPI1+0.241CVAPI2)(7b)
API+O3→fO3·(0.028CVAPI1+0.241CVAPI2)(7c)
LIM+OH→fOH·(0.163CVLIM1+0.247CVLIM1)(8a)
LIM+NO3→fNO3·(0.163CVLIM1+0.247CVLIM1)(8b)
LIM+O3→fO3·(0.163CVLIM1+0.247CVLIM1)(8c)
上述反应方程中，权重因子fOH＝0.228，fNO3＝0.0，fO3＝0.772；其他部分氧化产物因为具有较低的饱和蒸气压，因此能通过核化或凝结过程转换成气溶胶SOA。
进一步的，步骤2中所述八类半挥发有机气体SOA的反应产物中的半挥发性产物i的总质量浓度Ctot,i通过式(2)求取：
式(2)中，ai为半挥发性产物i的化学计量系数；mi和mROG分别为半挥发性产物i和活泼性有机气体ROG的摩尔质量，gmol-1；ΔROG为反应消耗掉的ROG的量。
进一步的，所述半挥发性产物i的总质量浓度，为气态部分浓度Cgas,i和颗粒态的浓度Caer,i之和，半挥发性产物i的气溶胶态浓度即为颗粒态浓度，所述半挥发性产物i的颗粒态的浓度Caer,i，通过式(3)求解：
式(3)中，Cgas,i和Caer,i分别为气态部分和颗粒态的浓度，Ctot,i为半挥发性产物i的总质量浓度，为半挥发性纯产物i的饱和浓度，Cinit为任何额外的气溶胶态的吸收物质浓度，minit为任何额外的气溶胶态的摩尔质量，为吸收性气溶胶的浓度，下角标n为气溶胶类型序号，下角标j为气溶胶类型总数。
进一步的，由于半挥发性产物之间相互影响，需要考虑混合气体的饱和蒸气压，在此将混合的半挥发性气体之间的相互影响，假定成颗粒态物质对气态产物的吸收过程，在满足热力学平衡时，半挥发性产物i的气态部分浓度Cgas,i等于饱和蒸气压Csat,i，从而通过式(4)求解半挥发性产物i的饱和蒸气压Csat,i：
式(4)中，Xiom为摩尔分数；γi为活性系数；R为气体常数，8.314J mol-1K-1；T为温度，K；mi为半挥发性产物i的摩尔质量；为其饱和蒸气压，为纯产物i的饱和浓度。
进一步的，式(4)中，各半挥发性气体的活性系数γi为1。
进一步的，式(4)中，半挥发性产物i的饱和蒸气压通过式(5)求解：
式中，为纯物种i在参考温度下的饱和蒸气压，ΔHvap,i为物种i的蒸发晗变。
进一步的，式(4)中，半挥发性产物i的摩尔分数Xi,om通过式(6)求解：
式中，Cinit为任何额外的气溶胶态的吸收物质浓度，为吸收性气溶胶的浓度，下角标j为气溶胶类型序号，下角标n为气溶胶类型总数。
进一步的，基于活泼有机气体ROG的消耗量、初始有机气溶胶浓度、环境温度和气压的影响因子，对各影响因子进行分档，并组合形成多种条件情形，并对每种条件情形下的有机气溶胶浓度变化进行计算，得到SOA的生成率，从而建立SOA查算表数据库。
图说明
图1为本发明实施例模式的嵌套网格设置示意图；
图2为本发明实施例CALGRID模式下模拟区域及其地形图；
图3为本发明实施例CALGRID模式下模拟的2006年珠江三角洲区域黑炭的颗粒物浓度图；
图4为本发明实施例CALGRID模式下模拟的2006年珠江三角洲区域有机颗粒物浓度图；
图5为本发明实施例中耦合法与数据查表法两种方式模拟的SOA浓度曲线比较图。
具体实施方式
下面结合附图并以具体实施例，进一步阐明本发明。应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
本研究主要是基于CALGRID模式，并增加有机气溶胶浓度的影响变化项，如下式(1)：
式(1)中在原CALGIRD模式的基础上，增加有机气溶胶浓度的变化项CHEMorg，CHEMorg的求取通过以下方法进行：
(一)有机气溶胶浓度SORGAM算法
假定活泼性有机气体(ROG)通过光化学反应生成8类半挥发性有机气体，它们通过气粒转换形成相应的SOA(见表1)。这8类半挥发性气体中4类来自人为活动排放的VOC物种的氧化(其中2类来自芳香烃化合物，1类来自烷烃，1类来自烯烃)，4类来自自然源贡献，也就是α蒎烯类(API)和柠檬烯类(LIM)，它们分别氧化生成2类半挥发性气体。
表1SORGAM中考虑的有机物和二次有机气溶胶种类
活泼性有机气体(ROG)在大气中主要被OH自由基，NO3自由基和臭氧氧化，氧化反应式如下：
TOL+OH→0.039CVARO1+0.108CVARO2 (1)
XYL+OH→0.039CVARO1+0.108CVARO2 (2)
CSL+OH→0.039CVARO1+0.108CVARO2 (3a)
CSL+NO3→0.039CVARO1+0.108CVARO2 (3b)
HC+OH→0.048CVALK1 (4)
OLT+OH→0.008CVOLE1 (5a)
OLT+NO3→0.008CVOLE1 (5b)
OLT+O3→0.008CVOLE1 (5c)
OLI+OH→0.008CVOLE1 (6a)
OLI+NO3→0.008CVOLE1 (6b)
OLI+O3→0.008CVOLE1 (6c)
API+OH→fOH·(0.028CVAPI1+0.241CVAPI2) (7a)
API+NO3→fNO3·(0.028CVAPI1+0.241CVAPI2) (7b)
API+O3→fO3·(0.028CVAPI1+0.241CVAPI2) (7c)
LIM+OH→fOH·(0.163CVLIM1+0.247CVLIM1) (8a)
LIM+NO3→fNO3·(0.163CVLIM1+0.247CVLIM1) (8b)
LIM+O3→fO3·(0.163CVLIM1+0.247CVLIM1) (8c)
上式中的权重因子fOH＝0.228,fNO3＝0.0,fO3＝0.772，其部分氧化产物具有较低的饱和蒸气压，因此能通过核化、凝结等过程转换为气溶胶(SOA)。
半挥发性产物的i总质量浓度Ctot,i(μg/m-3)可以表示为：
式中，αi为上述反应式中半挥发性产物i的化学计量系数；mi和mROG分别为半挥发性产物i和活泼性有机气体ROG的摩尔质量(g mol-1)，ΔROG为反应消耗掉的ROG的量(μg/m-3)。半挥发性产物经过气粒转换，部分变成颗粒态(即SOA)，部分保留气态，即有：
Ctot,i＝Cgas,i+Caer,i (3-1)
式中，Cgas,i为在气态部分的浓度，Caer,i为颗粒态的浓度，单位均为ug/m3。在SORGAM模式中，假定低半挥发性产物之间是相互影响的，因此，需要考虑混合气体的饱和蒸气压。当一旦有气溶胶存在时，SORGAM采用一种半理想的解决方案，假定气粒转换为颗粒态物质对气态产物的吸收过程。在这种半理想的方案中，低挥发性性物种i的饱和蒸气压Csat,i(μg/m-3)可以表示为：
上式中，Xi,om是摩尔分数，γi是活性系数，R是气体常数(8.314J mol-1K-1),T是温度(K)，mi是产物i的摩尔质量，pio是其饱和蒸气压，C*sat，i是纯产物i的饱和浓度，因子106是单位转换系数。对于活性系数γi，目前研究比较少，因此在该模式中暂时全部假定γi＝1。饱和蒸气压pio是温度T的函数，满足克劳修斯-克拉贝龙方程：
上式中，pioR是纯物种i在参考温度ToR下的饱和蒸气压，ΔHvap,i是蒸发晗变。一般而言，在吸收的有机气溶胶中有n种有机物质，因此摩尔分数Xi,om表示为：
上式中，Cinit表示任何额外的气溶胶态的吸收物质浓度，minit为其摩尔质量。当满足热力学平衡时，气态的浓度就等于其饱和时的浓度，
联合上式，低挥发性物种i气溶胶态的浓度可以通过一组n维的非线性方程组表示为：
由于吸收气溶胶的浓度是开始不知道的，因此求解该方程组只能通过反复迭代的方法。
另一方面，当没有气溶胶态的有机吸收物质时，模式假定低挥发性产物的浓度必须高于一定的浓度阈值以后才会转换为气溶胶态。该浓度阈值可以通过下式获得，
式中，如果左边的累积值小于1，则说明阈值没有被超过，因而没有气溶胶会生成，而当一旦超过阈值时，就假定气态和气溶胶态满足热力学平衡。值得注意的是，这样的阈值浓度只有当假定没有作为吸收物质的气溶胶时才存在，另外根据该式获得的阈值浓度必须假定只存在吸收过程，而没有吸附或其他的气粒转换过程。
纯的低挥发性产物的饱和时浓度的倒数刚好等于转换系数Kom,i(m3ug-1)：
通过推理可以计算得到SOA的质量产率Y，
(二)有机气溶胶浓度查算表法：
利用SORGAM模型分析影响SOA生成的重要因子，包括反应性有机气体的消耗量、初始有机气溶胶浓度、温度和气压。对这些因子进行分档，如表2所示，共计3686400种情形。因子的分档主要根据中国地区对流层低层的实际情况划分，其中气压下限取500hPa是由于SOA主要集中在6km以下的大气中，其环境气压一般高于500hPa。基于表2，利用SORGAM模型计算不同因子组合情况下SOA的生成率，最终建立了SOA查算表数据库。
表2SOA查算表分档设置
实施例应用
CALGRID中分别应用耦合SORGAM模式和有机气溶胶查算表，选取的模拟时段为2006年全年。选取的模拟区域(21.8N～24N 112.16E～114.82E)包括广东省广州、深圳、佛山、珠海、东莞、中山、惠州、清远、肇庆、江门等10个地区和香港、澳门两个特别行政区。气象场由中尺度气象模式WRF提供，共有四层嵌套，最内层的嵌套区域是空气质量模式的模拟区域。具体的网格设置见图1。所使用的源排放清单是Zhang Qiang和D.G.Streets等的2006年东亚地区源排放资料，空间分辨率为0.5°。
图2给出了模拟的含碳颗粒物的地面浓度分布。可见浓度较大的区域在广州佛山一带，大部分地区的含有机碳的颗粒物的浓度超过2.5μg/m3左右；最大值在广州附近，大值区浓度超过3.5μg/m3。与2006年珠江三角洲空气质量观测报告中的PM10的分布相似。CALGRID模型可以模拟常规一次污染物(SO2、NO2、PM10、PM2.5等)的浓度，可以模拟二次污染物(O3等)，硫酸盐(SO42-)、硝酸盐(NO3-)、铵盐(NH4+)、黑碳(BC)、有机碳(OC)等颗粒物，以及S、N沉降量，具有高分辨率、高时效性，能够方便快捷地给出源和受体之间的响应关系，结果可靠。
图3给出了通过2种方法的计算结果。可见，当SORGAM算法的SOA浓度小于10μg·m-3时，查表法和耦合法的计算结果总体上较为吻合；而当浓度高于10μg·m-3时，两种方法的模拟结果有一定差异。
造成这种结果的原因与查表法分档有关。在实际大气中人为源SOA的浓度一般不超过10μg·m-3，在此范围内本研究建立的查表法同耦合法的模拟结果十分接近。因此本研究建立的SOA查表法也能够反映SOA及其前体物的关系，基本具备SORGAM气溶胶模式的模拟能力。