CN111855915A - 一种大气污染物的溯源方法 - Google Patents

Abstract

本发明公开了一种大气污染物溯源方法，包括以下步骤：1)利用大气化学模式得到源汇敏感性矩阵；2)通过步骤1)得到的源汇敏感性矩阵计算大气化学模式每个网格点上释放污染物产生的模拟浓度与观测浓度或残差浓度的相关系数，并将最大相关系数所在位置确定为污染源位置；3)计算步骤2)中确定的污染源的排放强度；4)利用步骤2)得到的相关系数以及步骤3)得到的排放强度判断是否所有的污染源均已被识别，若是，则输出污染源的个数、每个污染源的位置以及排放强度；若否，则重复步骤2)‑3)，直到所有的污染源均被识别。本发明可应对突发性的大气污染事故，减少事故对生态环境及人民健康造成的危害。

Description

一种大气污染物的溯源方法

技术领域

本发明涉及污染物溯源技术领域，涉及一种大气污染物的溯源方法。

背景技术

目前的污染物溯源方法主要结合大气化学模式和观测浓度信息来反演污染源的特征，国内外的研究学者已提出较多的源项估计方法，例如较为简单的后向轨迹法以及较为复杂的概率方法。

后向轨迹法直接将大气化学模式的风场反向，计算所有观测数据的后向拉格朗日轨迹，并将所有轨迹的交点定为污染源所在位置。其优势在于简单易行，计算代价较小，但仅能做定性分析，不能估计污染源的强度，并且由于在计算后向轨迹时仅考虑了平流的作用，没有考虑扩散过程，因此估计结果在很大程度上依赖于风场的精确程度。

概率方法基于贝叶斯推论将可用的浓度观测数据与先验信息相结合，并考虑输入数据的误差(如观测误差、模式误差及先验误差)，最终通过采样估计得到后验参数的不确定性及置信区间。虽然该方法能够提供估计参数的不确定性信息，但是由于受到较为繁杂的采样过程的限制，这类方法通常较为耗时，且需要已知污染源的先验信息，这在大气应急响应中是难以实现的。

上述方法主要集中于单个污染源的估计，而对多污染源的溯源方法较少。但是在真实的状况中，如工厂的烟囱排放、毒性气体的突然泄露，一般都属于多点源排放事件。当存在多个点源时，观测到的浓度是多个污染源产生的浓度场的累加，难以区分每个污染源的贡献，并且观测点附近的弱排放源可能会与距观测点较远的强排放源对观测站点产生相同的影响，这使得多点源估计要比单点源估计更具挑战性。

现有技术中仅有少数几种方法可以实现对多个污染源的个数、位置及强度的精确估计，其中需要事先给定某些先验信息，或者计算代价会随着污染源个数的增加而显著增加，使其很难直接应用于大气污染的应急响应中。

因此，需要建立快速且准确的多点污染源估计方法，以应对突发性污染事故的发生，为应急决策的制定提供参考依据并保障人民的生命健康。

发明内容

针对现有技术存在的问题，本发明的目的在于提供一种较为快速且准确的大气污染物溯源方法，以应对突发性的大气污染事故，减少事故对生态环境及人民健康造成的危害。

为实现上述目的，本发明提供如下技术方案：

一种大气污染物溯源方法，包括以下步骤：

1)利用大气化学模式得到源汇敏感性矩阵；

2)通过步骤1)得到的源汇敏感性矩阵计算大气化学模式每个网格点上释放污染物产生的模拟浓度与观测浓度或残差浓度的相关系数，并将最大相关系数所在位置确定为污染源位置；

3)计算步骤2)中确定的污染源的排放强度；

4)利用步骤2)得到的相关系数以及步骤3)得到的排放强度判断是否所有的污染源均已被识别，

若是，则输出污染源的个数、每个污染源的位置以及排放强度；

若否，则重复步骤2)-3)，直到所有的污染源均被识别。

进一步，步骤1)中，以每个观测数据所在位置为起点，输入气象预报的风场信息，后向运行大气化学模式，可得到一组源汇敏感性矩阵H，矩阵中的每个元素H_ij代表第i个观测浓度数据对第j个模式格点上的污染排放通量的敏感性。

进一步，步骤2)中，在第一次计算时，设置潜在污染源的个数n＝1。

进一步，步骤2)中，通过步骤1)中得到的源汇敏感性矩阵H与污染源排放矩阵相乘，得到在大气化学模式的每个网格点上释放污染物产生的模拟浓度。

进一步，步骤2)中，将大气化学模式的所有网格点作为潜在污染源。

