CN104702685B

CN104702685B - 基于后向轨迹的污染源追踪方法及其系统

Info

Publication number: CN104702685B
Application number: CN201510106698.XA
Authority: CN
Inventors: 徐伟嘉; 刘靖伟; 黄慎; 李红霞; 梁健俊
Original assignee: Guangzhou Rising Sun Sincerity Breath Science And Technology Ltd; National Sun Yat Sen University
Current assignee: Guangzhou Rising Sun Sincerity Breath Science And Technology Ltd; National Sun Yat Sen University
Priority date: 2015-03-11
Filing date: 2015-03-11
Publication date: 2018-02-02
Anticipated expiration: 2035-03-11
Also published as: CN104702685A

Abstract

本发明公开了一种基于后向轨迹的污染源追踪方法，包括S1：根据污染物的后向轨迹计算后向轨迹经过的每一网格内污染物的输送概率、距离权重及时间权重；S2：根据源清单计算网格内污染物的排放强度；S3：根据输送概率、距离权重、时间权重及排放强度计算网格内污染物的输送强度；S4：根据各网格的输送强度计算后向轨迹上污染物的总输送强度；S5：取多个后向轨迹重复S1至S4以得到多个总输送强度；S6：将总输送强度最大的前N个后向轨迹经过的网格区域作为预测区域，根据预测区域及源清单确定污染源的地理位置。与现有技术相比，本发明实现了污染源的定位追踪，为解决空气污染问题提供了重要依据。本发明同时公开了一种基于后向轨迹的污染源追踪系统。

Description

基于后向轨迹的污染源追踪方法及其系统

技术领域

本发明涉及环境监测技术领域，更具体地涉及一种基于后向轨迹的污染源追踪方法及其系统。

背景技术

随着大气污染越来越严重，对污染源的分析也日趋重要。目前，计算污染物后向轨迹的模型有HYSPLIT，HYSPLIT-4模型是由美国国家海洋和大气管理局(NOAA)的空气资源实验室和澳大利亚气象局在过去20年间联合研发的一种用于计算和分析大气污染物输送、扩散轨迹的专业模型，该模型具有处理多种气象要素输入场、多种物理过程和不同类型污染物排放源功能的较为完整的输送、扩散和沉降模式,已经被广泛地应用于多种污染物在各个地区的传输和扩散的研究中。

然而，上述HYSPLIT模型采用的天气数据是111KM*111KM的数据，这个数据不能很好的模拟局部气流的运行状态，很有可能由于局部气流的影响而改变污染物的输送轨迹，因此通过HYSPLIT模型计算污染物后向轨迹时的精度不够；同时，HYSPLIT模型不能解决污染源追踪的问题，即不能准确的定位污染源的地理位置，从而不能为解决日益严重的空气污染提供依据。

因此，有必要提供一种基于后向轨迹的污染源追踪方法及其系统来来解决污染源追踪的问题。

发明内容

本发明的目的在于提供一种基于后向轨迹的污染源追踪方法，以通过后向轨迹结合污染源清单来确定污染源的地理位置，进而为解决空气污染提供依据。

为实现上述目的，本发明提供了一种基于后向轨迹的污染源追踪方法，包括：

S1：根据污染物的后向轨迹依次计算所述后向轨迹经过的每一网格内所述污染物的输送概率、距离权重以及时间权重；

S2：根据所述污染物的源清单计算所述网格内所述污染物的排放强度；

S3：根据所述输送概率、距离权重、时间权重以及排放强度计算所述网格内所述污染物的输送强度，计算公式为：

T_l(i,j)＝R_l(i,j)*E(i,j)*W_d(i,j)*W_tl(i,j)

其中，T_l(i,j)为网格(i,j)内所述污染物的输送强度，R_l(i,j)为网格(i,j)内所述污染物的输送概率，E(i,j)为网格(i,j)内所述污染物的排放强度，W_d(i,j)为网格(i,j)内所述污染物的距离权重，W_tl(i,j)为网格(i,j)内所述污染物的时间权重；

S4：根据各所述网格内的所述输送强度计算所述后向轨迹上所述污染物的总输送强度；

S5：取所述污染物在不同时刻的所述后向轨迹重复步骤S1至S4以得到多条所述后向轨迹对应的多个所述总输送强度；以及

S6：取所述总输送强度最大的前N个所述后向轨迹经过的网格区域作为预测区域，根据所述预测区域以及所述源清单确定污染源的地理位置，其中N为正整数。

