CN110427533B - 基于时序粒子跟踪的污染传播模式可视分析方法及系统 - Google Patents

Abstract

本发明提出了基于时序粒子跟踪的污染传播模式可视分析方法及系统。首先综合考虑气象因素采用Runge‑Kutta法跟踪空气污染物粒子的运动轨迹，构建空气污染传播模型，通过MIC相关性、影响频数、输送均值对污染源进行检查。然后采用SF‑kmediods算法对目标站点的空气质量特征向量进行聚类提取传播模式。最后设计并实现一个可视分析系统包括地图、热力图、雷达图等视图，从多个视角探索空气污染传播模式，辅助分析者制定更加可靠的空气联合治理方案。

Description

基于时序粒子跟踪的污染传播模式可视分析方法及系统

技术领域

本发明涉及基于时序粒子跟踪的污染传播模式可视分析方法及系统。

背景技术

随着经济、工业化的快速发展，中国的空气质量迅速恶化。空气污染已经成为中国最严重的环境问题之一，并且越来越受到公众和科学界的关注和重视。空气污染会导致各种呼吸疾病，增加交通事故风险，降低城市可持续发展的竞争优势。为了改善空气质量，研究者们对空气污染进行了各种各样的研究，国家和地方政府也制定了相关的法律法规。然而，由于真正污染源还没有被准确地检测到，政策的时空针对性不够，污染情况仍然很严重。城市空气质量的影响因素通常非常复杂，许多气象学家认为气象因素(如风向、湿度)是一个重要的原因，如果不考虑城市周边地区的气象属性，就不可能对城市空气污染进行有效的分析。目前，城市间空气质量的影响主要通过相关性来度量。然而，这些传统的研究往往忽略了气象因素对空气污染传播的影响，而空气污染的传播与风速、风向和湿度等气象因素紧密相关。不科学的建模方法甚至导致该方法的正确性以及结论的准确性受到大气科学领域的专家质疑。迄今为止，仍然缺少一种有效的方法识别具有重大影响的污染源及空气污染远程输送模式。

发明内容

(一)解决的技术问题

针对现有技术的不足，本发明提供了基于时序粒子跟踪的污染传播模式可视分析方法及系统，基于大规模的多源异构时空数据，提出了一种新的空气污染传播模型构建方法，设计并实现了一个全面的可视分析系统，以分析空气污染传播问题。针对多源异构数据带来的数据规模较大、关系模式复杂难于分析、空间分辨率不一致、时变性等问题，给出了解决方案，对空气污染传播模式进行有效的可视分析，能够帮助分析者直观、快速、全面地了解潜在污染源、复杂的空气污染传播模式，从而为分析者制定空气污染联合治理方案提供辅助，促进环境、经济可持续发展。

(二)技术方案

为实现上述目的，本发明提供如下技术方案：基于时序粒子跟踪的污染传播模式可视分析方法,包括以下步骤：

S1：数据收集与处理，采集气象数据和空气质量数据，并对其进行预处理；

S2：污染输送量化，将预处理过的数据与气象因素进行融合，采用高精度的Runge-Kutta法跟踪污染物粒子的运动，从而量化站点间的污染输送水平，进而检测重要污染源；

S3：污染源站点群检测，基于站点间时序上的污染输送信息，对于每一个目标站点，进行污染源检测；

S4：污染源站点群可视化，基于污染源站点群的检测结果，设计基于热力图和字形的地图可视化污染源站点群；

S5：空气污染传播模式提取，基于检测到的污染源站点群，构建时序特征向量，以表征目标站点空气质量及站点群对其输送的污染物信息，采用SF-kmediods聚类算法对时序上的特征向量进行时间聚类，提取空气污染传播模式；

S6：空气污染传播模式可视化，针对提取的空气污染传播模式提出两个指数Index1和Index2以衡量每个模式的重要性，并通过设计可视化视图交互式探索空气污染传播模式的时空多维特征。

优选的，所述步骤S1还包括以下子步骤：

S11：空气质量数据，通过PM25.in网站获取了中国375个城市1497个空气质量监测站点的数据，每条数据记录监测六种空气污染物的浓度：PM2.5、PM10、NO₂、SO₂、O₃、CO。

S12：采集气象数据，为了正确量化周边地区对监测站点所在地区的污染输送关系，拟构建空气污染传播模型来跟踪污染物的扩散，而气象属性是影响污染扩散的重要因素，从欧洲中期天气预报中心下载了全球的气象要素数据，包含比湿、风场U分量、风场V分量三种属性，风速可由风场U分量、风场V分量得到，计算公式为：

S13：数据时空插值，由于空气质量监测站点停电、传感器损坏等原因，采集的空气质量数据存在缺失问题。因此，采用在空间数据插值中广泛应用的反距离加权插值法填补缺失值。

插值的关键是合理地计算时空距离，设置一个缩放时间值的因子，因子计算方式如下：

其中X_max和Y_max表示插值站所在城市中所有站的最大经度和最大纬度，X_min和Y_min分别表示相应的最小值，t_num＝24表示时间范围。

将时空距离d_i定义为如下公式所示:

其中X_i和Y_i分别表示待插值站点的邻域站点的经度和纬度，x和y分别表示待插值站点的经度和纬度，t_i表示待插值站点i具有有效值的距离的最近时刻，t表示站点i插值的时刻，两个时间戳之间的相减等于小时的差。

对于IDW插值，将同一城市中的其他站视为插值站的邻域点，它需要插值点与其邻域点之间的距离来计算加权平均值，插值计算为：

其中V_i表示插值点的估计值，n表示插值点的邻域点数，V_j表示邻域点j的观测值，dj表示邻域点j与插值点之间的欧式距离。

S14：IAQI计算与空气质量等级划分，不同空气污染物，浓度单位不同，为了在同一可视空间中直观地比较各种污染物的信息，根据中华人民共和国国家环境保护标准《环境空气质量(AQI)技术规定》，计算6种污染物的空气质量分指数IAQI。进一步地，采用空气质量指数AQI衡量空气质量的好坏，其值为各污染物的IAQI中的最大值，范围为0-500，可划分为6个等级。

优选的，所述步骤S2还包括以下子步骤：

S21：污染输送量化，为了跟踪站点的空气污染物的扩散，通过双线性插值算法对站点的气象属性进行插值，同理，当粒子开始运动之后，其到达的位置是不确定的，即不一定在网格点上，在这种情况下，也需要进行气象属性插值。气象数据的网格点和空气质量数据的站点都是由经纬度表示，为了进行时序粒子轨迹跟踪，根据高斯投影把经纬度转换到平面直角坐标。当粒子运动到某一地点，为了进行气象属性插值，根据高斯投影把平面直角坐标转换成经纬度，由于风场是极其不稳定的速度场，为了提高粒子跟踪的准确性，根据Runge-Kutta法，以粒子所在位置s当前时间步t_i和下一个时间步t_i+1的速度均值作为粒子的实际运动速度，如下公式所示：

其中，

分别表示在时间步t_i，s处的风场U分量和风场V分量，

分别表示在时间步t_i，粒子从s处出发的实际风场U分量和风场V分量，则可根据步骤S12中的公式计算粒子的运动速率

并对其进行归一化得到

进一步地，根据

和湿度

计算粒子当前时刻的衰减因子，如下公式所示：

分别以站点s当前时刻的AQI和6种空气污染物的IAQI为粒子的初始生命值，则粒子从t_i-1运动t时刻时的现存生命值计算如下公式所示：

其中的

可以分别表示站点s处在t_i时刻的AQI和6中空气污染物的IAQI中的任一种，即AQI、IAQI_PM2.5、IAQI_PM10、

IAQI_CO、

基于粒子实际运动速度，可根据如下公式计算粒子的运动轨迹：

x_i＝x_i-1+V_i ^s′t                                (8)

优选的，所述步骤S3还包括以下子步骤：

S31：对于每种污染物，计算每个社区站点与目标站点之间的相关性，及该站点对目标站点的污染物输送均值，取每种污染物相关性较大且输送均值较大的站点群的并集，即强影响站点群，主要指那些对目标站点影响频次较多、至少有一种污染物输送相关性较强且输送均值较大的站点。关于相关性的计算，采用一个双变量关系的相关性性度量：最大信息指数MIC。MIC属于一类较大的基于最大信息的非参数勘探统计量，用于识别和分类关系，MIC能够捕获多种相关性关系，包括函数关联和非函数关联。

优选的，所述步骤S4还包括以下子步骤：

S41：污染源站点群可视化，在小倍数地图上，基于六边形的字形设计的地图，六边形每个顶点的圆分别表示该站点对目标站点的六种空气污染物PM2.5、PM10、NO₂、SO₂、O₃、CO。

优选的所述步骤S5还包括以下子步骤：

S51：构建特征向量，基于检测到的污染源站点群，收集其在时序上对目标站点的污染输送信息，并构建污染输送特征向量。

S52：模式提取，基于构建的时序特征向量，通过准确度高且效率好的SF-kmediods算法进行聚类，以提取污染传播模式。SF-kmediods算法计算某点到其它所有点的距离之和最小的点，通过距离之和最短的计算方式，可以减少离群数据对聚类过程的影响，从而使得最终效果更接近真实划分。SF-kmediod算法的难点是聚类个数k的确定，本发明采用手肘法，其核心指标是误差平方和SSE，公式如下：

其中，C_i表示第i个簇，p表示C_i中的样本点，m_i表示C_i的质心(C_i中所有样本的均值)，SSE表示所有样本的聚类误差，表征了聚类效果的好坏。

优选的，所述步骤S6还包括以下子步骤：

S61：模式重要性度量方法，提出两个指数Index1和Index2共同度量每个模式的重要性。基于粒子跟踪的结果，假设目标站点A空气污染指数超标(AQI>150)的超标时间步数为N_A，类C_i包含的时间步数为

类C_i中目标站点空气污染指数超标的时间步数为

可见

则对于类C_i，即模式C_i：

其中，Index1表示出现该模式时，目标站点AQI超标的可能性。Index1越大，说明出现这种模式时，目标站点AQI越容易超标。Index2表示目标站点AQI超标时，出现该模式的可能性。Index2越大，说明这种模式对目标站点AQI超标的贡献率越大。