进一步，步骤2)中，每个潜在污染源处释放污染物所产生的模拟浓度与观测浓度之间的相关系数为：

其中，j为大气化学模式网格点的编号，N为总的网格点数，i为观测浓度的编号，m为观测数据的总个数，

为第i个观测浓度数据，H_ij为步骤1)中得到的源汇敏感性矩阵中的元素，q^s为一个包含污染源排放强度的向量，H_ij与q^s的乘积为污染物在第j个网格点上释放所产生的与第i个观测数据所对应的模拟浓度。

进一步，步骤2)中，将得到的N个相关系数的最大值记为R_max，n，并将其所处的网格点确定为第n个污染源位置。

进一步，步骤2)中，还包括存储每次计算中的最大相关系数R_max，n及其所在位置(x_n，y_n)，并将所有已经计算过的最大相关系数的最大值记录为R_max，R_max＝MAX(R_max，p)，p＝1，2，...，n。

进一步，步骤3)中，通过最小化代价函数计算已识别的所有n个污染源的排放强度，代价函数定义为：

其中，H_e为m×n维矩阵，代表n个污染源的排放与m个观测浓度之间的源汇关系，q为待求的n个污染源的排放强度向量，C°为m×1维的观测浓度向量。

进一步，步骤3)与步骤4)之间还包括计算观测浓度与模拟浓度之间的残差浓度，残差浓度为：

其中，q_e为步骤3)中得到的n个污染源的排放强度。

进一步，步骤4)中，判断所有的污染源均被识别的判断条件为：

如果在第n次计算得到的N个相关系数中的最大值R_max，n小于之前所有计算过程中最大相关系数R_max的二分之一，或者已经计算得到的污染源排放强度中最大值

与最小值

相差一个量级。

进一步，步骤4)中，所有的污染源未被识别进行下一次计算时，设置污染源的个数n＝n+1，且在计算相关系数时利用残差浓度代替观测浓度。

本方法可以较为快速且准确地定位污染源，并估计污染源的个数及相对应的污染源的排放强度，且在估计过程中仅需要大气化学模式及观测数据，不需要任何的初猜信息，非常适用于突发性污染事故的应急响应，能够迅速定位出源头，帮助决策者制定合理的疏散措施，以保障人民的生命财产安全。

附图说明

图1为本发明大气污染物溯源方法的***流程图。

具体实施方式

下面结合具体实施方式对本发明的技术方案作进一步详细地说明。

本发明中的大气污染物溯源方法，主要的步骤包括：

1)利用大气化学模式得到源汇敏感性矩阵；

2)计算大气化学模式每个网格点上释放污染物产生的模拟浓度与观测浓度或残差浓度的相关系数，并将最大相关系数所在位置确定为估计的污染源的位置；

3)以最小化代价函数J计算确定的污染源的排放强度；

4)利用计算得到的相关系数以及污染源的排放强度判断是否所有主要污染源均已被识别，若是，则结束迭代计算，输出污染源的个数、污染源的位置以及各个污染源的排放强度；若否，则重复步骤2)-3)，继续下一次迭代计算，利用残差浓度估计未被识别的污染源，污染源的个数相应的增加，直到所有的污染源均被识别，输出，结束迭代计算。

本发明中上述步骤的***流程参见图1。

以下对上述步骤从第一次计算开始，进行展开说明：

利用大气化学模式得到源汇敏感性矩阵：

以每个观测数据所在位置为起点，输入气象预报的风场信息，后向运行大气化学模式，可得到一组源汇敏感性矩阵H，矩阵中的每个元素H_ij代表第i个观测浓度数据对第j个模式格点上的污染排放通量的敏感性。

在后续计算过程中，直接将该矩阵与污染源排放矩阵相乘，便可得到相应排放所产生的模拟浓度，避免了较为耗时的正向大气化学模式的运行，可大大提高计算效率。

计算大气化学模式每个网格点上释放污染物产生的模拟浓度与观测浓度或残差浓度的相关系数：

首先，设置初始污染源的个数n＝1，由于在应急响应中一般难以提供潜在污染源的位置信息，这里直接将大气化学模式的所有网格点作为潜在污染源，计算在每个网格点上释放污染物产生的模拟浓度与观测浓度之间的相关系数为：

其中，j为大气化学模式网格点的编号，N为总的网格点数，i为观测浓度的编号，m为观测数据的总个数，

为第i个观测浓度数据，H_ij为大气化学模式后向运行得到的源汇敏感性矩阵中的元素，q^s为一个包含污染源排放强度的向量，H_ij与q^s的乘积即为污染物在第j个网格点上释放所产生的与第i个观测数据所对应的模拟浓度。