与现有技术相比，本发明基于后向轨迹的污染源追踪方法通过后向轨迹计算污染物在每一网格内的输送概率、距离权重以及时间权重，并通过源清单计算相应网格内污染物的排放强度，最后计算得到每个网格内污染物的输送强度，进而得到后向轨迹上污染物的总输送强度，然后通过计算多条后向轨迹的总输送强度以及将总输送强度最大的N个后向轨迹经过的网格区域作为预测区域，再将预测区域结合源清单确定污染源的地理位置，实现了污染源的定位追踪，且结果较为精确，为解决空气污染问题提供了重要依据。

较佳地，步骤S5具体为：

取1个月内的720条所述后向轨迹重复步骤S1至S4以得到所述后向轨迹对应的720个所述总输送强度，其中1个月内每隔1小时为所述污染物绘制一条所述后向轨迹。

较佳地，根据污染物的后向轨迹计算所述网格内所述污染物的输送概率具体为：

R_(i,j)＝τ_l(i,j)/n

其中τ_l(i,j)为所述后向轨迹在所述网格(i,j)内的停留时间，n为所述后向轨迹的运行总时间。

较佳地，所述停留时间以所述后向轨迹在所述网格(i,j)内的轨迹点个数表示，所述运行总时间以所述后向轨迹上的所有轨迹点个数表示。

较佳地，根据污染物的后向轨迹计算所述网格内所述污染物的距离权重具体为：

W_d(i,j)＝1/(d(i,j)/5+1)

其中d(i,j)为所述网格(i,j)与所述后向轨迹上所述污染物的目的地之间的距离。

较佳地，根据污染物的后向轨迹计算所述网格内所述污染物的时间权重具体为：

W_tl(i,j)＝1/((t_l(i,j))/18+1)

其中t_l(i,j)为所述污染物由所述网格(i,j)移动至所述后向轨迹上所述污染物的目的地所需的时间。

较佳地，根据所述污染物的源清单计算所述网格内所述污染物的排放强度具体为:

从所述源清单中读取所述网格以及3个相邻网格内所述污染物在不同高度层处的排放强度；

依次计算每一所述网格内所述污染物的总排放强度，所述总排放强度为各个高度层相应的所述排放强度的总和；

将所述网格以及相邻网格的所述总排放强度取平均值以作为所述网格内所述污染物的排放强度。

较佳地，所述后向轨迹的计算方法包括：

(1)从MM5模式生成的NetCDF文件中获取所述污染物在初始点的初始速度；

(2)根据所述初始点的位置坐标和所述初始速度计算预测点的位置坐标，计算公式为：P'(t+Δt)＝P(t)+V(P，t)*Δt；