S62：模式空间分布可视化，将S61中的每种模式设计了一个基于南丁格尔玫瑰图的空间污染输送图，每个南丁格尔玫瑰图表示一个影响站点，站点的位置分布基于其地理位置，同时采用碰撞检测算法以解决图标遮挡问题。通过鼠标交互地悬停于玫瑰图标上可以观察到站点的细节信息，包括站点名称、站点所在城市、站点编号及影响频次。

S63：模式时间分布可视化，统计每个模式中每个月份包含的时间步数，并通过热力图展示模式的时间分布。

S64：模式多变量可视化，了解到模式的时空分布时，可通过雷达图进一步观察目标站点在相应模式下的多变量特征。

基于时序粒子跟踪的污染传播模式可视分析系统，包括：

数据收集与处理模块：用于采集气象数据和空气质量数据，并对其进行预处理；

污染输送量化模块：用于将所述数据收集与处理模块中经过预处理的空气质量数据与气象因素进行融合；

污染源站点群检测模块：基于所述污染输送量化模块中站点间时序上的污染输送信息，对于每一个目标站点，进行污染源检测；

污染源站点群可视化模块，基于所述污染源站点群检测模块中污染源站点群的检测结果，设计基于热力图和字形的地图可视化污染源站点群；

空气污染传播模式提取模块，基于所述污染源站点群可视化模块检测到的污染源站点群，构建时序特征向量，以表征目标站点空气质量及站点群对其输送的污染物信息，采用SF-kmediods聚类算法对时序上的特征向量进行时间聚类，提取空气污染传播模式；

空气污染传播模式可视化模块；针对所述空气污染传播模式提取模块提取的空气污染传播模式提出两个指数Index1和Index2以衡量每个模式的重要性，并通过设计可视化视图交互式探索空气污染传播模式的时空多维特征。

本发明具备以下有益效果：基于大规模的多源异构时空数据，提出了一种新的空气污染传播模型构建方法，设计并实现了一个全面的可视分析系统，以分析空气污染传播问题。针对多源异构数据带来的数据规模较大、关系模式复杂难于分析、空间分辨率不一致、时变性等问题，给出了解决方案。

附图说明

图1为本发明方法框架图；

图2为气象数据经纬网格；

图3为空气质量等级划分及其颜色编码；

图4为相关性与影响值关系散点图；

图5为污染源站点群可视化；

图6为污染源污染输送信息可视化；

图7为基于站点污染输送的字形设计；

图8为AQI超标指数及对其应的重要性表格；

图9为实施例中保定市1056A站点的污染站点群；

图10为实施例中保定市1056A站点10种传播模式时间分布图；

图11为实施例中保定市1056A站点模式10的污染输送分布图；

图12为实施例中不同模式重要指数视图；

图13为实施例中1056A站点的模式属性视图；

图14为本发明的系统流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1：

如图1所示，基于时序粒子跟踪的污染传播模式可视分析方法，包括以下步骤：

S1：采集气象数据和空气质量数据，并对其进行预处理；

S2：将预处理过的数据与气象因素进行融合，采用高精度的Runge-Kutta法跟踪污染物粒子的运动，从而量化站点间的污染输送水平，进而检测重要污染源；

S3：基于站点间时序上的污染输送信息，对于每一个目标站点，进行污染源检测；

S4：基于污染源站点群的检测结果，设计基于热力图和字形的地图可视化污染源站点群；

S5：基于检测到的污染源站点群，构建时序特征向量，以表征目标站点空气质量及站点群对其输送的污染物信息，采用SF-kmediods聚类算法对时序上的特征向量进行时间聚类，提取空气污染传播模式；

S6：针对提取的空气污染传播模式提出两个指数Index1和Index2以衡量每个模式的重要性，并通过设计可视化视图交互式探索空气污染传播模式的时空多维特征。

需要说明的是，本发明是为了帮助决策者有针对性地进行空气污染联合治理，所以准确检测污染传播模式的污染源、污染时间、污染物类型是至关重要的。

步骤S1还包括以下子步骤：

S11、空气质量数据：通过PM25.in网站获取了中国375个城市1497个空气质量监测站点的数据，每条数据记录监测六种空气污染物的浓度：PM2.5、PM10、NO₂、SO₂、O₃、CO。

S12、采集气象数据：为了正确量化周边地区对监测站点所在地区的污染输送关系，拟构建空气污染传播模型来跟踪污染物的扩散，而气象属性是影响污染扩散的重要因素，从欧洲中期天气预报中心下载了全球的气象要素数据，包含比湿、风场U分量、风场V分量三种属性，风速可由风场U分量、风场V分量得到，计算公式为：

S13、数据时空插值：由于空气质量监测站点停电、传感器损坏等原因，采集的空气质量数据存在缺失问题。因此，采用在空间数据插值中广泛应用的反距离加权插值法填补缺失值。

插值的关键是合理地计算时空距离，设置一个缩放时间值的因子，因子计算方式如下：

其中X_max和Y_max表示插值站所在城市中所有站的最大经度和最大纬度，X_min和Y_min分别表示相应的最小值，t_num＝24表示时间范围。

