基于遥感气溶胶和萤火虫群算法的空气污染源识别方法

技术领域

本发明属于空气污染源识别技术领域，更为具体地讲，涉及一种基于遥感 气溶胶和萤火虫群算法的空气污染源识别方法。

背景技术

大气污染是我国目前面临的严重环境污染问题，包括PM2.5在内的大气污 染物不仅降低空气能见度，影响出行，同时有研究表明PM2.5能直接进入人体 下呼吸道，与呼吸道疾病和心脏疾病有密切联系，呼吸道疾病目前在我国呈日 益上升趋势，因而日益受到人们的关注。

防止大气污染在于源头控制，找到大气污染物的排放源，摸清其排放量是 进行下一步管理的必要条件，对于决策者具有重要的意义。

目前，常用的大气污染物溯源的方法有两类：一类是基于地面采样的下推 法，其原理是通过地面实地采样和实验室化学组分分析的方法，通过污染物化 学成分特征分析污染物来源。这类方法中，最具代表性的是受体模型，也是目 前最广泛采用的污染物溯源方法；另一类是基于排放清单和扩散模型的上推法， 其原理是根据获取的排放清单，代入扩散模型中进行污染物运动模拟，其最大 的特点是能进行污染预测。这两类方法各有其优点，也存在明显的不足。对于 上推法，由于需要进行大量的实地采样和实验室化验分析，工作量很大、成本 昂贵，除了在一些重点区域间断性采用外，既难以大范围使用，也难以持续进 行；对于下推法，排放清单的建立是关键，但目前为止，我国还没有完整的排 放清单。另外，由于排放源分类不明确，缺少排放因子数据库等因素，使得构 建排放清单在某些情况下几乎不可能。

随着航天技术的快速发展，地表空间信息获取能力和水平大大提升。卫星 遥感数据在时间分辨率、空间分辨率和光谱分辨率上都有了极大的提升，利用 卫星遥感影像定量化获取地表或大气物理化学性质已被证明是可行的，随着研 究的深入和算法的改进，定量反演的精度不断提升。已有大量的研究利用遥感 影像估算大气污染物浓度以及分析大气污染物分布情况。如何利用成本相对低 廉、空间覆盖范围大、周期性获取的卫星遥感影像解决污染物溯源的问题是现 有技术需要解决的问题。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于遥感气溶胶和萤火 虫群算法的空气污染源识别方法，采用萤火虫群算法对卫星光谱遥感图像反演 得到的气溶胶厚度值进行污染溯源，高效准确地实现空气污染源的识别。

为实现上述发明目的，本发明基于遥感气溶胶和萤火虫群算法的空气污染 源识别方法包括以下步骤：

S1：获取卫星光谱遥感图像和对应的区域数字地图；

S2：获取卫星光谱遥感图像对应区域的风向和风速，得到风速矢量

S3：对卫星光谱遥感图像进行气溶胶光学厚度反演，得到卫星光谱遥感图 像中每个像元的气溶胶光学厚度值；

S4：对卫星光谱遥感图像进行网格化，网格边长根据实际情况设置，将卫 星光谱遥感图像划分为M×N个图像块，计算每个图像块中所包含像元的气溶胶 光学厚度平均值；同时对卫星光谱遥感图像对应的区域数字地图进行网格化， 划分为M×N个图像块，根据企业坐标得到企业k所对应的图像块坐标E_k， k＝1,2,…,K，K表示企业数量；

S5：采用基于GSO的污染溯源方法进行污染物溯源，其具体步骤包括：

S5.1：设置每只萤火虫的初始荧光素值L_i(0)相同，将各像素块的坐标作为 萤火虫的初始位置x_i(0)，将像素块所对应的气溶胶光学厚度平均值作为萤火虫 的属性值A_i，设置各萤火虫邻域初始值i＝1,2,…,Q，Q＝M×N；

S5.2：令迭代次数t＝1；

S5.3：根据以下公式计算第t代萤火虫的荧光素：

L_i(t)＝(1-ρ)L_i(t-1)+γF(x_i(t))

其中，t表示迭代次数，i表示萤火虫序号，L_i(t)表示第t代萤火虫i的荧光 素，L_i(t-1)表示第t-1代萤火虫i的荧光素，ρ表示荧光素的衰减率，ρ∈(0,1)， γ表示荧光素更新率，F(x_i(t))表示萤火虫i在当前位置x_i(t)的目标函数值，目标 函数的计算公式为：

