近年来,石化工业园区发展迅速,带来巨大的经济效益,但是污染问题也成为政府部门棘手的工作,有研究表明工业源是挥发性有机化合物(volatile organic compounds,VOCs)贡献的第一大源[1-2]。VOCs对大气环境和人体健康产生很大的危害,日益受到人们的关注,并成为国际研究的热点[3]。
对于监测难度较大的无组织排放,安装在园区内部的分布式监测方法是一种行之有效的方法[4]。该方法能实时监测园区内部的VOCs浓度情况,并且在发现浓度超标时进行报警并反演出泄漏源,如何基于报警的监测点锁定泄漏源是分布式监测网络的关键目标。虽然针对环境中VOCs溯源的研究有很多[5-9],但这些研究大都针对中尺度的溯源,基于分布式监测网络实现园区内部件级泄漏溯源的研究很少。由于VOCs无组织泄漏属于微量泄漏,其初始动量比较小,因此风场对扩散路径影响很大。但是由于园区内遍布工艺装置,因此泄漏点处局部风场同自然来流并不一致,在某些区域流动很通畅,而在有些设施密集的区域流动很混乱。研究园区内局部风场的特征,对于VOCs的溯源非常重要。
本文以大庆炼化芳烃抽提设备为对象,建立了基于质谱仪的VOCs无组织排放在线监测网络,采用计算流体力学仿真(computational fluid dynamics,CFD)对自然风下的设备区内流场进行了模拟,基于小波熵理论,对监测点局部流场的不稳定性进行了研究,并对多种流场参数进行了关联性分析。
1 系统介绍及方法 1.1 VOCs在线监测网络图 1为大庆炼化芳烃抽提装置的VOCs分布式监测网络示意图。该芳烃抽提设备区包括2个平台、2个泵房以及地面工艺单元,占地面积为100 m×60 m。图 2为某些工艺部件上标记的潜在泄漏点,可以发现单一设备上潜在泄漏点分布密集且主要存在于阀门、法兰以及密封连接件上,图 3则为红外摄像仪所捕捉到的图 2中泄漏源正在发生的无组织泄漏。为了对设备内的无组织泄漏进行实时监测,在生产装置内部多处布设60个监测点位,各监测点位负压抽取附近环境空气经过除尘除水后进入质谱仪分析。在线环境质谱仪放置于设备旁边的分析小屋,同时在装置高点放置气象参数传感器采集风速、风向和温湿度等气象参数。VOCs监测系统会实时分析各测点区域戊烷、己烷、庚烷、苯、甲苯、二甲苯和乙苯7种介质的组分情况,当发现浓度超标时,则需要进行泄漏报警并尝试推演出发生泄漏的区域。
由于无组织泄漏并非喷射流,而是以微量的羽流形式持续泄漏,其在外界的扩散路径受该处局部风场的影响极大,因此需要获得各监测点局部风场的信息。CFD是研究建筑群内流动的一种有效方法[10-11]。利用CFD获得监测点在不同来流下的局部风场,以作为溯源分析的依据。由于设备区内存在大量的工艺单元,分布在室内泵房、室外地面和平台等区域,且各区域的密集程度也不一致,障碍物分布极不规则。自然风进入设备区后产生的流场必然是不规则且不稳定的,而不同位置的监测点,其流场的不稳定性程度也不一致。例如,对于流动畅通的区域,其风场的波动较小,CFD获得的风场信息可以作为泄漏溯源的依据;若监测点所在区域对流不强或处于旋涡区,则可能无风或风向剧烈摆动, 则可能无法基于该点的风场进行溯源。因此需要对园区内局部风场的不稳定性进行研究,为溯源分析奠定基础。
1.2 CFD仿真模型计算域如图 4所示,对芳烃抽提设备采用激光扫描获得其点云图像,并进行了详细的三维建模,同时对设备周边大型建筑也进行了建模以考虑其对设备区风场的影响。采用Ansys workbench cutcell方法对计算域进行网格切分,网格数总量为1 400万。计算域为非稳态不可压缩流动,湍流模型为k-ε,压力耦合求解采用SIMPLE算法。计算域除了地面外的表面设为统一的速度入口,速度边界采用指数率风廓线[12]:
$ u = {u_0}{\left( {\frac{z}{{{z_0}}}} \right)^{0.22}}. $ | (1) |
其中:u为z高度下的风速,单位m/s;u0为z0高度下的风速,单位m/s。由于气象参数传感器安装的高度为18 m,因此z0=18。选择城市地区,幂指数为0.22。
18 m处风速设为1 m/s,大气来流风向为160°,流动模拟时间为3 600 s,时间步长为0.5 s,采用商业软件fluent求解。在仿真时,监测并输出50个位于室外的监测点处水平方向的速度分量。
1.3 小波理论和小波熵方法小波变换能将信号分解为各个不同尺度下的一系列基小波的组合,使得小波变换在时间和频率上具有变化的分辨率,即信号较高频率部分具有较低的时域分辨率,而在低频部分具有较高的时域分辨率[13-14]。本文使用二进离散小波变换的快速方法,利用正交基小波将信号分解成不同尺度下的各个分量。原始信号序列S可以表示为各分量之和
$ S = {a_n} + \sum\limits_{i = 1}^n {{d_i}} . $ | (2) |
其中,an为小波近似解a,di为各级细节解d,n为小波分解级数。
为了统一,将an表示为dn+1,则有
$ S = \sum\limits_{i = 1}^{n + 1} {{d_i}} . $ | (3) |
熵是用来描述系统复杂或混乱程度的一个物理量。最早来源于热力学熵,由Clausius[15]于1865年提出,用于阐明热力学第二定律。1948年Shannon[16]将熵的概念引入到信息论,定义了信息熵(information entropy)作为随机事件不确定性的量度,信息熵的增加意味着系统信息量的减少和复杂程度的增加。小波熵(wavelet entropy)的概念最早由Blanco等[17]基于小波变换提出,并用于ERPs(event related potentials)的分析。Rosso等[18-21]将其用于分析EEG (electro encephalogram)等非平稳信号的检测与分类,并取得了较好的效果。
小波变换是在不同时频域尺度上对信号的分解,而这一过程也将信号的能量划分开来。定义某一尺度下的小波分量的能量为该尺度下小波系数的平方和:
$ {E_i} = \sum\limits_{k = 1}^N {|{d_i}\left( k \right){|^2}} . $ | (4) |
其中,k为时间序列,k=1, 2, …, N,N为信号的采样点数。
各尺度分量的能量之和就是信号总能量,
$ E = \sum\limits_{i = 1}^{n + 1} {{E_i}} . $ | (5) |
归一化处理得到原始信号能量在不同小波尺度下的分布,
$ {p_i} = \frac{{{E_i}}}{E}. $ | (6) |
小波熵定义为
$ {H_E} = - \sum\limits_{i = 1}^{n + 1} {{P_i}{{\log }_2}{P_i}} . $ | (7) |
小波熵是对信号复杂程度和不确定性的一种度量,是衡量信号能量在各小波尺度上分布情况的物理量。小波熵值与小波变换系数的分布有关,与其值大小本身无关,也就是说小波熵反映的是信号各尺度之间的不变性(确定性)。小波熵越大,代表信号中出现的变化越多,信号越趋向于复杂不确定;反之,小波熵越小,说明信号变化越规律,周期性越明显。
1.4 关联性分析Pearson相关系数方法是一种统计学方法,可以定量地衡量变量之间的相关关系。2个变量之间的Pearson相关系数定义为2个变量之间的协方差和标准差的商[22]:
$ r\left( {X, Y} \right) = \frac{{\sum\limits_{i = 1}^n {\left( {{X_i} - \overline X } \right)\left( {{Y_i} - \overline Y } \right)} }}{{\sqrt {\sum\limits_{i = 1}^n {{{\left( {{X_i} - \overline X } \right)}^2}} } \sqrt {\sum\limits_{i = 1}^n {{{\left( {{Y_i} - \overline Y } \right)}^2}} } }}. $ | (8) |
Pearson相关系数的取值范围为-1≤r≤1。相关系数的绝对值越大,相关性越强;相关系数越接近于0,相关性越弱。
2 结果与讨论 2.1 风向小波熵的讨论对于CFD测得的监测点的风向,为了消除流场初始化的影响,在分析每个监测点时提取该点1 000~3 600 s的数据,并计算风向小波熵HD。由图 5可以发现各监测点的风向小波熵是不一致的。按照小波熵由小到大的顺序选取其中的32号、1号、3号、25号、34号、38号、4号,并绘制其风向Wd和风速Wv信号(图 6),可以发现各个点的风场均不稳定且一直变化。随着小波熵逐渐变大,风场的不稳定程度也随之增大。如小波熵最小的32号监测点,其风向一直维持在182°,而小波熵最大的4号监测点,其风向在0~360°范围内都有分布。这7个测点的风速和风向方差如图 7所示,可以发现小波熵小的点,其风场变化比较小,相反小波熵较大的点,其风场的变化非常大,尤其是风向方差的变化同风向小波熵表现出较好的一致性。
可以发现小波熵同监测点局部风场的不稳定性具有良好的相关性,风场的不稳定性会导致VOCs扩散的不确定性。小波熵可以反映风场的不稳定性,这是因为小波熵能揭示风场信号中的能量分布。如果一股正常气流的通道上没有障碍物,则其流动很畅通且不会改变。但由于园区存在很多形状和分布不规则的障碍物,气流流经障碍物时会发生绕流和附面层分离。绕流会改变流动的方向,在附面层脱落前,其流动是稳定的。而附面层脱落则会产生旋涡,并在障碍物后方流场产生尾涡区,旋涡会破坏流场的稳定性使流动非定常,典型的如圆柱绕流产生的Karman涡街。对于原本流动稳定的一点,如果附近出现了旋涡,则旋涡对该点的流动会施加作用并改变其流动状态,如流向上产生波动,这种作用带来了流动能量上的改变。小波熵能将信号在不同能量上进行分级,可以反映各影响因素之间能量相对分布情况,因此通过分析小波熵可以得知该点风场的影响因素。
为了进一步验证小波熵同流动状态的关系,提取了3、25、34和38号监测点处的风场,如图 8所示。其中3号和34号测点位于高层平台,25号和38号点位于地面,3号和25号监测点的小波熵小于34号和38号。3号点(图 8a)位于工艺部件的迎风面,气流在到达此处后沿表面流动,且未发生附面层脱离,因此虽然流向发生了偏转但流动稳定,小波熵较小;34号点(图 8b)处于高层平台,虽然整体遮挡物较少,但是在该风向下,其处于气流流经前方高塔产生的尾涡区中,同时其附近的工艺单元也产生了旋涡,因此多重旋涡作用导致该点的流动不稳定。25号点(图 8c)位于罐区围堰和泵房形成的通道中,虽然该区域位于设备区的中后方,但是前方混乱的气流进入通道后进行了整流而变得稳定;38号点(图 8d)虽然也位于泵房墙壁旁边,但是该区域设备相对密集,存在局部的旋涡,因此流速低且风向不稳定。
2.2 各变量相关性分析
通过上文分析可知,风向小波熵可以表征监测点局部风场的不稳定性,讨论变量与风向小波熵的相关性,还可以分析该变量是否适合做风场不稳定性的评价指标。对监测点风速和风向信号处理后得到流速小波熵Hv,风向方差δd,速度方差δv,风向平均偏转量:
$ \Delta \theta = \frac{{|\overline \theta - {\theta _0}|}}{{180}}. $ | (9) |
其中: θ为该监测点的平均风向,θ0=160°为来流风向。
速度变化量为
$ \Delta V = \frac{{\overline V }}{{{V_0}}} $ | (10) |
其中: V为监测点平均风速,V0=1 m/s为来流风速。
室外50个监测点的变量与风向小波熵的对比如图 9所示,图中也标明了2个变量的相关系数r。由图 9a和9d发现风向小波熵HD与流速小波熵Hv及风向方差δd的变化趋势比较一致,其相关系数(分别为0.736和0.768)也表明具有较强的正相关性。流速和流向是某点流场的2个自由度,并对该点流动唯一定义。当流动受到作用发生改变时,流速和流向都会发生相应的改变,但变化的幅度可能不一致。小波熵并不关注信号的绝对值,而是关注信号能量相对分布。从能量角度而言,流速和风向小波熵展示出来的都是导致流动变化的影响因素,所以二者具有相关性是合理的。对于图 9d风向方差,小波熵表征流场的不稳定性,设备区中存在大量尺度不同的旋涡是不稳定性的重要原因,旋转是旋涡的显著特征,因此熵值越大,流向的变化也会越大。但需要注意的是,虽然风向方差和小波熵变化趋势非常接近,但两者是不同的。方差表示的是风场的变化幅度,不能揭示流场的影响因素,而小波熵从能量分布角度进行分析,能揭示信号发生变化的原因,方差只是这些变化导致的结果。方差和小波熵并不总是一致,如单一正弦波,其方差可能会很大,但是其小波熵为0或者很小,因为只有一种能量分布。
图 9c表明风向小波熵与速度变化量呈现中等程度的负相关,这是因为旋涡的存在是导致园区局部流场不稳定的主要因素,旋涡区一般都为低速区,因此旋涡导致的不稳定性会导致流速减小。这种相关性不如风向方差强,是因为流动变化并不总是导致流速的减小,如附面层脱离前的绕流和渐缩通道会使流速增大,只不过在化工园区,这些过程相对旋涡发生的情形比较少。图 9b和9e表明风向小波熵与风向平均偏转量和风速方差不相关,这是因为流动中存在这样的例子,如附面层脱离前的绕流过程虽然流向随障碍物表面改变,但是流动是稳定的;旋涡虽然流向不稳定,但是风速大都维持在低速区,不稳定的流动可能流速变化不大。
综上所述,风向小波熵可以表征监测点局部风场的不稳定性,流速小波熵及风向方差与其具有强相关性,因此也可以表征局部风场的不稳定性,但小波熵可以揭示导致不稳定性的影响因素,是更接近流动本质的变量。速度变化量、风向平均偏转量和风速方差与风向小波熵关联性较弱或不相关,不适合用来评价局部流场的稳定性。
3 结论本文以大庆炼化芳烃抽提设备为对象,建立了基于质谱仪的VOCs无组织排放在线监测网络,采用CFD模拟获得了160°来流下监测点的风速和风向信号,基于小波熵理论,对监测点局部流场的不稳定性进行研究,并对多种流场参数进行了关联性分析。结果表明:
1) 化工园区监测点的局部流场是不稳定的,且各点的不稳定程度不同。
2) 小波熵能将信号在不同能量尺度上进行分级,揭示流场影响因素的能量相对分布情况,可以表征监测点局部流场的不稳定性。小波熵越大,流动的不稳定性越强。
3) 关联性分析表明风速小波熵以及风向方差与风向小波熵具有强正相关性,可以共同作为风场不稳定性的评价指标。风速变化量同风向小波熵呈现中等程度负相关,而风向平均偏转量和风速方差与不稳定性无关。
在进行溯源时,需要基于监测点处的流场推断能扩散到该处的泄漏源。本文研究表明各点局部流场的稳定性是不一致的。对于流动不稳定的区域,由于流动不能稳定在一个风向,在溯源时无法参照风向来反演泄漏源,只能依据其他原则(如组分和距离)来推断。在未来的工作中,需要基于本文得到的风场不稳定性的评价指标,对各种自然风向下的各区域的流场进行不稳定性评估,制定各区域相应的溯源策略。
[1] |
WANG H L, NIE J, LI J, et al. Characterization and assessment of volatile organic compounds (VOCs) emissions from typical industries[J]. Chinese Science Bulletin, 2013, 58(7): 724-730. |
[2] |
WEI W, WANG S X, CHATANI S, et al. Emission and speciation of non-methane volatile organic compounds from anthropogenic sources in China[J]. Atmospheric Environment, 2008, 42(20): 4976-4988. |
[3] |
DUMANOGLU Y, KARA M, ALTIOK H, et al. Spatial and seasonal variation and source apportionment of volatile organic compounds (VOCs) in a heavily industrialized region[J]. Atmospheric Environment, 2014, 98: 168-178. |
[4] |
HUTCHINSON M, OH H, CHEN W H. A review of source term estimation methods for atmospheric dispersion events using static or mobile sensors[J]. Information Fusion, 2017, 36: 130-148. |
[5] |
JAIN S, SHARMA S K, CHOUDHARY N, et al. Chemical characteristics and source apportionment of PM2.5 using PCA/APCS, UNMIX, and PMF at an urban site of Delhi, India[J]. Environmental Science and Pollution Research, 2017, 24(17): 14637-14656. |
[6] |
BROWN S G, FRANKEL A, HAFNER H R. Source apportionment of VOCs in the Los Angeles area using positive matrix factorization[J]. Atmospheric Environment, 2007, 41(2): 227-237. |
[7] |
WANG X L. Analysis of ambient VOCs levels and potential sources in Windsor[D]. Windsor, Caneda: University of Windsor, 2014.
|
[8] |
SHARAN M, SINGH S K, ISSARTEL J P. Least square data assimilation for identification of the point source emissions[J]. Pure and Applied Geophysics, 2012, 169(3): 483-497. |
[9] |
ALDEN C B, GHOSH S, COBURN S, et al. Bootstrap inversion technique for atmospheric trace gas source detection and quantification using long open-path laser measurements[J]. Atmospheric Measurement Techniques, 2018, 11(3): 1565-1582. |
[10] |
周莉, 席光. 高层建筑群风场的数值分析[J]. 西安交通大学学报, 2001, 35(5): 471-474. ZHOU L, XI G. Numerical analysis of the wind field on high buildings[J]. Journal of Xi'an Jiaotong University, 2001, 35(5): 471-474. (in Chinese) |
[11] |
ZHOU X H, MENG F K, JIANG Y L. Three-dimensional numerical simulation and analysis on wind environment of group buildings[J]. Science Technology and Engineering, 2007, 7(14): 3604-3606. |
[12] |
中华人民共和国住房和城乡建设部.建筑结构荷载规范: GB50009-2012[S].北京: 中国建筑工业出版社, 2012. Ministry of Housing and Urban-Rural Development of the People's Republic of China. Building Structure Load Specification: GB50009-2012[S]. Beijing: China Architecture & Building Press, 2012. (in Chinese) |
[13] |
UNSER M, ALDROUBI A. A review of wavelets in biomedical applications[J]. Proceedings of the IEEE, 1996, 84(4): 626-638. |
[14] |
MALLAT S G. A theory for multiresolution signal decomposition:The wavelet representation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989, 11(7): 674-693. |
[15] |
CLAUSIUS R. On a mechanical theorem applicable to heat[J]. Philosophical Magazine, Series 4, 1870, 40(265): 122-127. |
[16] |
SHANNON C E. A mathematical theory of communication[J]. Bell System Technical Journal, 1948, 27(3): 379-423. |
[17] |
BLANCO S, FIGLIOSA A, QUIAN Q R, et al. Time-frequency analysis of electroencephalogram series(Ⅲ):Information transfer function and wavelets packets[J]. Physical Review E, 1998, 57(1): 932-940. |
[18] |
ROSSO O A, BLANCO S, YORDANOVA J, et al. Wavelet entropy:A new tool for analysis of short duration brain electrical signals[J]. Journal of Neuroscience Methods, 2001, 105(1): 65-75. |
[19] |
QUIROGA R Q, ROSSO O A, BAŞAR E, et al. Wavelet entropy in event-related potentials:A new method shows ordering of EEG oscillations[J]. Biological Cybernetics, 2001, 84(4): 291-299. |
[20] |
ROSSO O A, MAIRAL M L. Characterization of time dynamical evolution of electroencephalographic epileptic records[J]. Physica A:Statistical Mechanics and Its Applications, 2002, 312(3-4): 469-504. |
[21] |
ROSSO O A, MARTIN M T, PLASTINO A. Brain electrical activity analysis using wavelet-based information tools[J]. Physica A:Statistical Mechanics and Its Applications, 2002, 313(3-4): 587-608. |
[22] |
EGGERS J J, BAUML R, TZSCHOPPE R, et al. Scalar costa scheme for information embedding[J]. IEEE Transactions on Signal Processing, 2003, 51(4): 1003-1019. |