本发明涉及污染物溯源技术领域，尤其涉及一种面向工业园区的大气污染物溯源方法。
背景技术：
大气污染物国内外对于溯源方法有很多的研究，主要集中在核扩散方面，大气污染物溯源一般是用于大范围区域内的常规污染物(尤其是PM2.5)的输送和来源分析，对于较小区域内大气污染物溯源分析的研究较少。Rao SK(2007)Source estimation methods for atmospheric dispersion.Atmos Environ 41(33):6964–6973.doi:10.1016/j.atmosenv.2007.04.064(大气扩散的源估计方法)中，将溯源算法分为前向与后向两大类。前向方法通常使用贝叶斯推断来估计扩散源，后向算法直接求解正向扩散模型的方程来得到结果。因此实际扩散中的不确定性，后向算法往往无法得到准确的结果。贝叶斯推断将求问题的精确解转化为求出解的后验概率分布，许多研究者在此基础上提出了相关的溯源算法。Huber MF(2014)On-line dispersion source estimation using adaptive gaussian mixture filter.IFAC(使用适应性高斯混合滤波器在线估计扩散源)中，使用适应高斯混合滤波器(AGMF)来求解在线贝叶斯问题。不仅限于此，其他的常用的滤波技术也能解决后验概率估计的问题，例如，Zhang XL,Su GF,Yuan HY,Chen JG,Huang QY(2014)Modified ensemble Kalman filter for nuclear accident atmospheric dispersion:prediction improved and sourceestimated.J Hazard Mater 280:143–155.doi:10.1016/j.jhazmat.2014.07.064(用于核事故大气扩散分析的改进的集合卡尔曼滤波器：改进预测结果并估计扩散源)中的卡尔曼滤波，Chin TM,Mariano AJ(2010)A particle filter for inverse lagrangian prediction problems.J Atmos Ocean Technol 27(2):371–384.doi:10.1175/2009JTECHO675.1(用于拉格朗日反问题的粒子滤波器)和Zhang Y,Wang L(2013)Particle filtering method for source localizationin wireless sensor network.advanced technology in teaching:selected papers from the 2012international conference on teaching and computational science(ICTCS 2012):517–523.doi:10.1007/978-3-642-29458-7_76(用于无线传感器网源定位的粒子滤波方法)中的粒子滤波。除了使用滤波技术，Hirst B,Jonathan P,Gonza′lez del Cueto F,Randell D,Kosut O(2013)Locating and quantifying gas emission sources using remotely obtained concentration data.Atmos Environ 74:141–158.doi:10.1016/j.atmosenv.2013.03.044(使用遥控获取的浓度数据对气体排放源进行定位并量化)提出了一种使用马尔可夫链-蒙特卡洛方法来估计多源的算法，这种算法具有较高的理论价值。Ma D,Wang S,Zhang Z(2014)Hybrid algorithm of minimum relative entropy-particle swarm optimization with adjustment parameters for gas source term identification in atmosphere.Atmos Environ 94:637–646.doi:10.1016/j.atmosenv.2014.05.034(用于大气扩散释放源参数辨识的MRE-PSO混合方法)使用了粒子群优化方法并结合高斯模型估计了扩散源的源项。So W,Koo J,Shin D,Yoon ES(2010)The estimation of hazardous gas release rate using optical sensor and neural network.Comput Aided Chem Eng 28:199–204.doi:10.1016/S1570-7946(10)28034-3(使用光学传感器和神经网络估计有害气体泄漏速率)中神经网络也可应用于扩散源的估计。Wang B,Chen B,Zhao J(2015)The real-time estimation of hazardous gas dispersion by the integration of gas detectors,neural network and gas dispersion models.J Hazard Mater 300:433–442.doi:10.1016/j.jhazmat.2015.07.028(基于气体传感器、神经网络和大气扩散模型的有害气体扩散实时预测方法)在神经网络的基础上提出了略过源项直接估计扩散趋势的方法。以上方法能够同时估计出扩散源的位置和释放速率，但是它们估计得到的释放速率都是一个常数，不能很好地反映释放速率不断变化的情况，而这种情况在实际中是很普遍的。
也有很多学者对可变释放速率的溯源问题进行了研究。Chino M,Nakayama H,Nagai H,Terada H,Katata G,Yamazawa H(2011)Preliminary Estimation of Release Amounts of 131I and 137Cs Accidentally Discharged from the Fukushima Daiichi Nuclear Power Plant into the Atmosphere.J Nucl Sci Technol 48(7):1129–1134.doi:10.1080/18811248.2011.9711799(福岛核电站事故泄漏的131I与137Cs总量初步估计)直接通过比例关系对福岛事件后131I和137Cs两种物质的释放速率和总量进行了分析。SchoM,Plastino W,Povinec PP,Wotawa G,Bella F,Budano A,De Vincenzi M,Ruggieri F(2012)Estimation of the time-dependent radioactive source-term from the Fukushima nuclear power plant accident using atmospheric transport modelling.J Environ Radioact 114:10–14.doi:10.1016/j.jenvrad.2011.11.008(使用大气传输模型估计福岛核事故中时变的放射性源项)使用类似的方法通过三个检测站测得的放射性活度数据分析了福岛事件后的核物质扩散情况及释放速率。Chai T,Draxler R,Stein A(2015)Source term estimation using air concentration measurements and a Lagrangian dispersion model–Experiments with pseudo and real cesium-137observations from the Fukushima nuclear accident.Atmos Environ 106:241–251.doi:10.1016/j.atmosenv.2015.01.070(使用大气测量数据与拉格朗日扩散模型估计源项：以福岛核事故中137Cs的仿真和真实数据为例)设计了一个代价函数，通过最小化代价函数得到释放速率的解。但是以上方法都是对一个位置已知的扩散源进行释放速率估计。
技术实现要素：
鉴于目前存在的上述不足，本发明提供一种面向工业园区的大气污染物溯源方法，能够利用扩散源园区周围监测站采集到浓度数据估算出源强大小和扩散源的具体位置，并且能够适应源强和风场数据可变的情景。
为达到上述目的，本发明的实施例采用如下技术方案：
一种面向工业园区的大气污染物溯源方法，所述面向工业园区的大气污染物溯源方法包括以下步骤：
通过监测点采集观测数据并进行超标判断；
根据超标的观测数据进行区域筛选；
获取筛选得到的区域内采集的观测数据根据源位置估计模型进行模型反算；
得到结果并进行验证。
依照本发明的一个方面，所述源位置估计模型为高斯多烟团模型，具体考虑雨水冲洗、干湿沉降以及变化风场的高斯烟团模型如下：
式中，q表示污染物释放速率,(xc,yc,zc)表示当前烟团的中心位置,σx,σy和σz分别表示x，y和z方向上的扩散系数。C(x,y,z,t)计算了t时刻该烟团在(x,y,z)产生的浓度值。。
依照本发明的一个方面，利用函数f来表示式(1)，则扩散方程表示为：
式(2)计算的是一个扩散源位于位置向量l处，瞬时释放量为qj，从tstart时刻开始扩散的烟团在ti时刻位于位置向量zi的观测点的浓度Ci。式中，qj为第j个烟团的释放量，也即烟羽的释放速率，n为烟团个数，δ为各个烟团的释放间隔时间。其中，各个烟团的释放时间间隔满足以下关系：
依照本发明的一个方面，按照公式(2)扩散函数可以使用矩阵表示，具体如下：
C≈G·q (4)
式中，G是已知的传递方程，包含的是正向高斯扩散模型计算信息；q是模型的未知参数向量；C包含了监测站的观测信息。用合理的估计值来代替q，满足上述方程，实现溯源最终的目标。
依照本发明的一个方面，所述根据超标的观测数据进行区域筛选包括：根据风速风向筛选区域和根据地块筛选区域。
依照本发明的一个方面，所述根据风速风向筛选区域具体为：研究区域设置为矩形，然后根据风速和风向确定研究区域范围；所述根据地块筛选区域具体为：根据风速风向确定计算区域后，需要根据企业地块信息进一步减小搜索的范围。
依照本发明的一个方面，所述观测数据为浓度数据根据源位置估计模型进行模型反算具体包括：根据观测数据进行源项估计，包括估计扩散源位置l和泄露速率q(t)。
依照本发明的一个方面，所述源项估计具体包括：估计污染源的位置及污染物释放速率。
数据d来自于位置集为z、风场向量为W的监测站，风场W随空间和时间变化。假设源项参数为θ＝{l,q(t)}的危害物质污染源在时间tstart开始释放，于时间tend停止泄漏。在时间t，位置z的浓度的理论测量值cz(t)可以通过上节的扩散方程cz(t)＝f(t,z,l,q,tstart,tend)计算得到。此外，我们有m个观测数据集，第ith数据集ωi包含了监测站的位置信息zi，时间信息ti以及观测浓度信息ci(ωi＝{zi,ti,ci|zi∈z})。因此，可以得到下面的公式：
其中，测量误差e＝{ei}且ei～N(0,σe2)。
利用贝叶斯推断方法，得到公式p(θ|d)∝p(d|θ)p(θ)，其中p(θ|d)是后验概率函数，p(θ)是先验概率函数，p(d|θ)是源项参数的似然函数，可以表示成如下公式：
式中σe表示仪器的测量误差。
依照本发明的一个方面，利用岭回归分析方法将最优q(t)先估算出来，用θ'＝{l}代替θ将问题转变为估算p(θ'|d)。根据最大似然估计方法，假设通过计算得到释放速率为qm(t)，似然函数进一步表示为：
p(d|θ')≈p(d|qm(t),θ') (7)
式中，p(d|θ')和p(d|qm(t),θ')分别表示在条件θ′和条件qm(t)∧θ′时产生测量值d的似然概率。。
依照本发明的一个方面，所述面向工业园区的大气污染物溯源方法包括：收集污染物超标排放历史事件、工业园区企业信息以及历史结果的多源数据为先验信息，以供进行结果验证。
本发明实施的优点：利用扩散源园区周围监测站采集到浓度数据估算出源强大小和扩散源的具体位置，并且能够适应源强和风场数据可变的情景；提出的源强估计算法，将一个烟羽划分多个烟团，使问题最终转变为多元线性回归问题，并采用岭回归分析方法进行了释放速率的估计，在提高计算效率的同时，保证了算法的精确度和稳定性，同时避免了最小二乘方法的多重共线性；对整个区域进行了网格划分，收集了企业排放的历史信息，保存系统每次运算的结果并根据企业的工作确定了不同情况下每个网格的先验概率，提高了计算结果的可信度。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明所述的一种面向工业园区的大气污染物溯源方法示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，一种面向工业园区的大气污染物溯源方法，所述面向工业园区的大气污染物溯源方法包括以下步骤：
步骤S1：通过监测点采集观测数据并进行超标判断；
在实际应用中，发现某污染物超标：通过安装的传感器检测到污染物的超标排放。当某一污染物连续检测到超标的时间超过某一阈值后，从空气质量监测站(检测仪)、传感器中获取所需的污染物和风速风向测量数据，从数据库中获取相关先验数据，启动污染物溯源反算。
步骤S2：根据超标的观测数据进行区域筛选；
所述根据超标的观测数据进行区域筛选包括：根据风速风向筛选区域和根据地块筛选区域。
根据风速风向筛选区域：研究区域设置为矩形，然后根据风速和风向确定研究区域范围。当风速较大时，风速对烟羽的覆盖范围没有明显的影响，因为监测站位于排放源的下风向，所以可以根据风向在监测站的上风向划定一个扇形区域，该区域内的可能出现超标排放源，而该区域外的部分可不参与反算计算。而当风速较小时，部分污染物可能会反向流动，因此风速较小时需要扩大搜索区域。
根据地块筛选区域：根据风速风向确定计算区域后，需要根据企业地块信息进一步减小搜索的范围。因为不同的企业生产不同的化工产品、需要不同的化工原料，所以不同企业地块排放的污染物种类是不同的。根据企业的排放信息，确定每种污染物的排放地块，并写入数据库。反算时通过查询数据库，将不排放待反算污染物的区域剔除，达到减少计算资源，提升计算效率的作用。
步骤S3：获取筛选得到的区域内采集的观测数据根据源位置估计模型进行模型反算；
在实际应用中，输入数据准备：采集多监测点得到的观测数据。采用之前步骤筛选得到的区域作为反算的估计空间。在后验概率计算公式(6)中引入仪器检测误差。并对估计区域进行网格划分，由历史排放信息、企业工作时间和系统历史运行结果等确定每个网格的先验概率。在模型反算中需要获取某点的先验概率时，可以使用该点所在网格的概率值。
在实际应用中，所述源位置估计模型可为高斯多烟团模型，具体考虑雨水冲洗、干湿沉降以及变化风场的高斯烟团模型如下：
利用函数f来表示上式，用于近似计算一个释放速率可变、风场可变的烟羽的浓度值，那么扩散方程还可以表示为：
式(2)计算的是一个扩散源位于位置向量l处，瞬时释放量为qj，从tstart时刻开始扩散的烟团在ti时刻位于位置向量zi的观测点的浓度Ci。式中，qj为第j个烟团的释放量，也即烟羽的释放速率，n为烟团个数，δ为各个烟团的释放间隔时间。其中，各个烟团的释放时间间隔满足一下关系：
当然，按照公式(2)扩散函数还可以使用矩阵表示，如下：
C≈G·q (4)
式中，G是已知的传递方程，包含的是正向高斯扩散模型计算信息；q是模型的未知参数向量(如扩散源位置和释放速率)；C包含了监测站的观测信息。用合理的估计值来代替q，满足上述方程，就是溯源最终的目标。
在实际应用中，所述观测数据为浓度数据根据源位置估计模型进行模型反算具体包括：根据观测数据进行源项估计，包括估计扩散源位置l和泄露速率q(t)。如何利用观测或测量的浓度数据来估计扩散源位置l和泄露速率q(t)，具体包括：
数据d来自于位置集为z、风场向量为W的监测站，风场W随空间和时间变化。假设源项参数为θ＝{l,q(t)}的危害物质污染源在时间tstart开始释放，于时间tend停止泄漏。在时间t，位置z的浓度的理论测量值cz(t)可以通过上节的扩散方程cz(t)＝f(t,z,l,q,tstart,tend)计算得到。此外，我们有m个观测数据集，第ith数据集ωi包含了监测站的位置信息zi，时间信息ti以及观测浓度信息ci(ωi＝{zi,ti,ci|zi∈z})。因此，可以得到下面的公式：
其中，测量误差e＝{ei}且ei～N(0,σe2)。
方程cz(t)＝f(t,z,θ)是一个非线性方程，难以直接计算求解得到源项参数的解析解，我们需要通过贝叶斯推断的后验概率来求解该类反问题。利用贝叶斯推断方法，我们得到公式p(θ|d)∝p(d|θ)p(θ)，其中p(θ|d)是后验概率函数，p(θ)是先验概率函数，p(d|θ)是源项参数的似然函数，可以表示成如下公式：
式中σe表示仪器的测量误差。准确计算后验概率分布函数的关键在于先验概率分布函数以及似然函数。先验概率函数可以从历史数据或经验信息中获得，并且在信息不充分的情况下可以用均匀分布函数U(θ)代替，如果有丰富的工厂历史排放信息，那么可以假设一个工厂超排的先验概率与其历史超排次数成比例关系，充分利用先验信息的作用。而对于似然函数的计算，考虑到释放速率q(t)是一个随时间变化的函数，所以直接估计是非常困难的。因此我们考虑利用岭回归分析方法将最优q(t)先估算出来，那么源项参数的估计就只剩下位置向量，可以用θ'＝{l}代替θ将问题转变为估算p(θ'|d)。根据最大似然估计方法(MLE)，假设通过计算得到释放速率为qm(t)，那么似然函数可以进一步表示为：
p(d|θ')≈p(d|qm(t),θ') (7)
式中，p(d|θ')和p(d|qm(t),θ')分别表示在条件θ′和条件qm(t)∧θ′时产生测量值d的似然概率。。
源项估计的结果实际上得到的是空间分布的后验概率分布函数，那么，一种简单的方式就是将需要估计的实际空间划分成一个个相等的网格，然后对每一个网格计算后验概率分布函数，这样就可以近似得到整个估计空间的PDF。这种方法的计算准确性依赖于网格划分的尺寸，如果网格尺寸大了就容易导致计算结果不准确；相反的，网格划分小了就会造成计算时间过长，不符合应急救援的要求。所以网格的划分也是实际实现时需要考虑的重要问题，另外使用蒙特卡洛-马尔可夫链(MCMC)方法求解后验概率分布函数时，可以对估计空间先进行采样，就会得到后验空间中可能释放源的样本。如果使用粒子群算法，也可以用很快的速度在估计空间中寻找到最优解。使用网格划分法能很好地覆盖整个估计空间，但是它的计算精度会严重影响运算速度，精度越高计算越缓慢；MCMC的随机性较强，理论效果最好，但是实际应用中由于扰动可能会产生比较大的误差；粒子群算法在计算效率和计算精度上都有巨大的优势，但是粒子群算法可能寻找到局部最优解而非全局最优解。提出了三种解决扩散源位置估计问题的算法，分别是网格化方法、MCMC算法和粒子群算法，可分别应用于不同的场景。应急场景下一般使用粒子群算法，常规反算可使用粒子群与网格划分综合的结果，而理论分析时可使用MCMC算法。
步骤S4：得到结果并进行验证。
得到计算结果：如使用粒子群算法，最后得到估计区域中后验概率最高的点，使用网格划分方法可以得到每个网格的后验概率，使用MCMC方法能得到整个估计空间后验概率的粒子分布表达形式。之后根据测量数据验证计算结果的准确性。
在实际应用中，所述面向工业园区的大气污染物溯源方法包括：收集污染物超标排放历史事件、工业园区企业信息以及历史结果的多源数据为先验信息，以供进行结果验证；对历史排放数据、企业工作时间以及过往系统运行结果等多源数据进行融合，为下一次计算提供先验信息。以往的反算算法都假设先验概率为均匀分布，也很少有研究者对先验信息进行收集和分析。本发明首先对整个区域进行了网格划分，收集了企业排放的历史信息，保存系统每次运算的结果并根据企业的工作确定了不同情况下每个网格的先验概率，提高了计算结果的可信度。
本发明实施的优点：利用扩散源园区周围监测站采集到浓度数据估算出源强大小和扩散源的具体位置，并且能够适应源强和风场数据可变的情景；提出的源强估计算法，将一个烟羽划分多个烟团，使问题最终转变为多元线性回归问题，并采用岭回归分析方法进行了释放速率的估计，在提高计算效率的同时，保证了算法的精确度和稳定性，同时避免了最小二乘方法的多重共线性；对整个区域进行了网格划分，收集了企业排放的历史信息，保存系统每次运算的结果并根据企业的工作确定了不同情况下每个网格的先验概率，提高了计算结果的可信度。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本领域技术的技术人员在本发明公开的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。