需要指出，由于这里考虑污染源是定常释放的，q^s为一常数，其取值不会影响两组变量之间的相关系数，因此在此方法的实施中可直接将q^s中的元素设为单位排放强度以简化计算。

寻找所有网格点上释放污染物产生的模拟浓度与观测浓度之间最大的相关系数，即寻找得到的N个相关系数的最大值R_maxn，并将其所处的网格点定为第1个污染源的估计位置，存储第一次迭代计算的最大相关系数R_max，1及其所在位置(x₁，y₁)。

计算确定的污染源的排放强度：

通过最小化代价函数计算第1个污染源的排放强度(源强)，在该计算步骤中，将代价函数定义为模拟浓度与观测浓度之间残差的平方和，即

其中，H_e为m×1维矩阵，代表第1个污染源的排放与m个观测浓度之间的源汇关系，q为待估计的第1个污染源的排放强度向量，C°为m×1维的观测浓度向量。

在最小化过程中，采用L-BFGS-B算法(具有约束的、在有限内存中进行的BFGS算法，BFGS为Broyden，Fletcher，Goldfarb，Shanno四位数学家提出的一种局部最小化算法)进行正定约束，可以避免估计的排放强度q中出现不合理的负值。

至此，已估计得到第1个污染源的位置和排放强度信息，再利用得到的源汇敏感性矩阵快速计算第1个污染源产生的模拟浓度，以及观测浓度与模拟浓度之间的残差浓度，即

其中q_e1为计算污染源排放强度时估计得到的第1个污染源的排放强度，

为残差浓度。

再进行第二次迭代计算，确定是否有第2个污染源信号，但是在计算相关系数以及N个相关系数的最大值时，利用残差浓度

代替观测浓度C°。

存储第二次迭代计算的最大相关系数R_max，2及其所在的位置(x₂，y₂)。

计算确定的第2个污染源的排放强度，在得到第2个污染源的排放强度后，即进行是否识别到所有污染源的判定。如果第二次迭代计算中得到的最大相关系数R_max，2小于第一次迭代计算中得到的最大相关系数R_max，1的二分之一，或者估计的污染源排放强度中q_e2与q_e1相差一个量级，则结束操作，输出污染源的个数为1，该污染源的位置及该污染源的排放强度等信息。

若未满足上述终止条件，说明第二个污染源也是重要污染，污染源数加1，同时根据新识别的污染源的位置，通过最小二乘优化算法重新估计所有被识别到的污染源的排放强度。

至此，已估计得到2个污染源的位置和排放强度信息，再利用得到的源汇敏感性矩阵快速计算这2个污染源产生的模拟浓度，以及观测浓度与模拟浓度之间的残差浓度，即

其中q_e为计算污染源排放强度时估计得到的2个污染源的排放强度，

为残差浓度。

利用得到的残差浓度

再进行第n次迭代计算，确定是否有第n个污染源信号，即进行第n个污染源的识别。

至此，已估计得到n个污染源的位置和排放强度，再利用得到的源汇敏感性矩阵快速计算这n个污染源产生的模拟浓度，以及观测浓度与模拟浓度之间的残差浓度，即

其中q_e为计算污染源排放强度时估计得到的n个污染源的排放强度，

为残差浓度。

如果此次迭代中的计算得到的最大相关系数R_max，n小于之前所有迭代过程中最大相关系数R_max的二分之一，或者估计的污染源排放强度中最大值

与最小值

相差一个量级，则结束操作，输出污染源的个数n，每个污染源的位置及污染源的排放强度等信息。

若未满足上述终止条件，说明仍有重要污染源未被成功识别，则令n＝n+1，重复计算，继续进行下一次迭代。在下一次迭代的位置估计中，利用残差

代替观测浓度C°以检验残差中是否包含其他的污染源信号，如有则确定新识别的污染源的位置。另外，为了更好地估计污染源特征，当有新的污染源被识别时，通过最小二乘优化算法重新估计所有被识别到的污染源的排放强度。

通过本发明中大气污染物的溯源方法，解决了目前研究较少的多点源估计问题，可同时估计污染源的个数、位置及源强信息，且由于使用源汇关系，其计算代价不会随着点源个数的增加而显著增加。

该方法采用了后向大气化学模式计算源汇关系，可大大提高计算效率。该方法不需要任何的初猜信息，非常适用于事故的应急响应。