将时空距离d_i定义为如下公式所示:

其中x_i和y_i分别表示待插值站点的邻域站点的经度和纬度，x和y分别表示待插值站点的经度和纬度，t_i表示待插值站点i具有有效值的距离的最近时刻，t表示站点i插值的时刻，两个时间戳之间的相减等于小时的差。

对于IDW插值，将同一城市中的其他站视为插值站的邻域点，它需要插值点与其邻域点之间的距离来计算加权平均值，插值计算为：

其中V_i表示插值点的估计值，n表示插值点的邻域点数，V_j表示邻域点j的观测值，d_j表示邻域点j与插值点之间的欧式距离。

S14、IAQI计算与空气质量等级划分：不同空气污染物，浓度单位不同，为了在同一可视空间中直观地比较各种污染物的信息，根据中华人民共和国国家环境保护标准《环境空气质量(AQI)技术规定》，计算6种污染物的空气质量分指数IAQI。进一步地，采用空气质量指数AQI衡量空气质量的好坏，其值为各污染物的IAQI中的最大值，范围为0-500，可划分为6个等级，具体划分如附图3所示。

所述步骤S2还包括以下子步骤：

S21、污染输送量化：为了跟踪站点的空气污染物的扩散，通过双线性插值算法对站点的气象属性进行插值，同理，当粒子开始运动之后，其到达的位置是不确定的，即不一定在网格点上，在这种情况下，也需要进行气象属性插值，气象数据的网格点和空气质量数据的站点都是由经纬度表示(如附图2所示)，为了进行时序粒子轨迹跟踪，根据高斯投影把经纬度转换到平面直角坐标。当粒子运动到某一地点，为了进行气象属性插值，根据高斯投影把平面直角坐标转换成经纬度，由于风场是极其不稳定的速度场，为了提高粒子跟踪的准确性，根据Runge-Kutta法，以粒子所在位置s当前时间步t_i和下一个时间步t_i+1的速度均值作为粒子的实际运动速度，如下公式所示：

其中，

分别表示在时间步t_i，s处的风场U分量和风场V分量，

分别表示在时间步t_i，粒子从s处出发的实际风场U分量和风场V分量，则可根据步骤S12中的公式计算粒子的运动速率

并对其进行归一化得到

进一步地，根据

和湿度

计算粒子当前时刻的衰减因子，如下公式所示：

分别以站点s当前时刻的AQI和6种空气污染物的IAQI为粒子的初始生命值，则粒子从t_i-1运动t时刻时的现存生命值计算如下公式所示：

其中的

可以分别表示站点s处在t_i时刻的AQI和6中空气污染物的IAQI中的任一种，即AQI、IAQI_PM2.5、IAQI_PM10、

IAQI_CO、

基于粒子实际运动速度，可根据如下公式计算粒子的运动轨迹：

x_i＝x_i-1+V_i ^s′t                                   (8)

所述步骤S3还包括以下子步骤：

S31、对于每种污染物：计算每个社区站点与目标站点之间的相关性，及该站点对目标站点的污染物输送均值，取每种污染物相关性较大且输送均值较大的站点群的并集，即强影响站点群，主要指那些对目标站点影响频次较多、至少有一种污染物输送相关性较强且输送均值较大的站点，关于相关性的计算，采用一个双变量关系的相关性度量：最大信息指数MIC，并采用散点图进行可视化，效果如附图4所示。MIC属于一类较大的基于最大信息的非参数勘探统计量，用于识别和分类关系，MIC能够捕获多种相关性关系，包括函数关联和非函数关联。

所述步骤S4还包括以下子步骤：

S41、污染源站点群可视化：采用热力图展示污染源的AQI输送强度，如图5所示，热力效果越明显，表示该站点对目标站点的AQI污染越强。在小倍数地图上，基于六边形的字形设计的地图，可视化污染站点向目标站点输送的各种污染物信息(如附图6所示)。六边形每个顶点的圆分别表示该站点对目标站点的六种空气污染物PM2.5、PM10、NO₂、SO₂、O₃、CO，具体如附图7所示。

所述步骤S5还包括以下子步骤：

S51、构建特征向量：基于检测到的污染源站点群，收集其在时序上对目标站点的污染输送信息，并构建污染输送特征向量。

S52、模式提取：基于构建的时序特征向量，通过准确度高且效率好的SF-kmediods算法进行聚类，以提取污染传播模式。SF-kmediods算法计算某点到其它所有点的距离之和最小的点，通过距离之和最短的计算方式，可以减少离群数据对聚类过程的影响，从而使得最终效果更接近真实划分。SF-kmediod算法的难点是聚类个数k的确定，本发明采用手肘法，其核心指标是误差平方和SSE，公式如下：

其中，C_i表示第i个簇，p表示C_i中的样本点，m_i表示C_i的质心(C_i中所有样本的均值)，SSE表示所有样本的聚类误差，表征了聚类效果的好坏，

所述步骤S6还包括以下子步骤：

S61、模式重要性度量方法：提出两个指数Index1和Index2共同度量每个模式的重要性。基于粒子跟踪的结果，假设目标站点A空气污染指数超标(AQI>150)的超标时间步数为N_A，类C_i包含的时间步数为

