基于三维元胞自动机的垃圾焚烧烟气扩散模拟

顾莎莎，王浩华

(海南大学 信息科学技术学院，海南 海口570228)

随着社会经济的不断发展，城市规模、人口的急剧增加，对城市的公共服务和社会管理提出了严峻挑战.数据显示，目前全国三分之二以上的城市面临“垃圾围城”问题，垃圾堆放累计侵占了土地50 000 ha.与此同时，我国的垃圾无害处理设施匮乏，尽管经过一系列处理，但垃圾对环境的污染依然比较严重.

目前垃圾处理的最好方法是进行焚烧，既可避免垃圾填埋侵占大量的土地，又可利用垃圾焚烧产生的能量进行发电和供暖，获得可观的经济效益.但由于建设投资规模、运行监管力度不到位，焚烧厂的建立多不成规模，在污染治理方面的投入也会受到影响，致使其污染物排放比较严重，难以达到国家的排放标准，对环境的危害也较大.此外，焚烧厂的建立往往伴随着争论、纠纷甚至冲突，例如2007 年北京海淀区六里屯焚烧厂事件、2009 年广州番禺垃圾发电厂事件等等.为了监控垃圾焚烧厂污染物对环境造成的影响，不仅在焚烧厂内进行监控，还有必要对焚烧厂的周边环境进行动态监控，从而形成为民众所信服的全方位垃圾焚烧厂环境监控体系，也为主管部门制定管理应变决策提供科学依据.

对于气体扩散的研究，大多采用风洞测试、现场试验以及数值模拟(如高斯烟羽模型)的方法.考虑到烟气的扩散具有时空特性，笔者从数值模拟角度出发，采用三维的元胞自动机方法建立模拟污染物扩散的CA 模型，并将风向、风速以及降雨等影响因子纳入考量，由此动态模拟焚烧厂烟气的扩散过程，实现对焚烧厂周围环境的动态监控.

1 三维元胞自动机模型

元胞自动机模型(Cellular Automata Model，CA 模型)，是一个时间和空间都离散的动力系统;通过抽象化来定义一个系统的状态，且没有一般的数学物理方法描述或仿真系统时的死板［1］，但其状态改变的规则在时间和空间上都是局部的.CA 模型被广泛应用于计算机图形学、生命系统的自我复制、森林火焰蔓延、城市交通流、城市空间变化和人员疏散等多个领域，是复杂性科学研究中相当活跃的前沿领域［2］.

CA 模型主要由元胞、状态、邻居以及规则4 个部分组成，其主要思想是将复杂系统分解为多个元胞，且在系统中各相邻元胞之间有一定的相互作用关系，通过制定一些简单的规则，大量元胞在这些规则下进行更新，不仅可以简化分析过程，而且可以得到较好的仿真效果(趋于稳态、逐渐膨胀、逐渐收缩和周期变化)，具有良好的稳定性与可靠性［3］.

具体来说，CA 模型即为A=(Ld，S，N，f)，其中，A 代表一个元胞自动机系统;Ld表示元胞空间，d 为正整数，表示元胞自动机内元胞空间的维数;S 是元胞的有限的、离散的状态集合;N 表示元胞邻居;f 表示状态转移函数(规则).

1.1 模型的建立 考虑烟气污染物的扩散是向空间的各个方向进行，笔者提出了三维烟气扩散演化的元胞自动机模型.

将焚烧厂烟气排放处周围离散成N×N×M 的长方体网格，每个网格代表一个元胞，则所有网格组成了一个三维元胞空间.取任一元胞其邻近的10 个元胞作为邻居，具体包括8 个水平面上的元胞和2 个垂直面上的元胞(如图1 所示)，则下一时刻元胞的状态只决定于当前时刻元胞及其邻居的状态.此外，在水平方以及垂直向上方向上的边界可以认为只出不进，而垂直向下方向由于当污染达到时处于接收状态，因此垂直下方边界采取反弹边界.

图1 三维元胞邻居示意