需要重点指出，本发明的目的是提供一种较为快速且准确的多点源估计方法以应对突发性的大气污染事故，减少事故对生态环境及人民健康造成的危害。该方法主要采用迭代过程逐步地估计每个污染源的位置及强度信息，在每个污染源的估计中，使用遍历法并利用大气化学模式计算所有潜在污染源处释放污染物产生的模拟浓度与观测浓度的之间的相关系数，将相关系数最大值所在位置定为污染源的位置，然后利用最小二乘的优化方法估计该位置上的源强。在每次迭代中，需要判断此次迭代中的最大相关系数及估计的污染源强度是否足够大，若小于预先设定的标准，则说明此时所有的主要污染源均已被识别，可结束循环迭代过程，输出估计得到的污染源个数、位置及强度信息。

需要注意的是，本发明中的设置的迭代结束准则可以根据实际状况进行调整，以方便决策者制定更为合理的应急措施。

上面对本发明的较佳实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域的普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (10)

1.一种大气污染物溯源方法，其特征在于，包括以下步骤：

1)利用大气化学模式得到源汇敏感性矩阵；

2)通过步骤1)得到的源汇敏感性矩阵计算大气化学模式每个网格点上释放污染物产生的模拟浓度与观测浓度或残差浓度的相关系数，并将最大相关系数所在位置确定为污染源位置；

3)计算步骤2)中确定的污染源的排放强度；

4)利用步骤2)得到的相关系数以及步骤3)得到的排放强度判断是否所有的污染源均已被识别，

若是，则输出污染源的个数、每个污染源的位置以及排放强度；

若否，则重复步骤2)-3)，直到所有的污染源均被识别。

2.根据权利要求1所述的大气污染物溯源方法，其特征在于，步骤1)中，以每个观测数据所在位置为起点，输入气象预报的风场信息，后向运行大气化学模式，得到源汇敏感性矩阵H，矩阵中的每个元素H_ij代表第i个观测浓度数据对第j个模式格点上的污染排放通量的敏感性。

3.根据权利要求2所述的大气污染物溯源方法，其特征在于，步骤2)中，通过步骤1)中得到的源汇敏感性矩阵H与污染源排放矩阵相乘，得到在大气化学模式的每个网格点上释放污染物所产生的模拟浓度。

4.根据权利要求2所述的大气污染物溯源方法，其特征在于，步骤2)中，每个网格点上释放污染物所产生的模拟浓度与观测浓度之间的相关系数为：

其中，j为大气化学模式网格点的编号，N为总的网格点数，i为观测浓度的编号，m为观测数据的总个数，

为第i个观测浓度数据，H_ij为步骤1)中得到的源汇敏感性矩阵中的元素，q^s为一个包含污染源排放强度的向量，H_ij与q^s的乘积为污染物在第j个网格点上释放所产生的与第i个观测数据所对应的模拟浓度。

5.根据权利要求4所述的大气污染物溯源方法，其特征在于，步骤2)中，将得到的N个相关系数的最大值记为R_max，n，并将其所处的网格点确定为第n个污染源位置。

6.根据权利要求5所述的大气污染物溯源方法，其特征在于，步骤2)中，还包括存储每次计算中的最大相关系数R_max，n及其所在位置(x_n，y_n)，并将所有已经计算过的最大相关系数的最大值记录为R_max，R_max＝MAX(R_max，p)，p＝1，2，...，n。

7.根据权利要求6所述的大气污染物溯源方法，其特征在于，步骤3)中，通过最小化代价函数计算已识别的所有n个污染源的排放强度，代价函数定义为：

其中，H_e为m×n维矩阵，代表n个污染源的排放与m个观测浓度之间的源汇关系，q为待求的n个污染源的排放强度向量，C^o为m×1维的观测浓度向量。

8.根据权利要求7所述的大气污染物溯源方法，其特征在于，步骤3)与步骤4)之间还包括计算观测浓度与模拟浓度之间的残差浓度，残差浓度为：

其中，q_e为步骤3)中得到的n个污染源的排放强度。

9.根据权利要求8所述的大气污染物溯源方法，其特征在于，步骤4)中，判断所有的污染源均被识别的判断条件为：

如果在第n次计算得到的N个相关系数中的最大值R_max，n小于之前所有计算过程中最大相关系数R_max的二分之一，或者已经计算得到的污染源排放强度中最大值

与最小值

相差一个量级。

10.根据权利要求9所述的大气污染物溯源方法，其特征在于，步骤4)中，所有的污染源未被识别进行下一次计算时，设置污染源的个数n＝n+1，且在计算相关系数时利用残差浓度代替观测浓度。

Images