(3)根据所述初始点和所述预测点的位置坐标对所述预测点进行修正，修正公式为：P(t+Δt)＝P(t)+0.5[V(P，t)+V(P'，t+Δt)]*Δt；

(4)将所述修正点作为所述初始点重复步骤(1)至步骤(3)以获取所述所述污染物的多个位置坐标进而得到所述所述污染物的所述后向轨迹；

其中，P'(t+Δt)为所述预测点的位置坐标，P(t)为所述初始点的位置坐标，V(P，t)为所述初始速度，Δt为所述预测点与所述初始点的时间差，P(t+Δt)为对所述预测点进行修正后得到的修正点的位置坐标，V(P'，t+Δt)为所述所述污染物在所述预测点的运行速度。

相应的，本发明还提供了一种基于后向轨迹的污染源追踪系统，包括：

第一参数计算模块，用于根据污染物的后向轨迹依次计算所述后向轨迹经过的每一网格内所述污染物的输送概率、距离权重以及时间权重；

第二参数计算模块，用于根据所述污染物的源清单计算所述网格内所述污染物的排放强度；

输送强度计算模块，用于根据所述输送概率、距离权重、时间权重以及排放强度计算所述网格内所述污染物的输送强度，计算公式为：

T_l(i,j)＝R_l(i,j)*E(i,j)*W_d(i,j)*W_tl(i,j)

总输送强度计算模块，用于根据各所述网格内的所述输送强度计算所述后向轨迹上所述污染物的总输送强度；

预测区域获取模块，用于将所述污染物在不同时刻的多个所述后向轨迹对应的多个所述总输送强度最大的前N个所述后向轨迹经过的网格区域作为预测区域，其中N为正整数；以及

地理位置确定模块，用于根据所述预测区域以及所述源清单确定污染源的地理位置。

通过以下的描述并结合附图，本发明将变得更加清晰，这些附图用于解释本发明的实施例。

附图说明

图1为本发明基于后向轨迹的污染源追踪方法一实施例的流程图。

图2为图1中网格的示意图。

图3为根据后向轨迹计算污染物的输送概率、时间权重时的参考示意图。

图4为根据后向轨迹计算污染物的距离权重时一实施例的参考示意图。

图5为根据后向轨迹计算污染物的距离权重时另一实施例的参考示意图。

图6为图1中步骤S102的子流程图。

图7为求解污染物的后向轨迹的方法一实施例的流程图。

图8为图7中步骤S201一优选实施例的子流程图。

图9为图8中计算初始速度的示意图。

图10为图7中步骤S201另一优选实施例的子流程图。

图11为图10中计算初始速度的示意图。

图12为图7中步骤S203计算修正点坐标的示意图。

图13为图7中步骤S204计算后向轨迹的示意图。

图14为本发明基于后向轨迹的污染源追踪系统一实施例的结构框图。

图15为图14中第一参数计算模块的结构框图。

图16为图14中第二参数计算模块的结构框图。

具体实施方式

现在参考附图描述本发明的实施例，附图中类似的元件标号代表类似的元件。

请参考图1，本发明基于后向轨迹的污染源追踪方法包括以下步骤：

步骤S101，根据污染物的后向轨迹l依次计算污染物轨迹l经过的每一网格内污染物的输送概率R_(i,j)、距离权重W_d(i,j)以及时间权重W_tl(i,j)；

需要说明的是，此处以及下面提及的污染物对应于同一污染物，即本发明是以某一特定污染物为例进行说明的，计算过程具体如下：

(1)根据污染物的后向轨迹l计算网格内污染物的输送概率R_(i,j)的公式为：

R_(i,j)＝τ_l(i,j)/n

其中τ_l(i,j)为污染物的后向轨迹l在网格(i,j)内的停留时间，n为后向轨迹l的运行总时间。由于后向轨迹l在网格内的停留时间获取较为困难，因此，可以通过网格内轨迹点的个数来表示时间，优选的，停留时间以后向轨迹l在网格(i,j)内的轨迹点个数表示，运行总时间以后向轨迹l上的所有轨迹点个数表示。举例说明如下：首先请参考图2，其中网格(i,j)表示点(i,j)的右上角区域(阴影部分)，在MM5模式生成的NetCDF文件中，网格(i,j)的大小为3Km*3Km，再请参考图3，图3中黑色点的连线表示污染物的后向轨迹l，此时：

τ_l(0,0)＝3，即在网格(0，0)内有3个轨迹点，

τ_l(0,1)＝0，即在网格(0，1)内有0个轨迹点，

τ_l(1,2)＝0，即在网格(1，2)内有0个轨迹点……

且n＝12，即后向轨迹l上共12个轨迹点，

因此，R_(0,0)＝τ_l(0,0)/n＝3/12，即当i＝0，j＝0时，网格(i,j)内污染物的输送概率为3/12，同理可知其他网格内污染物的输送概率。

(2)根据污染物的后向轨迹l计算网格内污染物的距离权重W_d(i,j)的公式为：

W_d(i,j)＝1/(d(i,j)/5+1)

其中d(i,j)为网格(i,j)与后向轨迹l上污染物的目的地之间的距离。下面参考图3至图5进行详细说明：d(1,1)为图4中粗线部分所示，即后向轨迹l在网格(1,1)内的第一个轨迹点至污染物的目的地(后向轨迹l的初始点)之间的距离；同理，d(2,1)为图5中粗线部分所示，即后向轨迹l在网格(2,1)内的第一个轨迹点至污染物的目的地之间的距离。

(3)根据污染物的后向轨迹l计算网格内污染物的时间权重W_tl(i,j)的公式为：

W_tl(i,j)＝1/((t_l(i,j))/18+1)

其中t_l(i,j)为污染物由网格(i,j)移动至后向轨迹l上污染物的目的地所需的时间，其单位为h(小时)。下面参考图3进行详细说明：t_l(i,j)是通过计算后向轨迹l上轨迹点的个数、然后乘以Δt计算得到的，如t_l(1,1)＝9*Δt，即污染物由网格(1,1)移动至后向轨迹l上污染物的目的地所需的时间为9*Δt。其中Δt为后向轨迹上相邻轨迹点之间的时间差，而网格内轨迹点的个数取决于Δt的大小，Δt越小，网格内的轨迹点个数越多，Δt越大，网格内轨迹点的个数越少。

步骤S102，根据污染物的源清单计算网格内污染物的排放强度E(i,j)；

下面参考图6以网格(i,j)为例说明排放强度的计算过程，包括以下步骤：

步骤S1021，从NetCDF格式的源清单中读取网格(i,j)以及3个相邻网格内不同高度层处污染物的排放强度(即层排放强度)，其中网格(i,j)的相邻网格为网格(i,j+1)、网格(i+1,j)、网格(i+1,j+1)，其中将大气分为14个高度层，每一网格点对应有14个高度层，每个高度层都有一个排放强度(存储于源清单中)，通过读取源清单就能获取排放强度；

具体的，本实施例中污染物的源清单采用世界上较为先进的污染源排放模型SMOKE(Sparse Matrix Operator Kernel Emission System)，该模型将通过统计、监测手段获得的各种污染源排放数据制作成网格化的污染源排放数据，具体的，SMOKE模型将污染源分为点源、面源、移动源、生物源四种，每种污染源都有相应的规范，满足这些规范的污染源排放数据才能作为SMOKE模型的输入数据，而SMOKE的输出文件大多是NetCDF格式的，源清单中存储有污染物在不同高度层的排放强度。

步骤S1022，依次计算每一网格内的总排放强度，总排放强度为网格内各个高度层排放强度的总和；如网格(i,j)内的总排放强度为14个高度层上网格(i,j)内的排放强度的总和，同理可知网格(i,j+1)、网格(i+1,j)、网格(i+1,j+1)处的总排放强度；

步骤S1023，将网格(i,j)以及相邻网格内的总排放强度取平均值以作为网格(i,j)内污染物的排放强度E(i,j)，即4个总排放强度的平均值即为网格(i,j)内污染物的排放强度E(i,j)。

步骤S103，根据输送概率R_l(i,j)、距离权重W_d(i,j)、时间权重W_tl(i,j)以及排放强度E(i,j)计算网格(i,j)内污染物的输送强度T_l(i,j)，计算公式为：

T_l(i,j)＝R_l(i,j)*E(i,j)*W_d(i,j)*W_tl(i,j)

由此可以得到后向轨迹l上点(i.j)对应的网格(i.j)内污染物的输送强度T_l(i,j)；

步骤S104，根据各网格内的输送强度计算后向轨迹l上污染物的总输送强度；具体的，总输送强度等于各个网格内输送强度的总和；

步骤S105，取污染物在不同时刻的后向轨迹重复步骤S101至步骤S104以得到多条后向轨迹对应的多个总输送强度；优选的，可以选取1个月或2个月内污染物的后向轨迹进行计算，本实施例中取1个月内污染物的后向轨迹执行步骤S101至步骤S104，其中没隔1小时为污染物绘制一条后向轨迹，从而1个月内共有720条污染物的后向轨迹，进而可以得到720条后向轨迹对应的720个总输送强度；

步骤S106，取总输送强度最大的前N个后向轨迹经过的网格区域作为预测区域，根据预测区域以及源清单确定污染源的地理位置；其中N为正整数，如取N＝10，即取总输送强度最大的10个后向轨迹所经过的网格区域作为预测区域，然后结合源清单查找预测区域内污染物排放强度最大且对应有污染源(如化工厂、排污厂等)的地理位置，即为污染源的地理位置，如预测区域内存在的某化工厂对应的位置等。

需要说明的是，本实施例中对多条后向轨迹进行计算来确定预测区域，相较于依靠单条后向轨迹来确定污染源的预测区域，能够有效解决一条后向轨迹线由于测量误差、内插误差、计算误差等因素而导致测量结果存在较大误差的问题。

下面结合图7说明求解污染物后向轨迹的方法，其包括以下步骤：

S201：从MM5模式生成的NetCDF文件获取污染物在初始点P的初始速度V(P，t)；

其中获取初始速度V(P，t)的方式有多种，由于NetCDF文件中的数据网格大小为3KM*3KM，因此在对初始点P进行速度合成时，可以简单地取初始点P周围的四个网格点处风速的速度平均值作为初始速度，该结果非常接近污染物的实际运行速度，同时，将风速的矢量平均值作为初始速度，运算时间极小，使得所有程序在运行时1秒内运算结束，大大缩短了运算时间。如图8所示，步骤S201包括：

S2011：从NetCDF文件中读取初始点P周围的四个网格点处的风速，如图9所示，初始点P周围的四个网格点处的风速分别记为V(P1)、V(P2)、V(P3)、V(P4)；

S2012：计算风速V(P1)、V(P2)、V(P3)、V(P4)的矢量平均值以作为初始速度V(P，t)，其中风速包括U、V、W方向上的三个分量，U代表东、西方向风速，V代表南、北方向风速，W代表垂直地面风速，从而计算得到的初始速度V(P，t)也包括三个方向上的分量，计算得到的初始速度V(P，t)如图9所示。

在另一优选实施例中，如图10所示，为了使获取的初始速度V(P，t)准确度更高，还可以通过反距离加权法计算初始速度V(P，t)。此时，步骤S201具体包括：

S2011’：从NetCDF文件中读取初始点P周围四个网格点处的风速；如图11所示，四个网格点对应的坐标分别为P1(x1,y1,z1)、P2(x2,y2,z2)、P3(x3,y3,z3)、P4(x4,y4,z4)，初始点P处的坐标为P(x,y,z)；

S2012’：分别计算四个网格点到初始点P的距离，如图11所示，计算公式如下：

其中，x、y、z表示初始点P的坐标，xi、yi、zi表示第i个网格点的坐标，di表示第i个网格点到初始点P的距离；

S2013’：分别计算四个风速对初始速度V(P，t)的权重，计算公式如下：

其中，di表示第i个网格点到初始点P的距离，k为影响权重的预设值，具体的，k的取值大于0，其表示距离的幂，k的取值直接影响权重Wi的值，进而影响初始速度V(P，t)；当k值选择较小(即接近0)时，此时初始点P的初始速度V(P，t)近似于四个风速V(P1)、V(P2)、V(P3)、V(P4)的平均值，当k取值较大(如大于10)时，则初始点P的初始速度V(P，t)将接近于距离初始点P最近的网格点的风速。在本实施例中，k选择标准参考最小平均绝对偏差，即k的取值为2。

S2014’：根据风速以及权重计算初始速度V(P，t)，计算公式如下：

其中，V(pi)表示第i个网格点处的风速。

由于与初始点P距离不同的网格点处的风速对污染物的运行速度的影响不同，故通过设置权重来计算初始速度V(P，t)，计算结果更为精确。

S202：根据初始点P的位置坐标P(t)(即上述P(x,y,z))和初始速度V(P，t)计算预测点P1的位置坐标P'(t+Δt)，计算公式为：

P'(t+Δt)＝P(t)+V(P，t)*Δt

其中P'(t+Δt)为根据实际观测数值P(t)、经过Δt时间计算得到的预测点P1的位置坐标，P(t)为污染物在初始点P的位置坐标，V(P，t)为污染物在初始点P、时刻为t时的初始速度，Δt为预测点P1与初始点P的时间差，其中Δt的取值满足：

Umax*Δt<0.75，

其中Umax为NetCDF文件中所有网格点对应的风速的绝对值的最大值，Umax的单位为网格/秒，Δt的单位为秒，0.75表示单位网格的0.75倍，理论上，Δt的取值越小，得到的坐标点越多，计算得到的后向轨迹越为精确，但数据处理量也越多，而通过该公式计算得到的Δt的取值在计算预测点P1的坐标时与污染物的实际坐标非常接近，准确度较高，且数据处理量维持在合理范围内。

S203：根据初始点P和预测点P1的位置坐标对预测点P1进行修正以得到修正点的位置坐标，修正公式为：

P(t+Δt)＝P(t)+0.5[V(P，t)+V(P'，t+Δt)]*Δt

初始点P、预测点P1以及修正点的位置关系如图12所示，即修正点的坐标P(t+Δt)等于初始点的坐标P(t)加上二分之一初始点和预测点坐标的矢量和，其中P(t+Δt)为对预测点P1进行修正后得到的修正点的位置坐标，V(P'，t+Δt)为污染物在预测点P1处的运行速度，该运行速度V(P'，t+Δt)的获取与初始速度V(P，t)的获取过程一致，此处不再详细说明；

S204：将修正点作为初始点P重复步骤S201至步骤S203以得到污染物的多个位置坐标，将多个位置坐标连接成直线后得到污染物的后向轨迹，如图13所示，其中图13中的虚线为进行修正时使用的线段，带箭头的线段表示污染物的后向轨迹；

S205：将后向轨迹显示于GIS地图中；其中后向轨迹上的每一个点都具有一个三维坐标，该坐标可以指示污染物在大气中的经度、纬度以及垂直地面的高度，从而可以在GIS地图中展示污染物的后向轨迹，更生动形象，方便观察。

在获取后向轨迹的过程中，计算初始速度时采用的数据来自MM5模式生成的NetCDF文件，而NetCDF文件中存储的为3*3KM的网格数据，因此该数据能很好的模拟局部气流的运行状态，同时，对计算得到的预测点进行修正并将修正点作为污染物的位置坐标，计算结果更为精确，从而得到的后向轨迹与污染物的实际运行轨迹更为相近。

再请参考图14，本发明基于后向轨迹的污染源追踪系统100包括：

第一参数计算模块10，用于根据污染物的后向轨迹依次计算后向轨迹经过的每一网格内污染物的输送概率、距离权重以及时间权重；

第二参数计算模块11，用于根据污染物的源清单计算网格内污染物的排放强度；

输送强度计算模块12，用于根据输送概率、距离权重、时间权重以及排放强度计算网格内污染物的输送强度，计算公式为：

T_l(i,j)＝R_l(i,j)*E(i,j)*W_d(i,j)*W_tl(i,j)

其中，T_l(i,j)为网格(i,j)内污染物的输送强度，R_l(i,j)为网格(i,j)内污染物的输送概率，E(i,j)为网格(i,j)内污染物的排放强度，W_d(i,j)为网格(i,j)内污染物的距离权重，W_tl(i,j)为网格(i,j)内污染物的时间权重；

总输送强度计算模块13，用于根据各网格内的输送强度计算后向轨迹上污染物的总输送强度；

预测区域获取模块14，用于将污染物在不同时刻的多个后向轨迹对应的多个总输送强度最大的前N个后向轨迹经过的网格区域作为预测区域，其中N为正整数，优选的，取1个月内的720条后向轨迹对应的720个总输送强度最大的前N个后向轨迹经过的网格区域作为预测区域，其中1个月内每隔1小时为污染物绘制一条后向轨迹；以及

地理位置确定模块15，用于根据预测区域以及源清单确定污染源的地理位置；以及

后向轨迹获取模块16，用于获取污染物的后向轨迹。

具体的，如图15所示，第一参数计算模块10包括：

输送概率计算单元101，用于根据污染物的后向轨迹计算网格内污染物的距离权重：R_(i,j)＝τ_l(i,j)/n

其中τ_l(i,j)为后向轨迹在网格(i,j)内的停留时间，n为后向轨迹的运行总时间；

距离权重计算单元102，用于根据污染物的后向轨迹计算网格内污染物的距离权重：W_d(i,j)＝1/(d(i,j)/5+1)

其中d(i,j)为网格(i,j)与后向轨迹上污染物的目的地之间的距离；以及

时间权重计算单元103，用于根据污染物的后向轨迹计算网格内污染物的时间权重：W_tl(i,j)＝1/((t_l(i,j))/18+1)

其中t_l(i,j)为污染物由网格(i,j)移动至后向轨迹上污染物的目的地所需的时间。

具体的，如图16所示，第二参数计算模块11包括:

读取单元111，用于从源清单中读取网格以及3个相邻网格内污染物在不同高度层处的排放强度；

总排放强度计算单元112，用于依次计算每一网格内污染物的总排放强度，总排放强度为各个高度层相应的排放强度的总和；以及

平均值计算单元113，用于将网格以及相邻网格的总排放强度取平均值以作为网格内污染物的排放强度。

优选的，后向轨迹获取模块16包括：

初始速度获取单元，用于从MM5模式生成的NetCDF文件中获取污染物在初始点的初始速度；

预测点位置计算单元，用于根据初始点的位置坐标和初始速度计算预测点的位置坐标，计算公式为：P'(t+Δt)＝P(t)+V(P，t)*Δt；

预测点修正单元，用于根据初始点和预测点的位置坐标对预测点进行修正，修正公式为：P(t+Δt)＝P(t)+0.5[V(P，t)+V(P'，t+Δt)]*Δt，其中Δt的取值满足：Umax*Δt<0.75，其中Umax为NetCDF文件中所有网格点的风速的绝对值的最大值；

后向轨迹绘制单元，用于根据获取到的污染物的多个位置坐标进行绘制以得到污染物的后向轨迹；

其中，P'(t+Δt)为预测点的位置坐标，P(t)为初始点的位置坐标，V(P，t)为初始速度，Δt为预测点与初始点的时间差，P(t+Δt)为对预测点进行修正后得到的修正点的位置坐标，V(P'，t+Δt)为污染物在预测点的运行速度。

其中，初始速度获取单元包括：

第一风速读取子单元，用于从NetCDF文件中读取初始点周围四个网格点处的风速；

第一初始速度计算子单元，用于计算四个风速的矢量平均值以作为初始速度。在其他实施例中，初始速度获取单元包括：

第二风速读取子单元，用于从NetCDF文件中读取初始点周围四个网格点处的风速；

距离计算子单元，用于分别计算四个网格点到初始点的距离：

权重计算子单元，用于分别计算四个风速对初始速度的权重：

第二初始速度计算子单元，用于根据风速以及权重计算初始速度：

其中，x、y、z表示初始点的坐标，xi、yi、zi表示第i个网格点的坐标，di表示第i个网格点到初始点的距离，Wi表示第i个网格点处的风速对初始速度的权重，V(pi)表示第i个网格点处的风速，V(p)表示初始速度，k为影响权重的预设值。

以上结合最佳实施例对本发明进行了描述，但本发明并不局限于以上揭示的实施例，而应当涵盖各种根据本发明的本质进行的修改、等效组合。

Claims

1.一种基于后向轨迹的污染源追踪方法，其特征在于，包括：

T_l(i,j)＝R_l(i,j)*E(i,j)*W_d(i,j)*W_tl(i,j)

2.如权利要求1所述的基于后向轨迹的污染源追踪方法，其特征在于，步骤S5具体为：

3.如权利要求1所述的基于后向轨迹的污染源追踪方法，其特征在于，根据污染物的后向轨迹计算所述网格内所述污染物的输送概率具体为：

R_(i,j)＝τ_l(i,j)/n

4.如权利要求3所述的基于后向轨迹的污染源追踪方法，其特征在于，所述停留时间以所述后向轨迹在所述网格(i,j)内的轨迹点个数表示，所述运行总时间以所述后向轨迹上的所有轨迹点个数表示。

5.如权利要求1所述的基于后向轨迹的污染源追踪方法，其特征在于，根据污染物的后向轨迹计算所述网格内所述污染物的距离权重具体为：

W_d(i,j)＝1/(d(i,j)/5+1)

6.如权利要求1所述的基于后向轨迹的污染源追踪方法，其特征在于，根据污染物的后向轨迹计算所述网格内所述污染物的时间权重具体为：

W_tl(i,j)＝1/((t_l(i,j))/18+1)

7.如权利要求1所述的基于后向轨迹的污染源追踪方法，其特征在于，根据所述污染物的源清单计算所述网格内所述污染物的排放强度具体为:

8.如权利要求1至7任一项所述的基于后向轨迹的污染源追踪方法，其特征在于，所述后向轨迹的计算方法包括：

(3)根据所述初始点和所述预测点的位置坐标对所述预测点进行修正以得到修正点的位置坐标，修正公式为：P(t+Δt)＝P(t)+0.5[V(P，t)+V(P'，t+Δt)]*Δt；

9.如权利要求8所述的基于后向轨迹的污染源追踪方法，其特征在于，步骤(1)具体包括：

从所述NetCDF文件中读取所述初始点周围四个网格点处的风速；

计算四个所述风速的矢量平均值以作为所述初始速度。

10.如权利要求8所述的基于后向轨迹的污染源追踪方法，其特征在于，步骤(1)具体包括：

分别计算四个所述网格点到所述初始点的距离：

分别计算四个所述风速对所述初始速度的权重：

<mrow> <mi>W</mi> <mi>i</mi> <mo>=</mo> <mfrac> <mrow> <mn>1</mn> <mo>/</mo> <msup> <mi>di</mi> <mi>k</mi> </msup> </mrow> <mrow> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>4</mn> </msubsup> <mn>1</mn> <mo>/</mo> <msup> <mi>di</mi> <mi>k</mi> </msup> </mrow> </mfrac> <mo>;</mo> </mrow>

根据所述风速以及所述权重计算所述初始速度：

<mrow> <mi>V</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>4</mn> </msubsup> <mi>W</mi> <mi>i</mi> <mo>*</mo> <mi>V</mi> <mrow> <mo>(</mo> <mi>P</mi> <mi>i</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

其中，x、y、z表示所述初始点的坐标，xi、yi、zi表示第i个所述网格点的坐标，di表示第i个所述网格点到所述初始点的距离，Wi表示第i个网格点处的风速对所述初始速度的权重，V(pi)表示第i个网格点处的风速，V(p)表示初始速度，k为影响所述权重的预设值。

11.如权利要求8所述的基于后向轨迹的污染源追踪方法，其特征在于，所述Δt的取值满足：Umax*Δt<0.75，其中Umax为所述NetCDF文件中所有网格点对应的风速的绝对值的最大值。

12.一种基于后向轨迹的污染源追踪系统，其特征在于，包括：

T_l(i,j)＝R_l(i,j)*E(i,j)*W_d(i,j)*W_tl(i,j)

13.如权利要求12所述的基于后向轨迹的污染源追踪系统，其特征在于，所述预测区域获取模块具体为：

取1个月内的720条所述后向轨迹对应的720个所述总输送强度最大的前N个所述后向轨迹经过的网格区域作为预测区域，其中1个月内每隔1小时为所述污染物绘制一条所述后向轨迹。

14.如权利要求12所述的基于后向轨迹的污染源追踪系统，其特征在于，所述第一参数计算模块包括：

输送概率计算单元，用于根据污染物的后向轨迹计算所述网格内所述污染物的距离权重：R_(i,j)＝τ_l(i,j)/n

其中τ_l(i,j)为所述后向轨迹在所述网格(i,j)内的停留时间，n为所述后向轨迹的运行总时间；

距离权重计算单元，用于根据污染物的后向轨迹计算所述网格内所述污染物的距离权重：W_d(i,j)＝1/(d(i,j)/5+1)

其中d(i,j)为所述网格(i,j)与所述后向轨迹上所述污染物的目的地之间的距离；以及

时间权重计算单元，用于根据污染物的后向轨迹计算所述网格内所述污染物的时间权重：W_tl(i,j)＝1/((t_l(i,j))/18+1)

15.如权利要求12所述的基于后向轨迹的污染源追踪系统，其特征在于，所述第二参数计算模块包括:

读取单元，用于从所述源清单中读取所述网格以及3个相邻网格内所述污染物在不同高度层处的排放强度；

总排放强度计算单元，用于依次计算每一所述网格内所述污染物的总排放强度，所述总排放强度为各个高度层相应的所述排放强度的总和；以及

平均值计算单元，用于将所述网格以及相邻网格的所述总排放强度取平均值以作为所述网格内所述污染物的排放强度。

16.如权利要求12至15任一项所述的基于后向轨迹的污染源追踪系统，其特征在于，还包括后向轨迹获取模块，所述后向轨迹获取模块包括：

初始速度获取单元，用于从MM5模式生成的NetCDF文件中获取所述污染物在初始点的初始速度；

预测点位置计算单元，用于根据所述初始点的位置坐标和所述初始速度计算预测点的位置坐标，计算公式为：P'(t+Δt)＝P(t)+V(P，t)*Δt；

预测点修正单元，用于根据所述初始点和所述预测点的位置坐标对所述预测点进行修正，修正公式为：P(t+Δt)＝P(t)+0.5[V(P，t)+V(P'，t+Δt)]*Δt；

后向轨迹绘制单元，用于根据获取到的所述污染物的多个位置坐标进行绘制以得到所述污染物的所述后向轨迹；

其中，P'(t+Δt)为所述预测点的位置坐标，P(t)为所述初始点的位置坐标，V(P，t)为所述初始速度，Δt为所述预测点与所述初始点的时间差，P(t+Δt)为对所述预测点进行修正后得到的修正点的位置坐标，V(P'，t+Δt)为所述污染物在所述预测点的运行速度。

17.如权利要求16所述的基于后向轨迹的污染源追踪系统，其特征在于，所述初始速度获取单元包括：

第一风速读取子单元，用于从所述NetCDF文件中读取所述初始点周围四个网格点处的风速；

第一初始速度计算子单元，用于计算四个所述风速的矢量平均值以作为所述初始速度。

18.如权利要求16所述的基于后向轨迹的污染源追踪系统，其特征在于，所述初始速度获取单元包括：

第二风速读取子单元，用于从所述NetCDF文件中读取所述初始点周围四个网格点处的风速；

距离计算子单元，用于分别计算四个所述网格点到所述初始点的距离：

权重计算子单元，用于分别计算四个所述风速对所述初始速度的权重：

第二初始速度计算子单元，用于根据所述风速以及所述权重计算所述初始速度：

19.如权利要求16所述的基于后向轨迹的污染源追踪系统，其特征在于，所述Δt的取值满足：Umax*Δt<0.75，其中Umax为所述NetCDF文件中所有网格点对应的风速的绝对值的最大值。