本发明公开了一种工业园区大气污染物溯源方法，包括：实时采集工业园区内站点监测数据，包括详细成分站、微型站及气象站，根据监测数据动态进行超标判断；将超标因子与园区企业进行污染物特征因子指纹库进行比对，初步筛选疑似污染源；根据气象数据与园区大小确定污染溯源时间及区间，在溯源追踪时段内，利用拉格朗日模型逐分钟动态计算疑似源对受体点的贡献情况；并结合企业在线工况信息以及企业排口在线监测数据，建立归一化分布滞后模型，进一步锁定污染源；最终统计模拟计算综合评价指数，综合评估园区企业对超标点位的贡献率。本发明能够在不同风向风速及源强排放情况下，实现对工业园区大气污染物的动态精准溯源。
1.一种工业园区大气污染物溯源方法，其特征在于，所述溯源方法包括以下步骤：S1，采集工业园区详细成分站及微型站的污染物监测数据，根据污染物监测数据实时动态进行超标判断；
S2，依据超标判断的结果，得到超标因子与对应的超标点位，将超标因子与所有园区企业的污染物特征因子指纹库比对，初步筛选疑似污染源，将筛选出的疑似污染源的排放贡献因子EF1赋值为1，其余不能匹配的排放源排放贡献因子EF1赋值为0；
S3，采集工业园区气象监测数据，根据气象监测数据与园区内企业分布情况，确定污染溯源的时间T和区间S；
S4，基于拉格朗日数值模型，在溯源追踪时段内，逐分钟动态计算疑似污染源对受体点的浓度贡献情况Mmn，Mmn为不同疑似污染源对受体点的浓度贡献系数；
S5，采集各企业在线工况信息及企业排口在线监测数据，在超标点位建立归一化分布滞后统计模型，计算得到各疑似污染源的排放贡献因子EF2和EF3；其中EF2和EF3分别为基于各企业在线工况信息及企业排口在线监测数据计算得到的长期乘数，表示疑似污染源变动一个单位时，由于滞后效应而形成的对受体污染物浓度总的影响大小；
S6，结合疑似污染源对受体点的浓度贡献情况Mmn和各疑似污染源的排放贡献因子EF1、EF2、EF3，计算得到各疑似污染源对受体点的浓度贡献综合指数Con；
步骤S6中，所述计算得到各疑似污染源对受体点的浓度贡献综合指数Con是指：根据下述公式计算得到各疑似污染源对受体点的浓度贡献综合指数Con：Con＝EF1×EF2×EF3×Mmn。
2.根据权利要求1所述的工业园区大气污染物溯源方法，其特征在于，步骤S1中，所述根据污染物监测数据实时动态进行超标判断包括以下步骤：S11，收集基础污染物监测数据，建立园区大气污染物基本浓度水平库；
S12，根据所有污染物的浓度实时变率和相关浓度标准限值，进行动态超标判断。
3.根据权利要求1所述的工业园区大气污染物溯源方法，其特征在于，步骤S2中，所述初步筛选疑似污染源的过程包括以下步骤：S21，对园区企业进行废气排放梳理，建立污染物特征因子指纹库；
S22，将超标因子与污染物特征因子指纹库进行比对，初步筛选疑似污染源；
S23，针对疑似污染源，定义其排放贡献因子EF1赋值为1，针对其余排放源，维持其排放贡献因子EF1赋值为0。
4.根据权利要求1所述的工业园区大气污染物溯源方法，其特征在于，步骤S3中，所述确定污染溯源的时间及区间的过程包括以下步骤：S31，收集园区基础气象数据；
S32，根据园区企业分布，选取疑似污染源到受体点的最大直线距离为S，根据实时变化的风速 按照如下公式进行积分求解，确定污染溯源的时间T及区间S：
5.根据权利要求1所述的工业园区大气污染物溯源方法，其特征在于，步骤S4中，所述逐分钟动态计算疑似污染源对受体点的浓度贡献情况Mmn的过程包括以下步骤：S41，建立拉格朗日数值模型；
S42，在溯源追踪时段内，在受体点位置释放大量粒子，基于三维空间和时间对疑似污染源输入的变化进行响应函数的计算，得到疑似污染源与受体点之间的传输关系矩阵以及随时间序列变化的网格点气团驻留时间；
6.根据权利要求5所述的工业园区大气污染物溯源方法，其特征在于，步骤S42中，在计算所述疑似污染源与受体点之间的传输关系矩阵时，采用下述公式计算网格点气团驻留时间τi：
式中，Ts为时间分辨率，单位是秒，N为Ts时间范围内采样的数量，J为释放的粒子总数，fijn代表对于指定个点有贡献的粒子数量。
7.根据权利要求5所述的工业园区大气污染物溯源方法，其特征在于，通过下述公式计算得到疑似污染源对受体点的浓度贡献情况Mmn：Ym＝MmnXn
式中，Ym代表模拟得到的浓度值矩阵；Xn为未知量，代表需要反演求解得到的排放源矩阵；Mmn矩阵代表不同网格点排放对受体点的浓度贡献系数矩阵；m代表受体点的个数，n代表疑似污染源的个数。
8.根据权利要求7所述的工业园区大气污染物溯源方法，其特征在于，步骤S5中，所述在超标点位建立归一化分布滞后统计模型计算得到各疑似污染源的排放贡献因子EF2和EF3的过程包括以下步骤：
S51，采集园区内企业的实时在线工况数据以及企业排口在线监测数据，建立园区在线生产排放信息库；
S52,将生产排放信息与超标点位建立归一化分布滞后统计模型进行分析：S521,采用反余切函数转换的非线性归一化方法，将步骤S51中采集到的所有数据进行归一化处理，将数据映射到0至1范围内，形成无量纲量；
S522,采用下述公式建立企业在线工况数据和企业排口在线监测数据与超标点位的归一化分布滞后统计模型：
At＝α+β0Bt+β1Bt‑1+β2Bt‑2+…+βsBt‑s+ut式中：
s为滞后长度，根据步骤S3中确定的溯源时间T进行动态调整；
β0为短期乘数，表示变量B变动一个单位时，对值A的平均影响大小；
βi为延迟乘数，表示过去各时期变量B变动一个单位对A值的平均影响大小；
为长期乘数，表示变量B变动一个单位时，由于滞后效应而形成的对A值总的影响大小；α为常量，ut为误差，t为当前时刻；
S523,结合归一化分布滞后统计模型，计算得到各排放点位的排放贡献因子EF2和EF3。
一种工业园区大气污染物溯源方法
技术领域
[0001] 本发明涉及污染物溯源技术领域，具体而言涉及一种工业园区大气污染物溯源方法。
背景技术
[0002] 随着经济的快速发展和城市化水平的不断提高，大气污染问题也日益突出，对人们的日常生产生活及身体健康具有较大的影响，而大气污染物溯源城市中的工业园区作为
主要的工业源排放，给当地的空气质量达标带来较大压力，也受到周边公众的日益关注甚
至是投诉。这时，一种综合利用园区监测数据及模型统计等方法的，针对工业园区大气污染
物的动态精准溯源是十分有必要的，不仅能够管理园区清洁安全生产服务，更能够改善当
地环境空气质量。
[0003] 目前主流的空气污染溯源方式主要有：方式一、利用稳定同位素释放，来实现污染物的溯源与示踪，该方法成本较高，且耗时较长；方式二、利用空气质量溯源模型，该方法包
含复杂的物理化学反应过程，且空间尺度较大，分辨率较低，相应计算量大，不适用于工业
园区层面；方式三、利用高斯反推模型，推断污染源位置，该方法并不能给出不同污染源对
受体点贡献比率。
[0004] 例如，专利号为CN201611041286.3的发明专利中提及一种面向工业园区的大气污染物溯源方法，包括以下步骤：通过监测点采集观测数据并进行超标判断；根据超标的观测
数据进行区域筛选；获取筛选得到的区域内采集的观测数据根据源位置估计模型进行模型
反算；得到结果并进行验证；利用扩散源园区周围监测站采集到浓度数据估算出源强大小
和扩散源的具体位置，并且能够适应源强和风场数据可变的情景。然而在该发明中，只考虑
了风速风向和企业分布因素，并未考虑到超标因子以及企业在线工况信息及企业排口在线
监测数据，另外，结合先验信息的基础上利用岭回归估计出最大后验概率的解，收集污染物
超标排放历史事件、工业园区企业信息以及历史结果的多源数据为先验信息，以供进行结
果验证。只能尽可能地提高污染源位置推算精度，无法准确定性每个污染源对受体点的贡
献比率。
发明内容
[0005] 本发明目的在于提供一种工业园区大气污染物溯源方法，能够在不同风向风速及源强排放情况下，实现对工业园区大气污染物的动态精准溯源。
[0006] 为达成上述目的，结合图1，本发明提出一种工业园区大气污染物溯源方法，所述溯源方法包括以下步骤：
[0007] S1，采集工业园区详细成分站及微型站的污染物监测数据，根据污染物监测数据实时动态进行超标判断；
[0008] S2，依据超标判断的结果，得到超标因子与对应的超标点位，将超标因子与所有园区企业的污染物特征因子指纹库比对，初步筛选疑似污染源，将筛选出的疑似污染源的排
放贡献因子EF1赋值为1，其余排放源的排放贡献因子EF1赋值为0；
[0009] S3，采集工业园区气象监测数据，根据气象监测数据与园区内企业分布情况，确定污染溯源的时间T和区间S；
[0010] S4，基于拉格朗日数值模型，在溯源追踪时段内，逐分钟动态计算疑似污染源对受体点的浓度贡献情况Mmn，Mmn为不同疑似污染源对受体点的浓度贡献系数；
[0011] S5，采集企业在线工况信息及企业排口在线监测数据，在超标点位建立归一化分布滞后统计模型对各疑似污染源的排放贡献因子EF2和EF3进行分析；
[0012] S6，结合疑似污染源对受体点的浓度贡献情况Mmn和各疑似污染源的排放贡献因子EF1、EF2、EF3，计算得到各疑似污染源对受体点的浓度贡献综合指数Con。
[0013] 进一步的实施例中，步骤S1中，所述根据污染物监测数据实时动态进行超标判断包括以下步骤：
[0014] S11，收集基础污染物监测数据，建立园区大气污染物基本浓度水平库；
[0015] S12，根据所有污染物的浓度实时变率和相关浓度标准限值，进行动态超标判断。
[0016] 进一步的实施例中，步骤S2中，所述初步筛选疑似污染源的过程包括以下步骤：
[0017] S21，对园区企业进行废气排放梳理，建立污染物特征因子指纹库；
[0018] S22，将超标因子与污染物特征因子指纹库进行比对，初步筛选疑似污染源；
[0019] S23，针对疑似污染源，定义其排放贡献因子EF1赋值为1，针对其余排放源，维持其排放贡献因子EF1赋值为0。
[0020] 进一步的实施例中，步骤S3中，所述确定污染溯源的时间及区间的过程包括以下步骤：
[0021] S31，收集园区基础气象数据；
[0022] S32，根据园区企业分布，选取疑似污染源到受体点的最大直线距离为S，根据实时变化的风速 按照如下公式进行积分求解，确定污染溯源的时间T及区间S：
[0023]
[0024] 进一步的实施例中，步骤S4中，所述逐分钟动态计算疑似污染源对受体点的浓度贡献情况Mmn的过程包括以下步骤：
[0025] S41，建立拉格朗日数值模型；
[0026] S42，在溯源追踪时段内，在受体点位置释放大量粒子，基于三维空间和时间对疑似污染源输入的变化进行响应函数的计算，得到疑似污染源与受体点之间的传输关系矩阵
以及随时间序列变化的网格点气团驻留时间；
[0027] S43，进行后向源受体关系和潜在源贡献分析。
[0028] 进一步的实施例中，步骤S42中，在计算所述疑似污染源与受体点之间的传输关系矩阵时，采用下述公式计算网格点气团驻留时间τi：
[0029]
[0030] 式中，Ts为时间分辨率，单位是秒，N为Ts时间范围内采样的数量，J为释放的粒子总数，fijn代表对于指定个点有贡献的粒子数量。
[0031] 进一步的实施例中，通过下述公式计算得到疑似污染源对受体点的浓度贡献情况Mmn：
[0032] Ym＝MmnXn
[0033] Xn＝[x1 … xn]
[0034] 式中，Ym代表模拟得到的浓度值矩阵；Xn为未知量，代表需要反演求解得到的排放源矩阵；Mmn矩阵代表不同网格点排放对受体点的浓度贡献系数矩阵；m代表受体点的个数，n
代表疑似污染源的个数。
[0035] 进一步的实施例中，步骤S5中，所述在超标点位建立归一化分布滞后统计模型对各疑似污染源的排放贡献因子EF2和EF3进行分析的过程包括以下步骤：
[0036] S51，采集园区内企业的实时在线工况数据以及企业排口在线监测数据，建立园区在线生产排放信息库；
[0037] S52,将生产排放信息与超标点位建立归一化分布滞后统计模型进行分析：
[0038] S521,采用反余切函数转换的非线性归一化方法，将步骤S51中采集到的所有数据进行归一化处理，将数据映射到0至1范围内，形成无量纲量；
[0039] S522,采用下述公式建立企业在线工况数据和企业排口在线监测数据与超标点位的归一化分布滞后统计模型：
[0040] At＝α+β0Bt+β1Bt‑1+β2Bt‑2+…+βsBt‑s+ut
[0041] 式中：
[0042] s为滞后长度，根据步骤S3中确定的溯源时间T进行动态调整；
[0043] β0为短期乘数，表示变量B变动一个单位时，对值A的平均影响大小；
[0044] βi为延迟乘数，表示过去各时期变量B变动一个单位对A值的平均影响大小；
[0045] 为长期乘数，表示变量B变动一个单位时，由于滞后效应而形成的对A值总的影响大小；α为常量，ut为误差，t为当前时刻。
[0046] S523,结合归一化分布滞后统计模型，计算得到各排放点位的排放贡献因子EF2和EF3。
[0047] 进一步的实施例中，步骤S6中，所述计算得到各疑似污染源对受体点的浓度贡献综合指数Con是指：
[0048] 根据下述公式计算得到各疑似污染源对受体点的浓度贡献综合指数Con：
[0049] Con＝EF1×EF2×EF3×Mmn。
[0050] 以上本发明的技术方案，与现有相比，其显著的有益效果在于：
[0051] (1)本发明通过分布式部署微型站，精确获取每个区域的监测数据，通过详细成分站分析出其中的超标因子和对应的受体点，以及对应的疑似污染源，减少后续污染源溯源
计算量，提高污染物溯源结果精确性。
[0052] (2)结合气象站和企业在线工况信息及企业排口在线监测数据，计算得到每个疑似污染源对受体点的贡献情况，并且从疑似污染源中进一步锁定真正的污染源，计算每个
污染源或园区企业对超标点位的贡献率，在不同风向风速及源强排放情况下，实现对工业
园区大气污染物的动态精准溯源。
[0053] (3)基于拉格朗日数值模型，在溯源追踪时段内，逐分钟动态计算疑似污染源对受体点的浓度贡献情况Mmn，计算结果精确度高，溯源效率高。
[0054] 应当理解，前述构思以及在下面更加详细地描述的额外构思的所有组合只要在这样的构思不相互矛盾的情况下都可以被视为本公开的发明主题的一部分。另外，所要求保
护的主题的所有组合都被视为本公开的发明主题的一部分。
[0055] 结合附图从下面的描述中可以更加全面地理解本发明教导的前述和其他方面、实施例和特征。本发明的其他附加方面例如示例性实施方式的特征和/或有益效果将在下面
的描述中显见，或通过根据本发明教导的具体实施方式的实践中得知。
附图说明
[0056] 附图不意在按比例绘制。在附图中，在各个图中示出的每个相同或近似相同的组成部分可以用相同的标号表示。为了清晰起见，在每个图中，并非每个组成部分均被标记。
现在，将通过例子并参考附图来描述本发明的各个方面的实施例，其中：
[0057] 图1是本发明的工业园区大气污染物溯源方法的流程图。
[0058] 图2是本发明的某个时刻各个网格点对受体点A的浓度贡献空间分布情况示意图。
[0059] 图3是本发明的某个时刻各个网格点对受体点B的浓度贡献空间分布情况示意图。
具体实施方式
[0060] 为了更了解本发明的技术内容，特举具体实施例并配合所附图式说明如下。
[0061] 结合图1，本发明提出一种工业园区大气污染物溯源方法，所述溯源方法包括以下步骤：
[0062] S1，采集工业园区详细成分站及微型站的污染物监测数据，根据污染物监测数据实时动态进行超标判断。
[0063] S2，依据超标判断的结果，得到超标因子与对应的超标点位，将超标因子与所有园区企业的污染物特征因子指纹库比对，初步筛选疑似污染源，将筛选出的疑似污染源的排
放贡献因子EF1赋值为1，其余排放源的排放贡献因子EF1赋值为0。
[0064] S3，采集工业园区气象监测数据，根据气象监测数据与园区内企业分布情况，确定污染溯源的时间T和区间S。
[0065] S4，基于拉格朗日数值模型，在溯源追踪时段内，逐分钟动态计算疑似污染源对受体点的浓度贡献情况Mmn，Mmn为不同疑似污染源对受体点的浓度贡献系数。
[0066] S5，采集企业在线工况信息及企业排口在线监测数据，在超标点位建立归一化分布滞后统计模型对各疑似污染源的排放贡献因子EF2和EF3进行分析。
[0067] S6，结合疑似污染源对受体点的浓度贡献情况Mmn和各疑似污染源的排放贡献因子EF1、EF2、EF3，计算得到各疑似污染源对受体点的浓度贡献综合指数Con。
[0068] 步骤一、根据污染物监测数据实时动态进行超标判断
[0069] 步骤S1中，所述超标判断包括：收集基础监测数据，建立园区大气污染物基本浓度水平库，根据污染物浓度的实时变率及相关浓度标准限值进行动态超标判断等。
[0070] 优选的，根据实际需求，在园区内分布部署多个微型站，以对各个重点关注区域或全部区域进行实时监控。
[0071] 应当理解，在本发明中，超标判断可以在微型站中进行，也可以由微型站将监测数据传输回详细成分站对监测数据做进一步的详细成分分析和超标判断，锁定其中的超标因
子，再根据监测数据所对应的微型站的安装位置确定超标点位(受体点)。
[0072] 步骤二、初步筛选疑似污染源
[0073] 首先对园区企业进行废气排放梳理，建立污染特征因子指纹库，其次将超标因子与污染物特征因子指纹库进行比对，初步筛选疑似污染源，针对疑似污染源，定义其排放贡
献因子EF1赋值为1，针对其余排放源，维持其排放贡献因子EF1赋值为0，以作区别。
[0074] 污染特征因子指纹库中录入有每个园区企业的名称、排放源(排放点)位置和每个排放点对应的污染物因子。当发现存在超标因子时，将超标因子与污染物特征因子指纹库
的污染物因子作对比，筛选出比对成功的排放源定义成疑似污染源，将其排放贡献因子EF1
赋值更新为1，其他排放源仍维持初始赋值为0。
[0075] 当存在不止一个超标因子时，可以针对每个超标因子逐一分析，筛选出若干组疑似污染源，再分别进行下述步骤的污染溯源分析过程，找出每个超标因子对应的污染物溯
源分析结果。
[0076] 步骤三、确定污染溯源的时间及区间
[0077] 具体的，首先收集园区基础气象数据，根据园区企业分布，选取疑似污染源到受体点的最大直线距离为S，根据实时变化的风速 按照如下公式进行积分求解，进而确定污
染溯源的时间(T)及区间(S)： 由此还可以删除疑似污染源中不符合时间和区
间的部分，进一步缩小疑似污染源的范围。
[0078] 步骤四、逐分钟动态计算疑似污染源对受体点的浓度贡献情况Mmn
[0079] 在本发明中，基于拉格朗日数值模型对疑似污染源对受体点的浓度贡献情况Mmn进行计算。所述拉格朗日数值模型是利用统计学方法跟踪大量粒子的运动轨迹，通过计算虚
拟粒子在空间中的运动，得到其在时空中的轨迹，从而开展其他计算。核心部分是如何描述
气团的运动。对于某个粒子，假设有它的初始位置P(t)，得到第一猜测位置P'(t+Δt)：
[0080] P′(t+Δt)＝P(t)+V(P,t)Δt。
[0081] 所以最终位置为：
[0082] P(t+Δt)＝P(t)+0.5[V(P,t)+V(P′,t+Δt)]Δt。
[0083] 精确到一阶，所用到的积分轨迹方程为：
[0084]
[0085] 式中t为时间，Δt为时间步长，V为风速。
[0086] 在溯源追踪时段内，首先在受体点位置释放大量粒子，然后基于三维空间和时间对排放源输入的变化进行响应函数的计算，得到排放源与受体之间的传输关系矩阵以及随
时间序列变化的网格点气团驻留时间，进而进行后向源受体关系和潜在源贡献分析。
[0087] 其中，在计算所述疑似污染源与受体点之间的传输关系矩阵时，采用下述公式计算网格点气团驻留时间τi：
[0088]
[0089] 式中，Ts为时间分辨率，单位是秒，N为Ts时间范围内采样的数量，J为释放的粒子总数，fijn代表对于指定个点有贡献的粒子数量。
[0090] 最终通过下述公式计算得到疑似污染源对受体点的浓度贡献情况Mmn：
[0091] Ym＝MmnXn
[0092] Xn＝[x1 … xn]
[0093] 式中，Ym代表模拟得到的浓度值矩阵；Xn为未知量，代表需要反演求解得到的排放源矩阵；Mmn矩阵代表不同网格点排放对受体点的浓度贡献系数矩阵；m代表受体点的个数，n
代表疑似污染源的个数。
[0094] 步骤五、在超标点位建立归一化分布滞后统计模型对各疑似污染源的排放贡献因子EF2和EF3进行分析
[0095] 该步骤包括以下步骤：采集园区内企业的实时在线工况数据以及企业排口在线监测数据，建立园区在线生产排放信息库；将生产排放信息与超标点位建立归一化分布滞后
统计模型进行分析；最后结合归一化分布滞后统计模型，计算得到各排放点位的排放贡献
因子EF2和EF3。其中，建立归一化分布滞后统计模型的过程包括：
[0096] S521,采用反余切函数转换的非线性归一化方法，将步骤S51中采集到的所有数据进行归一化处理，将数据映射到0至1范围内，形成无量纲量，便于数据处理。采用的归一化
公式如下：
[0097]
[0098] S522,采用下述公式建立企业在线工况数据和企业排口在线监测数据与超标点位的归一化分布滞后统计模型：
[0099] At＝α+β0Bt+β1Bt‑1+β2Bt‑2+…+βsBt‑s+ut
[0100] 式中：
[0101] s为滞后长度，根据步骤S3中确定的溯源时间T进行动态调整；
[0102] β0为短期乘数，表示变量B变动一个单位时，对值A的平均影响大小；
[0103] βi为延迟乘数，表示过去各时期变量B变动一个单位对A值的平均影响大小；
[0104] 为长期乘数，表示变量B变动一个单位时，由于滞后效应而形成的对A值总的影响大小；α为常量，ut为误差，t为当前时刻。
[0105] 步骤六、计算得到各疑似污染源对受体点的浓度贡献综合指数Con
[0106] 具体的，根据模型溯源及统计溯源方法结果，建立综合评价指数，以实现工业园区动态精准溯源。其中，综合评价指数的建立包括：一方面基于统计溯源方法，根据步骤二及
步骤五中求解得到的排放贡献因子EF1、EF2和EF3，另一方面在步骤4中通过模型溯源计算得
到的疑似点源逐分钟动态对受体点的浓度贡献情况Mmn，最终建立各排放源对受体点的浓度
贡献综合指数Con，实现工业园区动态精准溯源，计算公式如下：
[0107] Con＝EF1×EF2×EF3×Mmn。
[0108] 图2和图3分别是同一个工业园区中，某个时刻各个网格点对受体点A和受体点B的浓度贡献空间分布情况示意图。
[0109] 在本公开中参照附图来描述本发明的各方面，附图中示出了许多说明的实施例。本公开的实施例不必定义在包括本发明的所有方面。应当理解，上面介绍的多种构思和实
施例，以及下面更加详细地描述的那些构思和实施方式可以以很多方式中任意一种来实
施，这是因为本发明所公开的构思和实施例并不限于任何实施方式。另外，本发明公开的一
些方面可以单独使用，或者与本发明公开的其他方面的任何适当组合来使用。
[0110] 虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因
此，本发明的保护范围当视权利要求书所界定者为准。