类C_i中目标站点空气污染指数超标的时间步数为

可见

则对于类C_i，即模式C_i：

其中，Index1表示出现该模式时，目标站点AQI超标的可能性。Index1越大，说明出现这种模式时，目标站点AQI越容易超标。Index2表示目标站点AQI超标时，出现该模式的可能性。Index2越大，说明这种模式对目标站点AQI超标的贡献率越大。

图8说明了这两种指数及其对应的4种事件，对于事件3，当Index1很小时，由

可知Index2比Index1更小。因此，当Index1很小时，不可能出现Index2很大，即事件3是不可能事件，而事件1是最值得用户关注的模式，因为它对于空气联合治理具有较强的参考意义。

以下，举个例子说明采用两个指数共同度量每个模式的重要性的必要性，假设目标站点AQI超标的时间步数为200个：

模式C₁：总时间步:100个，目标站点超标时间步：80个

模式C₂：总时间步:10个，目标站点超标时间步：8个

当仅仅采用Index1时，难以区分事件1和事件2。而当仅仅采用Index2时，难以发现目标站点AQI超标的事件，因此，采用Index1和Index2来共同度量每个模式的重要性。针对Index1和Index2，如图12所示，采用折线图在一个视图中同时展示这两个指数。

S62、模式空间分布可视化：当关注重要模式C_i时，需要进一步细致地确定目标站点地区的本地污染源，因此如图11所示，将S61中的每种模式设计了一个基于南丁格尔玫瑰图的空间污染输送图，每个南丁格尔玫瑰图表示一个影响站点，站点的位置分布基于其地理位置，同时采用碰撞检测算法以解决图标遮挡问题。通过鼠标交互地悬停于玫瑰图标上可以观察到站点的细节信息，包括站点名称、站点所在城市、站点编号及影响频次；每个玫瑰花瓣编码一种属性，花瓣半径表示在C_i模式下该站点向目标站点输送的对应污染物水平，花蕊的颜色表示在C_i模式下该站点对目标站点的影响频次，颜色越深，影响的频次越高。

S63、模式时间分布可视化：统计每个模式中每个月份包含的时间步数，并通过热力图展示模式的时间分布。

S64、模式多变量可视化：了解到模式的时空分布时，可通过雷达图进一步观察目标站点在相应模式下的多变量特征。

如图9-图13所示，使用本发明提供的基于时序粒子跟踪的污染传播模式可视分析方法对2015年河北保定市莲池区五四东路的监测站点1056A进行分析的过程。

采用Runge-Kutta法跟踪空气污染物粒子的运动轨迹，构建空气污染传播模型，并通过MIC相关性、影响频数以及输送均值检测污染站点群。结合中国所有地区的气象数据中的风场U分量和风场V分量采用Runge-Kutta分别计算AQI和6种空气污染物的实际运动速度；此外利用风速Vel与比湿S计算衰减因子和生命值，若污染物经过站点，则当时的生命值作为其输送值。当污染物生命值衰减为0时，通过高斯投影将其轨迹投影到平面直角坐标。最后取每种污染物MIC相关性较大且输送均值较大的站点群的并集，找出对1056A影响较大的站点群，并通过地图、热力图进行可视化，如图9所示。

采用SF-kmediods算法对目标站点的空气质量特征向量进行聚类提取传播模式。2015年河北保定市莲池区五四东路的监测站点1056A为目标站点对全国其他的1496个站点构建污染输送特征向量，以小时为粒度，因此共有8760个污染输送特征向量。采用SF-kmediods算法对8760个污染输送特征向量进行聚类，采用手肘法确定最佳聚类个数。通过计算监测站点1056A的污染输送特征向量可以划分为10个簇，每个簇则表示一种空气污染传播模式。

分析不同空气污染传播模式的时间、空间分布特征。统计每个模式中每个月份包含的时间步数，并通过热力图展示模式的时间分布。统计每个模式下污染站点群对目标站点的污染输送均值，采用南丁格尔玫瑰图可视化。从附图10中可以发现模式0、1、2、4、5、8主要集中出现在夏天(4-7月)；模式3、6、7、9主要集中出现在冬天(1-2月、12月)，少量出现在3月、10月和11月。此外，模式10主要出现在10月，而在其他月份，特别是1月、12月都没有出现。因此可推测模式10可能是一种异常情况。进一步观察附图11模式10的污染输送分布图，通过系统提供的影响频次过滤功能，发现向站点1056A输送污染物的主要站点((a)中圆圈包含的站点)。对比图(b)和(c)，发现目标站点1056A该模式下PM2.5、PM10和NO₂的浓度较高，O₃、CO、SO₂的浓度较低，而其周边站点输送的各种污染物也呈现相似的趋势。用户可以通过鼠标悬停于站点玫瑰图上方，获取这些站点的相关信息，或者导出为文件。因此，推断保定市空气质量监测站点1056A的空气污染传播模式10，其主要影响站点分布在本市、石家庄市、衡水市、沧州市、德州市、邢台市、聊城市。