$F\left(x_i\right(t\left)\right)=\underset{k=1}{\overset{K}{\Sigma }}\frac{1}{d_{ik}^2\left(t\right)}$

其中，d_ik(t)表示萤火虫i在当前位置x_i(t)与企业k坐标E_k的距离；

S5.4：对于每只萤火虫i，分别计算向其邻域N_i(t)内另一只萤火虫j的移动 概率p_ij(t)；采用轮盘赌规则选择本次萤火虫i所要向其运动的萤火虫j，将其序 号记为j′，然后根据下式更新萤火虫i的位置：

$x_i(t+1)=x_i\left(t\right)+s\left(\frac{x_{j^{\prime }}\left(t\right)-x_i\left(t\right)}{\left|\right|x_{j^{\prime }}\left(t\right)-x_i\left(t\right)\left|\right|}\right)+{\Delta }_{attibute}+{\Delta }_{wind}$

其中，x_j′(t)表示第t代萤火虫j′的位置，s表示移动步长，||||表示求欧几里 德范数；

Δ_attibute表示萤火虫相似性修正因子，其计算方法为：首先在萤火虫i的邻域 N_i(t)内，搜索与萤火虫i的属性值A_i最接近的萤火虫j^*，然后计算Δ_attibute：

${\Delta }_{attibute}=s_a·\left(\frac{x_{j^{*}}\left(t\right)-x_i\left(t\right)}{\left|\right|x_{j^{*}}\left(t\right)-x_i\left(t\right)\left|\right|}\right)$

其中，s_a表示预设的步长，x_j*(t)是第t代萤火虫j^*的位置；

Δ_wind表示风速风向修正因子，其计算公式为：其中s_w表 示缩放因子；

S5.5：如果迭代次数t＝T，T表示预设的最大迭代次数，将每只萤火虫的当 前位置作为其源头位置停止迭代；否则进入步骤S5.6；

S5.6：更新每只萤火虫的邻域半径然后令t＝t+1，返回步骤S5.3：

S6：根据企业所对应的图像块坐标E_k，得到企业污染覆盖范围的半径R_k；

S7：根据污染覆盖范围的半径R_k与各只萤火虫的源头位置将企业k作 为其污染覆盖范围内的萤火虫所对应的污染物的产生企业，统计其污染覆盖范 围内的萤火虫所对应的气溶胶厚度值水平，作为企业k的污染指标。

本发明基于遥感气溶胶和萤火虫群算法的空气污染源识别方法，首先根据 卫星光谱遥感图像反演得到气溶胶光学厚度值，得到对应区域的风速矢量，对 卫星光谱遥感图像和对应的区域数字地图网格化，得到每个图像块的气溶胶光 学厚度平均值和企业对应的图像块坐标；将每个图像块的坐标作为GSO算法中 萤火虫的初始位置，气溶胶厚度平均值作为萤火虫的属性，然后迭代，在每次 更新萤火虫位置时，引入由属性值计算得到的萤火虫相似性修正因子和风速矢 量计算得到的风速风向修正因子，多次迭代后得到萤火虫的源头位置；然后计 算企业的污染覆盖范围半径，将企业作为其污染覆盖范围内的萤火虫所对应的 污染物的产生企业，实现污染源识别。

本发明具有以下有益效果：

(1)采用萤火虫群算法，准确高效地迭代得到污染物的源头位置；

(2)萤火虫群算法中更新萤火虫位置时，引入由属性值计算得到的萤火虫 相似性修正因子和风速矢量计算得到的风速风向修正因子，可以使萤火虫群算 法适应本发明应用环境的实际情况，使其迭代结果更符合实际，提高污染源识 别的准确度；

(3)本发明只需要远程光谱遥感信息以及企业的地理信息系统去量化每个 企业的污染情况，而不需要像以往那样通过极为昂贵的样本检测和低效的排放 物清单方法，从而高效准确、低成本地实现空气污染源识别。

附图说明

图1是本发明基于遥感气溶胶和萤火虫群算法的空气污染源识别方法的具 体实施方式流程图；

图2是本发明中基于GSO的污染溯源方法的流程图；

图3是成都市区域图；

图4是都江堰区域的气溶胶光学厚度值网格化图像；

图5是新津区域的气溶胶光学厚度值网格化图像；

图6是都江堰区域的污染溯源图像；

图7是新津区域的污染溯源图像。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员 更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和 设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图1是本发明基于遥感气溶胶和萤火虫群算法的空气污染源识别方法的具 体实施方式流程图。如图1所示，本发明基于遥感气溶胶和萤火虫群算法的空 气污染源识别方法包括以下步骤：