对于具体规则的确定，只考虑重力效应以及平流输送和湍流扩散对污染气体的影响，不考虑污染物之间的化学反应、沉降等作用，因此对于t+1 时刻用t 时刻的状态可描述为其中表示t 时刻元胞ijk 的污染物浓度表示水平上8 个元胞在平流作用下对元胞ijk 的影响;表示垂直方向上2 个元胞在大气湍流的作用下对元胞ijk 的影响.而10 个邻居对元胞ijk 的影响则利用加权的形式确定.对于中心元胞的相邻单元，污染物传递系数为对于非相邻单元(水平面的4 个斜角方向)，其传递量应小于相邻单元，故利用斜角系数d 作为补正，有故有:由于元胞流出量不能大于原来该元胞的污染总量，则有:wa+wad≤0.25，由于在无风、无流的时候实际扩散形状最接近圆形，而由Karafyllidis 的研究［18］得出wa=0.084，d=0.16.

考虑到对流传输对元胞ijk 的影响来自8 个方向，因而对流传输中元胞ijk 在t+1 时刻的浓度变化为其中表示风速度在上一时间步t 时刻沿北、南、东、西、东北、西北、东南、西南方向风速的影响程度.

湍流、重气效应以及降雨量对元胞ijk 的影响主要表现在垂直方向上，有其中tk'z，tk″z为垂直向下与向上元胞在湍流、初始速度及降雨量作用下的综合系数.而降雨也对污染物有吸附作用，故需在上述模型中加入吸附项，有其中为吸附系数.

综合上述有下列规则

1.2 参数的确定 风作用于污染烟气使得烟气随风漂移，其漂移速率同风速成一定的比例关系，类似于水体污染带的扩散漂移且当风向确定后其他水平方向上有twn=0.

垂直方向上的扩散主要由湍流效应与雨量共同确定，其中湍流对烟气的扩散有上下2 个方向的影响，而雨量只产生垂直向下的影响，当雨量的作用不明显(如没有降雨)时，湍流将占向下影响的主导地位.

垂直方向湍流扩散系数tkz可通过由大气稳定度得到的垂直扩散系数σz［5］确定，本文取向上因子tk″z=min σz，向下因子的确定同时考虑湍流与降雨的影响.不同的雨速对烟气的稀释作用也不同，文献［6 －7］指出了雨滴半径与雨速、降雨量的关系，经过整理可得到小雨、中雨、大雨及暴雨4 个等级下的雨速，如下表1所示.

表1 降雨强度与雨速的关系

由此确定向下的降雨因子tkzy' =tvzy/ vzmax，其中vzmax表示雨速的最大值.当没有有降雨时湍流效应起主导地位，因此有

2 元胞自动机实现对环境影响的监控

2.1 元胞空间与尺度 考虑污染物扩散的范围，以烟囱的出口点为中心将元胞空间划分为120 m(排烟口高80 m，由文献［8］可知，烟气抬升高度平均约40 m)，长和宽分别为4 000 m 的长方体.

此外，元胞作为CA 模型最基本的空间单元，其不同的空间尺度和时间尺度也会影响到CA 模型的模拟结果.

文献［9］阐明了随着元胞空间尺度的增大，孤立元胞将增多且在运行中很难发生演化，故会使得CA模型的模拟精度下降.事实上，元胞空间尺度越精细，模拟结果的精度越高，运行所花费的时间越长.根据文献［9］中模拟得到的结果，可以发现相对于100 m×100 m、150 m ×150 m、200 m ×200 m 的元胞尺度，50 m×50 m 的精度达到了最高.本文同样采用50 m×50 m 的尺度，而在竖直方向上由于扩散较慢，取每一元胞高10 m.对于时间尺度，直接利用根据已知数据确定迭代步长对应的真实时间，对一年内的扩散进行迭代，每天以小时更新.

2.2 计算说明 利用前文建立的CA 模型来描述水平与垂直方向上污染物浓度的变化，其具体过程如图2 所示.

以深圳市宝安区一垃圾焚烧厂污染烟气作为算例，根据深圳市气象局网站［10］可得到2011 年4 月～2012 年3 月的气象资料(风速、风向和雨量等，如图3 和图4 所示).