进一步地，对上述提取到的传播模式，利用空气污染传播模式的AQI超标指数Index1和Index2度量其重要性。采用折线图进行可视化展示，同时采用雷达图观察目标站点在相应模式下的对变量特征。对于监测站点1056A在B中可以得到10种模式，其重要性指数如图12所示：图中可以发现，模式3、模式8和模式10的Index1都很高，表明出现这三种模式时，目标站点1056A站点AQI超标的可能性大。然而，模式10的Index2很小，表明相对于2015年站点1056A超标频次，该模式出现AQI超标频次较少。进一步查看附图13的雷达图，发现模式10的PM2.5和PM10浓度都比较高。

实施例2：如图14所示，基于时序粒子跟踪的污染传播模式可视分析系统，包括：

数据收集与处理模块：用于采集气象数据和空气质量数据，并对空气质量数据进行预处理；

通过PM25.in网站获取了中国375个城市1497个空气质量监测站点的数据，每条数据记录监测六种空气污染物的浓度：PM2.5、PM10、NO₂、SO₂、O₃、CO；从欧洲中期天气预报中心下载了全球的气象数据，包含比湿、风场U分量、风场V分量三种属性，风速可由风场U分量、风场V分量得到；由于空气质量监测站点停电、传感器损坏等原因，采集的空气质量数据存在缺失问题，因此，采用在空间数据插值中广泛应用的反距离加权插值法填补缺失值。

污染输送量化模块：通过站点间的污染输送量化算法将所述数据收集与处理模块中经过预处理的空气质量数据和气象数据进行融合，具体算法框架如下：

输入：网格气象数据：中国地区的网格经纬度gridLongLat，及相应网格的风场U分量、风场V分量、风速Vel、比湿S；监测站点空气质量数据：站点的经纬度stationLongLat，及相应站点的AQI、IAQI_PM2.5、IAQI_PM10、

IAQI_CO、

粒子运动时间阈值φ，站点邻域半径r。

输出：时序上站点间的各种污染物输送值

算法描述：通过双线性插值方法、高斯投影和Runge-Kutta法，计算站点间的污染输送。

1.对风速Vel数据进行归一化，得到Vel^Nor。

2.对于每个站点s，每个时间步t_i：

2.1根据双线性插值算法，计算站点s的气象属性值；

2.2基于Runge-Kutta法计算粒子的运动速度，如公式(5)所示；

2.3根据公式(6)计算粒子生命值此时的衰减因子；

2.4根据高斯投影把粒子当前位置的经纬度转换成平面直角坐标；

2.5粒子在风场作用下运动，其生命值随着运动时间衰减，如公式(7)所示。

2.6粒子每运动一个时间步，根据公式(7)计算其剩余生命值，并根据公式(8)计算其运动轨迹。如果粒子生命值衰减为零，或者运动一定时间阈值φ，则停止运动；否则，如上述更新粒子的运动速度、生命值、衰减速度，粒子继续运动。

2.7重复2.6直到粒子停止运动。

2.8如果粒子经过站点S_k的邻域(以r为半径的邻域)，则计算粒子运动到S_k邻域时的时间t_j，然后计算此时的粒子生命值作为站点s在时间步t_j对站点S_k的各种污染物输送值。

3.算法结束。

污染源站点群检测模块：基于所述污染输送量化模块中站点间时序上的污染输送信息，对于每一个目标站点，进行污染源检测；

基于站点间时序上的污染输送信息，对于每一个目标站点，进行污染源检测。首先，对目标站点影响次数较少但输送污染物较多可能是一种突发情况，空气污染联合治理政策很难针对这种异常情况制定，而输送污染物较少的也不是关注的重点，因此，过滤掉那些影响次数较少的站点。其次，进行空气联合治理时，分析者更关注那些与目标站点相关性大且负面影响强的站点，如图4右上角所示。如果只考虑一种污染物，则其他相关性大的污染物可能被忽略，从而忽略了一些相关性大的站点，因此，同时考虑多种污染物。综合上述考虑，对于每种污染物，计算每个社区站点与目标站点之间的相关性，及该站点对目标站点的污染物输送均值。最后，取每种污染物相关性较大且输送均值较大的站点群的并集，即强影响站点群，主要指那些对目标站点影响频次较多、至少有一种污染物输送相关性较强且输送均值较大的站点。

污染源站点群可视化模块，基于所述污染源站点群检测模块中污染源站点群的检测结果，设计基于热力图和字形的地图可视化污染源站点群；

基于污染源站点群的检测结果，设计基于热力图和字形的地图可视化污染源站点群，用以帮助客户不同详细程度地了解污染源站点群的污染信息。在较高倍数地图上，采用热力图展示污染源的AQI输送强度，如图5所示，热力效果越明显，表示该站点对目标站点的AQI污染越强。