S101：获取图像：

首先获取卫星光谱遥感图像和对应的区域数字地图。光谱遥感图像一般采 用多光谱和高光谱遥感图像。区域数字地图中包含需要监测的各个企业的位置 信息。

S102：获取风速矢量：

获取卫星光谱遥感图像对应区域的风向和风速，得到风速矢量本 发明中风速矢量的指向方向为风的来向。

S103：气溶胶光学厚度反演：

对卫星光谱遥感图像进行气溶胶光学厚度反演，得到卫星光谱遥感图像中 每个像元的气溶胶光学厚度值。本发明中将气溶胶光学厚度值作为大气污染物 的量化参照物。目前气溶胶光学厚度反演算法有暗目标法、结构函数法、多角 度偏振法等，可以根据实际需要来进行选择，其具体反演过程在此不再赘述。 一般来说，在反演后的卫星光谱遥感图像中还存在空白区域，采用插值方法即 可得到空白区域的气溶胶光学厚度值。

S104：图像网格化：

对卫星光谱遥感图像进行网格化，网格边长根据实际情况设置，将卫星光 谱遥感图像划分为M×N个图像块，计算每个图像块中所包含像元的气溶胶光学 厚度平均值，记为A(m,n)，其中m＝1,2,…,M，n＝1,2,…,N。因为原始卫星光谱 遥感图像所包含的像元数量较多，网格化的目的主要是为了减少计算量。

同时对卫星光谱遥感图像对应的区域数字地图也进行网格化，同样地，将 区域数字地图也划分为M×N个图像块，根据企业坐标得到企业k所对应的图像 块坐标E_k，k＝1,2,…,K，K表示企业数量。

S105：基于GSO的污染溯源：

GSO即萤火虫群优化算法(GlowwormSwarmOptimization)是由 krishnanand和Ghose提出来的一种新颖的自然启发式优化算法，GSO是第一 个使用目标函数值相等或不等的多模态函数优化算法。在GSO中，每只萤火虫 被视为解空间的一个解，每只人工萤火虫携带了一定量的荧光素，发出一定量 的荧光并拥有各自的视野，即决策域。荧光素的大小和自己所在位置的目标函 数有关，越亮的萤火虫表示它所在的位置越好，所对应的目标函数值也更优。 每只萤火虫向邻域内较亮的萤火虫移动，最终，使得萤火虫聚集在较优的位置 上，找到多个极值点，从而达到种群寻优的目的。GSO算法是目前一种常用的 寻优算法，其具体流程在此不再赘述。

在原始的GSO算法，一只萤火虫的运动取决于荧光素，而荧光素由目标函 数值和萤火虫位置决定，原始算法没有考虑每只萤火虫的属性。然而在本发明 的具体应用中，将每个图像块的气溶胶厚度值A(m,n)作为萤火虫的属性，该属 性是影响萤火虫运动的一个重要因素，因为有着相同气溶胶水平的点更有可能 来自同一个污染源。此外，风也是一个重要因素。这是因为风影响了污染物的 扩散，那么在进行污染溯源的时候，就需要考虑风力的影响。基于这两方面的 考虑，本发明在传统GSO算法的基础上提出了一种基于GSO的污染溯源方法。

图2是本发明中基于GSO的污染溯源方法的流程图。如图2所示，基于 GSO的污染溯源方法具体步骤包括：

S201：初始化萤火虫：

与传统GSO算法中随机初始化萤火虫不同，本发明将每个像素块作为一只 萤火虫，设置每只萤火虫的初始荧光素值L_i(0)相同，将各像素块的坐标(m,n)作 为萤火虫的初始位置x_i(0)，将像素块所对应的气溶胶光学厚度平均值A(m,n)作 为萤火虫的属性值A_i。此外，还要设置各萤火虫邻域初始值以及其他GSO 算法需要的参数。显然本发明中萤火虫的数量Q＝M×N，i＝1,2,…,Q。

S202：令迭代次数t＝1；

S203：荧光素更新：

根据以下公式计算第t代萤火虫的荧光素：

L_i(t)＝(1-ρ)L_i(t-1)+γF(x_i(t))

其中，t表示迭代次数，i表示萤火虫序号，L_i(t)表示第t代萤火虫i的荧光 素，L_i(t-1)表示第t-1代萤火虫i的荧光素，ρ表示荧光素的衰减率，ρ∈(0,1)， γ表示荧光素更新率，F(x_i(t))表示萤火虫i在当前位置x_i(t)的目标函数值，目标 函数的计算公式为：