图2 元胞自动机模拟流程图

图3 风速、风向频率玫瑰图

图4 宝安区2011 年4 月~2012 年3 月降雨量

计算时，主要考虑颗粒物、HCL、SO2、NOX，重金属和二噁英等污染物的影响.且根据文献［11 －12］可知不同降雨对PM10，PM2.5 的稀释作用整理后如表2 所示(其中二噁英按重金属处理).

表2 不同降雨量下各污染物吸附因子

2.3 模拟结果 根据上述演化规则和条件，利用Matlab 编程可得各污染物浓度的分布图，其中部分图如图5 所示.

图5 元胞自动机下部分污染物浓度分布

再根据熵权法确定污染度，可得到污染度分布，最后取适当的阀值对周围区域划分污染度划分，其结果如图6 和7 所示.

图6 污染度分布图

图7 CA 模型污染度等级分布图

根据上述结果，可以发现各污染物浓度分布大致相同，只存在一些小的差别，这是由于各污染物本身扩散的差别引起的.由此也得到了2011 年4 月～2012 年3 月焚烧厂周围污染严重程度的分布，受污染的主要方向为西南方，与这一年内的主要风向是相符合的，且在扩散到一定程度时，污染程度会降低到正常水平.

对本文CA 模型在不同尺度下(未考虑改变竖直尺度)的效率进行了简单测试，结果如图8 所示.可以发现，在相同迭代次数下，运行时间随着元胞空间尺度的减小而呈指数型增大，且本文选择50 m×50 m×10 m 的元胞也是比较合适的.

图8 元胞尺度与运行时间的关系

3 结束语

由于元胞自动机在时间和空间上的离散作用，只需确定简单的规则就可以模拟复杂的系统. 本文利用元胞自动机对污染烟气的扩散进行建模，通过风速、风向、降雨等因素的影响来确定了元胞的演化规则，从而使元胞自动机模型更加有效地模拟了污染烟气的扩散，这为焚烧厂周围制定污染补偿计划等政策都提供了有效的依据.

另外，在模拟烟气的扩散时，本文只重点考虑了影响烟气扩散的主要因素，而实际中还有影响扩散的其他不确定因素，如不同种类污染物在扩散过程中的化学反应、土壤对污染物的吸收等等，这都需在今后做进一步的研究.

［1］Li X . A constrained CA model for the simulation and planning of sustainable urban forms by using GIS［J］.Environment and Planning B，2001，28(9):733 －753.

［2］Crowl D A，Louvar J F. Chemical Process Safety:Fundamentals with Applications［M］.NewJersy:Prentice Hall PTR，2002.

［3］王安麟，王炬香.基于CA 的结构拓扑设计的研究—平面薄板加强筋拓扑形态创成的仿真［J］.机械强度，2001，23(2):181 －186.

［4］王璐，谢能刚，李锐，等.基于元胞自动机的水体污染带扩散漂移仿真［J］.水利学报，2009(4):481 －485.

［5］吴琦. 环境风险评价中危险气体扩散模式计算方法研究［D］.哈尔滨:东北林业大学，2008.

［6］盛裴轩. 大气物理学［M］.北京:北京大学出版社，2003.

［7］瑞语. 计量雨滴［J］.中国计量，2009(7):52.

［8］谷清. 烟气抬升公式计算对比［J］.环境科学研究，1991(3):25 －32.

［9］柯新利，邓祥征，何书金. 地理元胞自动机模型的尺度敏感性及原因［J］.地理研究，2010，29(5):863 －872.

［10］深圳市气象局.气象资料.［2014 －07 －17］.http://www.szmb.gov.cn/article/QiHouYeWu/ .

［11］岳毅，马千里，赵伟. 液态降水对PM10 质量浓度的影响分析［C］. //中国气象学会.第28 届中国气象学会年会论文集.出版地不详，中国气象学会，2011.

［12］李芳，张承中. 西安市大气颗粒物PM_(2.5)与降水关系的探讨［J］. 中国环境监测，2013(6):22 －28.