在小倍数地图上，设计基于字形的地图，可视化污染源站点向目标站点传输的各种污染物信息，如图6所示。基于六边形的字形设计如图7所示，六边形每个顶点的圆分别表示该站点对目标站点的六种空气污染物PM2.5、PM10、NO₂、SO₂、O₃、CO，圆的颜色编码该站点对目标站点对应污染物的输送均值(颜色映射如图7图c渐变色条所示),圆的大小编码该站点对目标站点输送对应污染物的MIC相关性。字形中心的六边形的颜色表示该站点对目标站点的AQI输送均值(颜色映射如图7图c渐变色条所示)，六边形的边框颜色表示该站点对目标站点的影响频次(颜色映射如图7图d渐变色条所示)。中心与顶点连线上的柱形图：柱形高度表示对应属性的输送频率分布，柱形颜色表示对应属性的输送频次分布(颜色映射如图7图d渐变色条所示)。同时考虑频数和频率是出于以下考虑，不同站点对目标站点的影响频数不同，甚至相差很大，而计算得到的频率可能是相同的。如A站点影响次数为10次，SO₂输送小于50的有8次，而B站点影响次数为100次，SO₂输送小于50的有80次，计算得到的频率都是80％。为了区分这种情况，用柱形颜色来编码该输送值的输送频数。为了清晰显示输送信息，翻转柱形图对称展示，如图7图b所示。

空气污染传播模式提取模块，基于所述污染源站点群可视化模块检测到的污染源站点群，构建时序特征向量，以表征目标站点空气质量及站点群对其输送的污染物信息，采用SF-kmediods聚类算法对时序上的特征向量进行时间聚类，提取空气污染传播模式。

空气污染传播模式可视化模块；针对所述空气污染传播模式提取模块提取的空气污染传播模式提出两个指数Index1和Index2以衡量每个模式的重要性，并通过设计可视化视图交互式探索空气污染传播模式的时空多维特征。

本发明的工作原理为：提出了基于时序粒子跟踪的污染传播模式可视分析方法及系统，首先综合考虑气象因素采用Runge-Kutta法跟踪空气污染物粒子的运动轨迹，构建空气污染传播模型。通过MIC相关性、影响频数、输送均值对污染源进行检查，然后采用SF-kmediods算法对目标站点的空气质量特征向量进行聚类提取传播模式，最后设计并实现一个可视分析系统包括地图、热力图、雷达图等视图，从多个视角探索空气污染传播模式，辅助分析者制定更加可靠的空气联合治理方案。

需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (8)

1.基于时序粒子轨迹跟踪的空气污染传播模式可视分析方法,其特征在于，包括以下步骤：

S1：采集气象数据和空气质量数据，并对其进行预处理；

S2：将步骤S1中预处理过的气象数据与空气质量数据进行融合，采用高精度的Runge-Kutta法跟踪污染物粒子的运动，从而量化站点间的污染输送水平；

S3：基于步骤S2中站点间时序上的污染输送信息，对于每一个目标站点，进行污染源检测；

S4：基于步骤S3中污染源站点群的检测结果，设计基于热力图和字形的地图可视化污染源站点群；

S5：基于步骤S4中检测到的污染源站点群，构建时序特征向量，以表征目标站点空气质量及站点群对其输送的污染物信息，采用SF-kmediods聚类算法对时序上的特征向量进行时间聚类，提取空气污染传播模式；

S6：基于步骤S5中提取的空气污染传播模式提出两个指数Index1和Index2以衡量每个模式的重要性，并通过设计可视化视图交互式探索空气污染传播模式的时空多维特征。

2.根据权利要求1所述的基于时序粒子轨迹跟踪的空气污染传播模式可视分析方法，其特征在于，所述步骤S1还包括以下子步骤：

S11、空气质量数据：获取空气质量监测站点的数据，每条数据记录监测六种空气污染物的浓度：PM2.5、PM10、NO₂、SO₂、O₃、CO；

S12、采集气象数据：从天气预报中心下载了全球的气象要素数据，包含比湿、风场U分量、风场V分量三种属性，风速由风场U分量、风场V分量得到，计算公式为：

S13、数据时空插值：计算时空距离，设置一个缩放时间值的因子，因子计算方式如下：

其中X_max和Y_max表示插值站所在城市中所有站的最大经度和最大纬度，X_min和Y_min分别表示相应的最小值，t_num＝24表示时间范围，

将时空距离d_i定义为如下公式所示：

其中x_i和y_i分别表示待插值站点的邻域站点的经度和纬度，x和y分别表示待插值站点的经度和纬度，t_i表示待插值站点i具有有效值的距离的最近时刻，t表示站点i插值的时刻，两个时间戳之间相减等于小时的差，

对于IDW插值，将插值站周围的站点视为插值站的邻域点，插值计算为：

其中V_i表示插值点的估计值，n表示插值点的邻域点数，V_j表示邻域点j的观测值，d_j表示邻域点j与插值点之间的欧式距离，

S14：IAQI计算与空气质量等级划分，计算六种空气污染物的空气质量分指数IAQI，采用空气质量指数AQI衡量空气质量的好坏，并划分为6个等级。

3.根据权利要求2所述的基于时序粒子轨迹跟踪的空气污染传播模式可视分析方法，其特征在于，所述步骤S2的具体内容为：

通过双线性插值算法对站点的气象属性进行插值，气象数据的网格点和空气质量数据的站点都是由经纬度表示，根据高斯投影把经纬度转换到平面直角坐标，根据Runge-Kutta法，以粒子所在位置s当前时间步t_i和下一个时间步t_i+1的速度均值作为粒子的实际运动速度，如下公式所示：