$F\left(x_i\right(t\left)\right)=\underset{k=1}{\overset{K}{\Sigma }}\frac{1}{d_{ik}^2\left(t\right)}$

其中，d_ik(t)表示萤火虫i在当前位置x_i(t)与企业k坐标E_k的距离，可以表示 为d_ik(t)＝||x_i(t)-E_k||。

S204：萤火虫移动：

对于每只萤火虫i，分别计算向其邻域N_i(t)内另一只萤火虫j的移动概率 p_ij(t)：

$p_{ij}\left(t\right)=\frac{L_j\left(t\right)-L_i\left(t\right)}{{\Sigma }_{k\in N_i\left(t\right)}[L_k\left(t\right)-L_i\left(t\right)]}$

其中，邻域集$N_i\left(t\right)=\{j:d_{ij}\left(t\right)<r_d^i\left(t\right);L_i\left(t\right)<L_j\left(t\right)\},$d_ij(t)表示第t代萤火虫i与 萤火虫j当前位置的距离，表示第t代萤火虫i的邻域半径，r_s表 示萤火虫的感知半径，L_j(t)表示第t代萤火虫j的荧光素。

采用轮盘赌规则选择本次萤火虫i所要向其运动的萤火虫j，将其序号记为 j′，然后根据下式更新萤火虫i的位置：

$x_i\left(t+1\right)=x_i\left(t\right)+s\left(\frac{x_{j^{\prime }}\left(t\right)-x_i\left(t\right)}{\left|\right|x_{j^{\prime }}\left(t\right)-x_i\left(t\right)\left|\right|}\right)+{\Delta }_{attibute}+{\Delta }_{wind}$

其中，x_j′(t)表示第t代萤火虫j′的位置，s表示移动步长，||||表示求欧几里 德范数。

Δ_attibute是本发明中考虑萤火虫属性所引入的萤火虫相似性修正因子，代表了 由两只萤火虫相似属性决定的位移，其计算方法为：首先在萤火虫i的邻域N_i(t) 内，搜索与萤火虫i的属性值A_i最接近的萤火虫j^*，即j^*＝{j:minA_j-A_i|}，然后 计算萤火虫相似性修正因子：

${\Delta }_{attibute}=s_a·\left(\frac{x_{j^{*}}\left(t\right)-x_i\left(t\right)}{\left|\right|x_{j^{*}}\left(t\right)-x_i\left(t\right)\left|\right|}\right)$

其中，s_a表示预设的步长，x_j*(t)是第t代萤火虫j^*的位置。

Δ_wind是本发明中考虑风力影响所引入的风速风向修正因子，代表由风导致 的移动距离，其计算公式为：

${\Delta }_{\mathrm{wind}}=s_w·{\overrightarrow{\mathrm{speed}}}_{\mathrm{wind}}$

其中，s_w表示缩放因子，即风速单位与网格化光谱遥感图像分辨率的比率。 例如当网格化后光谱遥感图像的分辨率为5km，即每个网格的边长为5km，风 速单位为m，则缩放因子s_w＝0.0002。

S205：判断是否迭代次数t＝T，T表示预设的最大迭代次数，如果是，进 入步骤S206，否则进入步骤S207。

S206：得到萤火虫源头位置：

将每只萤火虫的当前位置作为其源头位置停止迭代。

S207：更新邻域半径：

在GSO中，每只萤火虫的邻域半径不是恒定的，在每次迭代都会更新。 其更新公式为：

$r_d^i(t+1)=min\{r_s,max\{0,r_d^i\left(t\right)+\beta (n_t-|N_i\left(t\right)\left|\right)\left\}\right\}$

其中，β是常量参数，用于控制更新邻域的速度，n_t表示领域数量控制参数， |N_i(t)|表示邻域集中萤火虫的数量。

S208：令t＝t+1，返回步骤S203。

根据基于GSO的污染溯源方法可以看出，本发明考虑了气溶胶厚度值的相 似度以及风力作用，借助于GSO算法，使代表污染物的萤火虫逐步向其源头移 动，从而实现对污染溯源。

S106：计算企业污染覆盖范围：

根据企业所对应的图像块坐标E_k，得到企业污染覆盖范围的半径R_k。本实 施例中所采用的计算方法为：对于企业k，首先搜索与其距离最近的企业k′，计 算这两个企业之间距离D_kk′，R_k＝0.5D_kk′。

S107：判定污染物产生企业：

根据污染覆盖范围的半径R_k与每只萤火虫的源头位置将企业k作为其 污染覆盖范围内的萤火虫所对应的污染物的产生企业，统计其污染覆盖范围内 的萤火虫所对应的气溶胶厚度值水平，作为企业k的污染指标。根据污染指标， 就可以对该企业的污染情况进行评价。本实施例中所采用的污染指标有三种： 污染总量、污染强度和面积归一化污染。

污染总量表示污染源企业周围的污染点的气溶胶厚度值总和。企业k所对应 的污染总量PG_k的计算公式为：

${\mathrm{PG}}_k=\underset{w=1}{\overset{W_k}{\Sigma }}{x}_w^0$

其中，W_k表示企业k污染覆盖范围内的萤火虫数量，表示企业k污染覆 盖范围内萤火虫w的气溶胶厚度值。

污染强度表示污染源企业周围的污染点的气溶胶厚度值平均值。企业k所对 应的污染强度PI_k的计算公式为：

${\mathrm{PI}}_k=\frac{1}{W_k}\underset{w=1}{\overset{W_k}{\Sigma }}{x}_w^0=\frac{1}{W_k}{\mathrm{PG}}_k$

考虑到企业的分布是不均匀的，有些企业是独立的，有些是聚集的。有时 孤立的企业有可能吸引污染点。为了避免这种影响，提出了面积归一化污染指 标。企业k所对应的面积归一化污染ANP_k的计算公式为：

${\mathrm{ANP}}_k=\frac{1}{{\mathrm{area}}_k}\underset{w=1}{\overset{W_k}{\Sigma }}{x}_w^0=\frac{1}{{\mathrm{area}}_k}{\mathrm{PG}}_k$

其中，area_k表示企业k污染覆盖范围的面积，即

为了说明本发明的技术效果，选择成都市对本发明进行了实验验证。由于 其显著的气候特征(多云多雾，日照时间短，空气湿度高)和封闭式的地理环 境(位于四川盆地的中间)，空气污染物不能快速地散开，很容易聚集在城市和 郊区。图3是成都市区域图。如图3所示，本次实验验证中选择了成都都江堰 (图3(b))和新津(图3(c))两个城市作为实验区域。根据GIS数据的更新日 期，获取2009.3.17、基于6MODISL1B产品(mod021km)、h04/V30的遥 感数据作为实验数据。研究区域的企业的GIS信息采用2009年的信息。

首先采用常规的气溶胶反演算法DDV算法来反演得到气溶胶光学厚度 值，检索气溶胶的分辨率为1km，然后进行网格化。企业的位置可以由GIS矢 量点转换为TIF格式的栅格点后得到。图4是都江堰区域的气溶胶光学厚度值 网格化图像。图5是新津区域的气溶胶光学厚度值网格化图像。如图4和图5 所示，每个小圆点表示了一个气溶胶光学厚度值，每个小方格代表了一个企业。

然后采用GSO算法进行污染溯源。GSO算法参数设置如下：各萤火虫的荧 光素初始值均为L_i(0)＝2，感知范围r_s＝2，自适应决策域初始值邻域 阈值n_t＝2，参数β＝0.08，移动步长s＝s_a＝0.03，荧光素衰减率ρ＝0.2，荧光素 更新率γ＝0.6，最大迭代次数T＝200。

然后计算各个企业的污染覆盖范围，然后得到其污染覆盖范围内的萤火虫 所对应的污染物，作为其源头。

图6是都江堰区域的污染溯源图像。图7是新津区域的污染溯源图像。如 图6和图7所示，经过基于GSO的污染溯源，污染物向各个企业进行了集中， 从而可以很容易对污染物的源头进行判断。

然后统计其污染覆盖范围内的萤火虫所对应的气溶胶厚度值水平，计算各 个企业的污染指标。表1是都江堰区域11家企业的污染指标。

表1

表2是新津区域11家企业的污染指标。

表2

如表1和表2可知，都江堰区域的企业A11的污染程度最高，其污染总量、 污染强度最高和面积归一化污染都是最高的，新津区域的企业B10的污染程度 最高。可见，设置三种污染指标的阈值，根据本发明计算得到的各企业对应的 污染指标，可以很容易地对各个企业的污染程度进行判别，筛选出重点污染企 业。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域 的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对 本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定 的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发 明创造均在保护之列。