其中，

分别表示在时间步t_i，s处的风场U分量和风场V分量，

分别表示在时间步t_i，粒子从s处出发的实际风场U分量和风场V分量，则可根据步骤S12中的公式计算粒子的运动速率

并对其进行归一化得到

根据

和湿度

计算粒子当前时刻的衰减因子，如下公式所示：

分别以站点s当前时刻的AQI和六种空气污染物的IAQI为粒子的初始生命值，则粒子从t_i-1运动t时刻时的现存生命值计算如下公式所示：

其中的

表示站点s处在t_i时刻的AQI和6种空气污染物的IAQA中的任一种，即AQI、IAQA_PM2.5、IAQA_PM10、

IAQA_CO、

基于粒子实际运动速度，可根据如下公式计算粒子的运动轨迹：

x_i＝x_i-1+V_i ^s′t。

4.根据权利要求1所述的基于时序粒子轨迹跟踪的空气污染传播模式可视分析方法，其特征在于，所述步骤S3的具体内容为：

对于每种污染物计算每个社区站点与目标站点之间的相关性，及该社区站点对目标站点的污染物输送均值，取每种污染物相关性较大且输送均值较大的站点群的并集，即强影响站点群，计算相关性采用一个双变量关系的相关性度量：最大信息指数MIC。

5.根据权利要求4所述的基于时序粒子轨迹跟踪的空气污染传播模式可视分析方法，其特征在于，所述步骤S4的具体内容为：

在小倍数地图上，基于六边形的字形设计的地图，六边形每个顶点的圆分别表示该社区站点对目标站点的六种空气污染物PM2.5、PM10、NO₂、SO₂、O₃、CO。

6.根据权利要求1所述的基于时序粒子轨迹跟踪的空气污染传播模式可视分析方法，其特征在于，所述步骤S5还包括以下子步骤：

S51、构建特征向量：基于检测到的污染源站点群，收集其在时序上对目标站点的污染输送信息，并构建污染输送特征向量；

S52、模式提取：基于步骤S51构建的时序特征向量，通过准确度高且效率好的SF-kmediods算法进行聚类，以提取污染传播模式，SF-kmediod算法的难点是聚类个数k的确定，采用手肘法确定k，其核心指标是误差平方和SSE，公式如下：

其中，C_i表示第i个簇，p表示C_i中的样本点，m_i表示C_i的质心，C_i的质心为C_i中所有样本的均值，SSE表示所有样本的聚类误差，表征了聚类效果的好坏。

7.根据权利要求1所述的基于时序粒子轨迹跟踪的空气污染传播模式可视分析方法，其特征在于，所述步骤S6还包括以下子步骤：

S61、模式重要性度量方法：提出两个指数Index1和Index2共同度量每个模式的重要性，基于粒子跟踪的结果，假设目标站点A空气污染指数超标的超标时间步数为N_A，类C_i包含的时间步数为

类C_i中目标站点空气污染指数超标的时间步数为

可见

则对于类C_i，即模式C_i：

Index1表示出现该模式时，目标站点AQI超标的可能性，Index1越大，说明出现这种模式时，目标站点AQI越容易超标，Index2表示目标站点AQI超标时，出现该模式的可能性，Index2越大，说明这种模式对目标站点AQI超标的贡献率越大；

S62、模式空间分布可视化：将S61中的每种模式设计了一个基于南丁格尔玫瑰图的空间污染输送图，每个南丁格尔玫瑰图表示一个影响站点，站点的位置分布基于其地理位置，同时采用碰撞检测算法以解决图标遮挡问题，通过鼠标交互地悬停于玫瑰图标上可以观察到站点的细节信息，包括站点名称、站点所在城市、站点编号及影响频次；

S63、模式时间分布可视化：统计每个模式中每个月份包含的时间步数，并通过热力图展示模式的时间分布；

S64、模式多变量可视化：了解模式的时空分布，通过雷达图进一步观察目标站点在相应模式下的多变量特征。

8.基于时序粒子轨迹跟踪的空气污染传播模式可视分析系统，其特征在于，包括：

数据收集与处理模块：用于采集气象数据和空气质量数据，并对其进行预处理；

污染输送量化模块：用于将所述数据收集与处理模块中经过预处理的空气质量数据与气象因素进行融合；

污染源站点群检测模块：基于所述污染输送量化模块中站点间时序上的污染输送信息，对于每一个目标站点，进行污染源检测；

污染源站点群可视化模块，基于所述污染源站点群检测模块中污染源站点群的检测结果，设计基于热力图和字形的地图可视化污染源站点群；

空气污染传播模式提取模块，基于所述污染源站点群可视化模块检测到的污染源站点群，构建时序特征向量，以表征目标站点空气质量及站点群对其输送的污染物信息，采用SF-kmediods聚类算法对时序上的特征向量进行时间聚类，提取空气污染传播模式；

空气污染传播模式可视化模块；针对所述空气污染传播模式提取模块提取的空气污染传播模式提出两个指数Index1和Index2以衡量每个模式的重要性，并通过设计可视化视图交互式探索空气污染传播模式的时空多维特征。

Images