CN109325616A - 一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测与溯源方法 - Google Patents

Abstract

本发明属于细微颗粒物探测领域，具体涉及一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测与溯源方法，包括以下步骤：对萤火虫算法进行优化，得到萤火虫优化算法；创建模型的样本训练集并为高斯过程回归模型选择协方差函数；利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化；根据优化的超参数得出模型，输出预测。本发明提出一种基于高斯过程回归和萤火虫算法的细颗粒物预测与溯源方法，该方法将高斯过程回归和萤火虫算法相结合，既能细颗粒物空间分布及短时变化态势进行预测，也能对细颗粒物多污染源动态回溯，实现对细颗粒物多污染源的准确定位。

Description

一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测 与溯源方法

技术领域

本发明属于细微颗粒物探测领域，具体涉及一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测与溯源方法。

背景技术

近年来，经济发展带来的环境问题日益凸显，空气质量已成为公众关注的焦点。其中细颗粒物可以深入到人体肺泡和血液循环系统，对人体健康带来极大危害。面对严峻的细颗粒物污染问题，各级政府部门积极采取综合治理措施，并取得一定的成效。然而雾霾治理是一个漫长的过程，相关部门不但需要长期的治理计划，更需要科学化、定量化的预防和实时监测机制，实现标本兼治。对细颗粒物进行预测和溯源，及时有效地发现细颗粒物变化趋势和污染源分布，为有关部门的科学决策和事后追责提供有力支撑，为公众的及时预防提供准确指导，具有很重要的现实意义。

现阶段在细颗粒物预测方面，利用经典的扩散模型如高斯羽流扩散模型或者社区多尺度空气质量模型，可以得到比较准确的结果。但共同的缺点是很难确定应用的条件，很多关键参数的设置都很繁琐，需要的现实环境建模数据输入较多。对于空气中的气体或者细颗粒物等物质的预测应用较多的是人工神经网络模型。然而，面对室外复杂环境下的预测，人工神经网络模型输入等条件的限制可能导致模型产生复杂的网络结构进而导致预测估计时间过长。尽管人工神经网络在非线性建模方面有较强的能力，但面对小样本数据显得有些无力。另一方面，在细颗粒物空间维度分布预测数据需要支撑后续溯源算法的实现，要实现准确的溯源特征，就必须保证数据在突变状态下的敏感性，而人工神经网在处理非平稳数据方面跟踪性误差较大。基于以上考虑选择应用高斯过程回归模型来完成细颗粒物空间和时间维度上的预测，其在处理小样本、非线性等复杂问题时具有很好的效果，相较于人工神经网络具有容易实现、超参数自适应获取、泛化能力强等特点，并且高斯过程回归模型在处理非平稳性问题时，跟踪性更好误差更小。

多模态函数又叫多峰值函数，多模态函数优化问题的一种求解是寻找目标函数全局最优，同时避免局部最优解。另一种优化问题的目标是找到具有相等或不等函数值的多个最优解。

现实环境中雾霾天气的出现，并不是单污染源造成的，其主要来自多污染源多空气污染物复杂作用的结果，通过现有研究可以知道城市雾霾污染的来源主要包括：交通密集区尾气排放、工业密集区废气排放和监测区域以外传入等情况。与源定位相关的大部分工作都致力于寻找、检测和跟踪单个源位置点，当研究多个源位置点时问题就变得复杂化。由于扩散过程由于受到室外环境的多种因素影响，细颗粒物在空间中的分布是动态变化的，这就使得在函数模型创建过程中并不能应用静态的多模态函数建模。而在创建动态气体模型过程中，模型需要设定的环境假设过多，并且对现实环境的空气变化状态有各种要求。对于多模态函数最优值的求解，萤火虫优化算法是一个很好的选择。

综上所述，现有技术存在应用条件不确定、参数设置繁琐、耗时长、误差大，不能处理小样本等问题。

发明内容

本发明是针对空气污染监测领域提出的一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测与溯源方法，旨在对细颗粒物空间分布及短时变化态势进行预测，并对细颗粒物多污染源动态回溯，实现对细颗粒物多污染源的准确定位。

一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测与溯源方法，包括以下步骤：

(1)对萤火虫算法进行优化，得到萤火虫优化算法；

(2)创建模型的样本训练集并为高斯过程回归模型选择协方差函数；

(3)利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化；

(4)根据优化的超参数得出模型，输出预测。

所述对萤火虫算法进行优化，得到萤火虫优化算法，包括：

(1.1)荧光素更新阶段

在荧光素更新阶段期间，每个萤火虫向其之前的荧光素水平增加目标位置在目标函数中的适应性成正比的荧光素量，此外，荧光素值的一部分被减去以模拟荧光素随时间的衰减，荧光素更新规则由下式给出：

l_i(t+1)＝(1-ρ)l_i(t)+γJ(x_i(t+1))

其中，l_i(t)表示与萤火虫i在t时刻的荧光素大小，ρ是荧光素衰减常数(0＜ρ＜1)，γ是荧光素增强常数，J(x_i(t))表示t时i位置处目标函数的值；

(1.2)运动阶段

在运动阶段，每个萤火虫使用概率机制决定朝着一个荧光素值高于其自身的邻居移动，对于每个萤火虫i，向邻居j移动的概率由下式给出：

萤火虫根据上式给出的p_ij(t)选择一个j∈N_i(t)的萤火虫；

其中，j∈N_i(t)，表示时间t处的萤火虫i的邻居的集合，d_ij(t)表示时间t处的萤火虫i和j之间的欧几里得距离，并且表示在时间t与萤火虫i相关的邻域范围

萤火虫运动的离散时间模型为：

其中，x_i(t)∈R^m是在时间t时萤火虫i在m维实空间R^m中的位置，||.||表示欧几里得范数算子，s(s＞0)表示是步长；

(1.3)邻域范围更新阶段

令r0为每个萤火虫的初始邻域范围，为了自适应地更新每个萤火虫的邻域范围，应用以下规则：

其中，β是常数参数，n_t是用于控制邻居数量的参数；

通过下式，计算得到邻域集合中个体与目标萤火虫个体的平均距离：

的范围：在初始化时令当时，令q＝0.96；

(1.4)萤火虫自探索方面

当萤火虫个体在其感知范围内找不到优秀个体时，在此范围内没有找到优秀的个体，那就以现在的位置为起点，进行对数螺旋搜索，使得搜索面积变大，并且在移动过程中移动两圈，如果移动后，新位置的荧光素水平比原来低，而且周围没有优秀的个体，萤火虫会以下式跳出该位置：

x_i(t+1)＝x_i(t)+(1-2rand())·s_i(t)

如果更新后荧光素值仍然小于原来位置水平，则取消位置更新

其中，rand()是(0,1)的随机数。

所述创建模型的样本训练集并为高斯过程回归模型选择协方差函数，包括：

模型训练样本集为：

D＝{(x⁽ⁱ⁾,y⁽ⁱ⁾)|i＝1,2,3,L,n}＝(X,y)

其中，x⁽ⁱ⁾∈R^d为d维输入矢量，表示输入的特征向量，X＝[x⁽¹⁾,x⁽²⁾,L,x⁽ⁿ⁾]表示为d×n维输入矩阵，y⁽ⁱ⁾∈R表示第i个样本对应的输出标量，而y则表示由每个样本对应的输出标量组成的输出矢量，回归问题就要通过训练样本集学习到输入矩阵X与输出矢量y的映射关系：(f(·):R^d→R)，从而得到新的测试点x_*对应的输出f(x_*)

平方指数协方差函数形式如下：

其中，来控制先验方差的，最大允许协方差为当某个函数在y轴的数据范围很大时，其数值相应增大，M＝diag(l²)，其中l∈R^d被称为方差尺度，用来控制协方差衰减率的各向同性长度尺度参数，其中ω_i∈l表示输入的第i个特征的权重，其中集合为超参数。

所述利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化，包括：

(4.1)设置PSO参数：首先选定粒子群大小N，设置最大迭代次数iter_max，精度误差为ε，c₁＝c₂＝2，ω_end＝0.9，ω_ini＝0.4，CR＝0.9,F＝0.9，并初始化每个粒子的x_i＝(x_i1,x_i2,x_i3,L,x_id,L,x_iD)和v_i＝(v_i1,v_i2,v_i3,L,v_id,L,v_iD)；

(4.2)通过下式估计每个粒子的适应度，并找出每个粒子的对应的pbest和gbest；

(4.3)在[0,1]内生成一个随机数r，在[1,D]生成一随机数k，并根据下式对粒子进行变异操作；

Trail_d＝gbest_d+F(pbest_1d+pbest_2d-pbest_3d-pbest_4d)(if r＜CR or d＝k)

其中，CR是交叉常数，缩放因子F是[0,2]内的常数，pbest_1d、pbest_2d、pbest_3d、pbest_4d是从pbest中随机选择的；

(4.4)根据下式更新x_i和v_i；

v_i(t+1)＝ωv_i(t)+c₁r₁[pbest_i-x_i(t)]+c₂r₂[gbest_i-x_i(t)]

x_i(t+1)＝x_i(t)+v_i(t+1)

其中，公式中i＝1,2,3,L,N，N表示粒子个数；r₁,r₂代表两个从U(0,1)中采样的随机序列；c₁,c₂是两个学习因子，c₁＝c₂＝2；ω称为关惯性因子，并且ω≥0

ω(t)＝(ω_ini-ω_end)(G_K-g)/G_k+ω_end

其中G_k为最大迭代次数，ω_ini权值为0.4，ω_end设定为0.9；

(4.5)重复以上过程，达到最大迭代次数或gbest的适应度满足精度误差的要求，该过程将结束。

本发明的有益效果在于：

(1)本发明提出一种基于高斯过程回归和萤火虫算法的细颗粒物预测与溯源方法，该方法将高斯过程回归和萤火虫算法相结合，既能细颗粒物空间分布及短时变化态势进行预测，也能对细颗粒物多污染源动态回溯，实现对细颗粒物多污染源的准确定位。

(2)本发明利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化，提高了细颗粒物的预测精度，并对萤火虫算法在萤火虫移动步长和和萤火虫自探索两方面进行优化设计，提高了萤火虫算法的搜索效率和搜索精度。

附图说明

图1是细颗粒物预测方法架构图；

图2是基于改进PSO算法的改进GPR算法流程图；

图3是改进的萤火虫算法实现流程图；

图4是细颗粒物溯源方法架构图；

具体实施方式

下面结合附图对本发明做进一步描述。

本发明主要是用于细颗粒物的预测与溯源，对细颗粒物空间分布及短时变化态势进行预测，并对细颗粒物多污染源动态回溯，实现对细颗粒物多污染源的准确定位。

近年来，经济发展带来的环境问题日益凸显，空气质量已成为公众关注的焦点。其中细颗粒物(PM2.5)可以深入到人体肺泡和血液循环系统，对人体健康带来极大危害。面对严峻的细颗粒物污染问题，各级政府部门积极采取综合治理措施，并取得一定的成效。然而雾霾治理是一个漫长的过程，相关部门不但需要长期的治理计划，更需要科学化、定量化的预防和实时监测机制，实现标本兼治。对细颗粒物进行预测和溯源，及时有效地发现细颗粒物变化趋势和污染源分布，为有关部门的科学决策和事后追责提供有力支撑，为公众的及时预防提供准确指导，具有很重要的现实意义。

现阶段在细颗粒物预测方面，利用经典的扩散模型如高斯羽流扩散模型或者社区多尺度空气质量模型，可以得到比较准确的结果。但共同的缺点是很难确定应用的条件，很多关键参数的设置都很繁琐，需要的现实环境建模数据输入较多。对于空气中的气体或者细颗粒物等物质的预测应用较多的是人工神经网络模型(ANN)。然而，面对室外复杂环境下的预测，人工神经网络模型输入等条件的限制可能导致模型产生复杂的网络结构进而导致预测估计时间过长。尽管人工神经网络在非线性建模方面有较强的能力，但面对小样本数据显得有些无力。另一方面，在细颗粒物空间维度分布预测数据需要支撑后续溯源算法的实现，要实现准确的溯源特征，就必须保证数据在突变状态下的敏感性，而人工神经网在处理非平稳数据方面跟踪性误差较大。基于以上考虑选择应用高斯过程回归模型来完成细颗粒物空间和时间维度上的预测，其在处理小样本、非线性等复杂问题时具有很好的效果，相较于人工神经网络具有容易实现、超参数自适应获取、泛化能力强等特点，并且高斯过程回归模型在处理非平稳性问题时，跟踪性更好误差更小。

多模态函数又叫多峰值函数，多模态函数优化问题(Multimodal FunctionOptimization Problems)的一种求解是寻找目标函数全局最优，同时避免局部最优解。另一种优化问题的目标是找到具有相等或不等函数值的多个最优解。

现实环境中雾霾天气的出现，并不是单污染源造成的，其主要来自多污染源多空气污染物复杂作用的结果，通过现有研究可以知道城市雾霾污染的来源主要包括：交通密集区尾气排放、工业密集区废气排放和监测区域以外传入等情况。与源定位相关的大部分工作都致力于寻找、检测和跟踪单个源位置点，当研究多个源位置点时问题就变得复杂化。由于扩散过程由于受到室外环境的多种因素影响，细颗粒物在空间中的分布是动态变化的，这就使得在函数模型创建过程中并不能应用静态的多模态函数建模。而在创建动态气体模型过程中，模型需要设定的环境假设过多，并且对现实环境的空气变化状态有各种要求。对于多模态函数最优值的求解，萤火虫优化算法是一个很好的选择。

综上所述，本发明提出一种基于高斯过程回归和萤火虫算法的细颗粒物预测与溯源方法。该方法首先通过高斯过程回归对细颗粒物空间分布和短时间变化态势进行预测，然后通过萤火虫算法对预测数据进行多模函数的求解，进而完成对细颗粒物的溯源。该方法可以对细颗粒物不稳定变化的准确跟踪，并实现了对细颗粒物多污染源的准确定位。

本发明是针对空气污染监测领域提出的一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测与溯源方法，旨在对细颗粒物空间分布及短时变化态势进行预测，并对细颗粒物多污染源动态回溯，实现对细颗粒物多污染源的准确定位。

假设预测模型训练样本集为：

D＝{(x⁽ⁱ⁾,y⁽ⁱ⁾)|i＝1,2,3,L,n}＝(X,y) (1)

其中，x⁽ⁱ⁾∈R^d为d维输入矢量，表示输入的特征向量，X＝[x⁽¹⁾,x⁽²⁾,L,x⁽ⁿ⁾]表示为d×n维输入矩阵。y⁽ⁱ⁾∈R表示第i个样本对应的输出标量，而y则表示由每个样本对应的输出标量组成的输出矢量。回归问题就要通过训练样本集学习到输入矩阵X与输出矢量y的映射关系：(f(·):R^d→R)，从而得到新的测试点x_*对应的输出f(x_*)。

高斯过程(GP)是具有联合高斯分布的任意个随机变量的集合。高斯过程的性质可以通过其平均函数和协方差函数完全指定。可以表示为：

其中x、x′∈R^d表示输入的特征向量，因此高斯过程可以有以下表示形式：

f(x)～GP(m(x),k(x,x′)) (3)

在实际建模过程中，通过预处理使得m(x)＝0，并且考虑数据噪声的存在，如设备误差等等，并且进一步设定噪声为因此对于回归问题，得到以下模型：

其中y为观测到的带有噪声的数值，于是进一步得到：

对于带噪声的观测训练数据集D＝{(x⁽ⁱ⁾,y⁽ⁱ⁾)|i＝1,2,L,n}，通过高斯过程回归模型(GPR)从训练数据集中学习函数来预测点集该函数涉及假设的函数高斯先验分布：

其中，K(X,X)＝K_n＝(k_ij)为n×n阶的协方差矩阵，并且由高斯过程定义可得，此矩阵是对称正定的，而且协方差矩阵中的元素k_ij＝k(x_i,x_j)表示x_i,x_j所处状态的线性联系的密切程度，若其绝对值较大时，则表示两者的线性关系较密切；若其绝对值较小，则表示两者的线性关系不密切，K(X_*,X)＝K_*，K_*为测试点集输入X_*与训练数据集X的协方差矩阵。最后K(X_*,X_*)＝K_**则是测试点输入X_*的方差，I_n为n阶单位阵。

根据贝叶斯法则计算新点的分布，可以得到高斯过程输出的后验概率分布满足以下公式：

其中，对y_*的最好估值是这个分布的均值：

并且估值的不确定性被记录在其差异中：

到这为止，完成了相应基本模型的构建。

由创建的高斯过程回归模型可以得到，模型的协方差选择可以有多种。一般而言，半正定核是有效的协方差函数。平方指数协方差函数形式如下：

其中，来控制先验方差的，最大允许协方差为当某个函数在y轴的数据范围很大时，其数值相应增大。M＝diag(l²)，其中l∈R^d被称为方差尺度，用来控制协方差衰减率的各向同性长度尺度参数，其中ω_i∈l表示输入的第i个特征的权重，其中集合为超参数。

本发明考虑的特征参数有温度和湿度(时间维度短时间的变化态势预测特征参数包括：时间，温度，湿度；空间维度的分布预测包括：位置坐标，温度，湿度)，由于特征参数对细颗粒物的相关性并不一致，因此需要对核函数中选择不同的权重ω_i来实现对不同特征参量的加权。训练样本能够反应真实的样本分布，并且符合独立同分布，为了能够让模型从以上训练数据中优化超参数，本发明使用极大似然估计进行求解：

从公式(11)中可以看到求解超参数时，要对求逆运算，这样计算的时间复杂度为O(n³)，而在进行数据预测时，时间复杂度则为O(n²)，这种复杂度的情形下，超参数的优化求解结果由于数据的不稳定性容易满足不了模型精度的要求，最后由于迭代次数到而跳出训练。

为了提高细颗粒物的预测精度，本发明应用一种改进的粒子群优化算法对高斯过程回归方法进行改进。基于样本学习优化参数是高斯过程回归算法的一个关键特性。在高斯过程回归模型的超参数求解过程中利用改进的粒子群优化算法提高超参数的准确度。粒子群优化算法是利用群体中粒子间的协作和信息共享产生的群体智能，来实现最后的全局寻优。在粒子群优化算法中，每个粒子代表候选解决方案，根据自己及其粒子同伴的历史经验活动在D维空间S中。解的质量好坏取决于自适应度函数。第i个粒子的速度和位置由以下式子表示：

x_i＝(x_i1,x_i2,x_i3,L,x_id,L,x_iD) (12)

v_i＝(v_i1,v_i2,v_i3,L,v_id,L,v_iD) (13)

粒子群优化算法(PSO)随机初始化一组粒子，并通过随机搜索找到最优解。在迭代过程中，每个粒子根据其迄今为止已经找到的两个“适应度值”来指导其更新速度与位置状态，这两个“适应度值”分别为：单粒子i的最佳位置被存储为：

p_i＝(p_i1,p_i2,L,p_iD) (14)

表示为pbest，并且粒子群中的最佳位置也被记录存储：

p＝(p_g1,p_g2,L,p_gD) (15)

表示为gbest，其中gbest就是pbest中最好的。当获得gbest和pbest后，第i个粒子通过下式来更新自己的状态：

v_i(t+1)＝ωv_i(t)+c₁r₁[pbest_i-x_i(t)]+c₂r₂[gbest_i-x_i(t)] (16)

x_i(t+1)＝x_i(t)+v_i(t+1) (17)

公式中i＝1,2,3,L,N，其中N表示粒子个数；r₁,r₂代表两个从U(0,1)中采样的随机序列；c₁,c₂是两个学习因子，通常情况下c₁＝c₂＝2；ω称为关惯性因子，并且ω≥0，通过公式(16)可以知道，ω越大，其全局寻优能力会很强，局部寻优能力就会较差，同理当其越小时，效果相反。本发明应用线性递减权值策略来获得动态的ω，用以得到更好的寻优效果：

ω(t)＝(ω_ini-ω_end)(G_K-g)/G_k+ω_end (18)

其中G_k为最大迭代次数，ω_ini权值一般设定为0.4，ω_end一般设定为0.9。

为了保持多样性，避免早熟和容易陷入局部最优的问题，引入差分进化(Differential Evolution)算子。使得改进的粒子群优化算法中，PSO保持粒子群的动态性，同时通过差分进化保持种群的多样性。

在粒子群已经获得gbest和pbest后，我们将对gbest进行差分进化变异操作。在[0,1]内产生一个随机的r，在[1,D]内产生一个随机的整数值k。当r小于给定的CR或者维数d等于k时，可以得到第d维的变异点Trail：

Trail_d＝gbest_d+F(pbest_1d+pbest_2d-pbest_3d-pbest_4d)(if r＜CR or d＝k) (19)

其中CR是交叉常数，并且缩放因子F是[0,2]内的常数，pbest_1d，pbest_2d，pbest_3d，pbest_4d是从pbest中随机选择的。

然后可以计算出Trail_d的适应度，只有当它比gbest_d好时，Trail_d才会取代gbest_d。否则，以前的gbest_d将被保留。

接下来把以上改进的粒子群优化算法应用到高斯过程回归模型的超参数优化的求解过程中。基于改进粒子群优化算法的改进高斯过程回归模型在细颗粒物预测的过程主要包括三个步骤：第一步，初始化；第二步，利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化；第三步，预测。在步骤1中构建合适的训练集并选择用于高斯过程回归的协方差函数；在步骤2中，使用基于改进的粒子群优化算法对协方差函数的超参数进行优化，其中把超参数的极大似然概率估计设置为粒子群优化算法的适应度函数为：

在步骤3中，计算预测分布并输出预测。详细步骤如下：

步骤1(初始化)：创建模型的样本训练集并为高斯过程回归模型选择协方差函数。

步骤2(超参数优化)：

步骤2.1：设置PSO参数：首先选定粒子群大小N，设置最大迭代次数iter_max，精度误差为ε，c₁＝c₂＝2，ω_end＝0.9，ω_ini＝0.4，CR＝0.9,F＝0.9，并初始化每个粒子的x_i＝(x_i1,x_i2,x_i3,L,x_id,L,x_iD)和v_i＝(v_i1,v_i2,v_i3,L,v_id,L,v_iD)；

步骤2.2：通过(20)估计每个粒子的适应度，并找出每个粒子的对应的pbest和gbest；

步骤2.3：在[0,1]内生成一个随机数r，在[1,D]生成一随机数k。并根据公式(19)对粒子进行变异操作；

步骤2.4：根据公式(16)，(17)，(18)更新x_i和v_i；

步骤2.5：重复以上过程，达到最大迭代次数或gbest的适应度满足精度误差的要求，该过程将结束。

步骤3(预测)：根据优化的超参数得出模型，并输入x_*得出

由于扩散过程由于受到室外环境的多种因素影响，细颗粒物在空间中的分布是动态变化的，这就使得在函数模型创建过程中并不能应用静态的多模态函数建模。针对以上问题，本发明基于改进的萤火虫算法设计了细颗粒物溯源算法。

萤火虫优化算法(Glowworm swarm optimization)，它是一种同时计算多模态函数最优值的新算法。算法首先将n个萤火虫的种群随机放置在设定好范围的搜索空间中，以使它们很好地分散。最初，所有萤火虫都含有等量的荧光素l ₀。每次迭代包括一个荧光素更新阶段，然后是基于转换规则的移动阶段。整个算法主要分为一下三个阶段：第一阶段，荧光素更新阶段；第二部分，运动阶段实现位置更新；第三阶段，邻域范围更新阶段。

荧光素更新阶段：在荧光素更新阶段期间，每个萤火虫向其之前的荧光素水平增加目标位置在目标函数中的适应性成正比的荧光素量。此外，荧光素值的一部分被减去以模拟荧光素随时间的衰减。荧光素更新规则由下式给出：

l_i(t+1)＝(1-ρ)l_i(t)+γJ(x_i(t+1)) (21)

其中l_i(t)表示与萤火虫i在t时刻的荧光素大小，ρ是荧光素衰减常数(0＜ρ＜1)，γ是荧光素增强常数，J(x_i(t))表示t时i位置处目标函数的值。

运动阶段：在运动阶段，每个萤火虫使用概率机制决定朝着一个荧光素值高于其自身的邻居移动。对于每个萤火虫i，向邻居j移动的概率由下式给出：

其中j∈N_i(t)，表示时间t处的萤火虫i的邻居的集合，d_ij(t)表示时间t处的萤火虫i和j之间的欧几里得距离，并且表示在时间t与萤火虫i相关的邻域范围。让萤火虫根据公式(3-2)给出的p_ij(t)选择一个j∈N_i(t)的萤火虫。那么，萤火虫运动的离散时间模型可以表述为：

其中x_i(t)∈R^m是在时间t时萤火虫i在m维实空间R^m中的位置，||·||表示欧几里得范数算子，s(s＞0)表示是步长。

邻域范围更新阶段：本发明将每个萤火虫与一个邻域相关联，这些邻域的范围是动态的如果在邻域范围设定中选择固定邻域范围，当萤火虫仅依赖于局部信息来决定它们的移动时，预计捕获到的目标函数峰值数量将是径向感知范围的函数。由于本发明假设关于目标函数(例如：峰值数量和峰间距离)的先验信息不可用，很难找到合适的邻域范围。因此，萤火虫算法使用自适应邻域范围来检测多峰函数格局中多峰的存在。

令r₀为每个萤火虫的初始邻域范围(即)。为了自适应地更新每个萤火虫的邻域范围，应用以下规则：

其中β是常数参数，n_t是用于控制邻居数量的参数。

GSO在应用过程中存在一定的缺点，其主要表现为搜索精度低、在后期迭代过程中收敛缓慢、而且峰值发现率低的问题。本发明在萤火虫移动步长和和萤火虫自探索两方面进行优化设计。在移动步长s改进方面，考虑针对每一个萤火虫个体进行自适应步长设计，让每一个萤火虫个体拥有属于自己的步长。每个萤火虫个体的邻域集合N_i(t)中的优秀个体数量n_t是一定的，如果萤火虫个体周围优秀个体数比较密集时，会通过迭代操作，减小邻域半径，使n_t大小满足参数设定值。当萤火虫个体周围优秀萤火虫个体比较稀疏时，算法通过迭代增加的大小，进而增加邻域内个体的个数。并且每次萤火虫活动的目标位置只会在N_i(t)中选择，所以N_i(t)中的个体直接决定着萤火虫个体步长的合适与否。所以通过下式，先计算得到邻域集合中个体与目标萤火虫个体的平均距离：

的范围：在初始化时另当时，令其中q＝0.96。

萤火虫自探索方面，萤火虫个体的移动是完全因为邻域内优秀个体，对优秀个体依赖性过高，当感应范围内没有找到优秀个体时将停止移动，这就这种搜索策略使得萤火虫个体本身的局部搜索能力得不到充分发挥，这也使得整体算法的搜索性能得到严重影响。

当萤火虫个体在其感知范围内找不到优秀个体时，在此范围内没有找到优秀的个体，那就以现在的位置为起点，进行对数螺旋搜索，使得搜索面积变大，并且在移动过程中移动两圈，如果移动后，新位置的荧光素水平比原来低，而且周围没有优秀的个体，萤火虫会以下式跳出该位置：

x_i(t+1)＝x_i(t)+(1-2rand())·s_i(t) (26)

如果更新后荧光素值仍然小于原来位置水平，则取消位置更新。在式(26)中，rand()是(0,1)的随机数，至此完成火虫自探索机制改进。

通过以上两方面对萤火虫算法的改进，本发明实现了萤火虫算法自适应步长和自探索两种机制的优化，提高了萤火虫算法的搜索效率和搜索精度。

本专利提出一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测与溯源方法：

1.本专利提出一种基于高斯过程回归和萤火虫算法的细颗粒物预测与溯源方法，该方法将高斯过程回归和萤火虫算法相结合，既能细颗粒物空间分布及短时变化态势进行预测，也能对细颗粒物多污染源动态回溯，实现对细颗粒物多污染源的准确定位。

2.本方法利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化，提高了细颗粒物的预测精度，并对萤火虫算法在萤火虫移动步长和和萤火虫自探索两方面进行优化设计，提高了萤火虫算法的搜索效率和搜索精度。

图1细颗粒物预测方法架构图

图2基于改进PSO算法的改进GPR算法流程图

图3改进的萤火虫算法实现流程图

图4细颗粒物溯源方法架构图

本发明提出一种基于高斯过程回归和萤火虫算法相结合5的细颗粒物预测与溯源方法，该方法将高斯过程回归和萤火虫算法相结合，既能细颗粒物空间分布及短时变化态势进行预测，也能对细颗粒物多污染源动态回溯，实现对细颗粒物多污染源的准确定位。下面结合附图进行详细说明。

图1展示了细颗粒物预测方法架构图，其中在单节点时间维度预测中，本发明在线下利用MongoDB数据库中历史统计数据样本完成对模型的训练，并通过训练后模型预测未来系统设定短时间内某些时间节点的数值，通过把每个节点的训练数值输入到空间维度分布情况预测模型中，获取到未来时间细颗粒物空间分布情况，用于渲染百度地图上层热力图。此外，在空间维度分布情况预测中，利用实时数据来在线完成模型训练，并得到空间分布预测值，用于渲染热力图，其中，在线训练模型需要利用每一帧数据训练一个新的模型。图2是基于改进的粒子群算法的改进高斯过程回归算法流程图。图3是改进的萤火虫算法实现流程图。图4是细颗粒物溯源方法架构图，细颗粒物溯源方法实施如下：

实际部署后整个网络监测数据更新周期为1分钟，系统每1个周期得到的监测数据作为细颗粒物溯源模块的1帧数据，对于每1帧数据会通过细颗粒物预测模块对其空间维度分布情况进行预测，获取到整个监测环境下细颗粒物在此时刻的分布情况；

每1帧静态的预测数据就是一个多模态函数，本发明应用改进的萤火虫算法进行多模态函数的求解，获取到此帧数据的局部、全局最优解，并标记为本帧数据的目标源，在得到每1帧数据的目标源后，会记录本帧数据目标源在百度地图上的经纬度坐标信息；

设定细颗粒物的污染溯源阈值，当细颗粒物浓度大于这一阈值时，开始对每1帧数据进行溯源分析与源点标记，得到每1帧数据的目标源快照，并把模型的输入输出数据读入MySQL数据库，其中包括每1帧空间维度预测数据和百度地图标记点坐标；

当监测区域出现细颗粒物超标时系统会通过动态刷新的形式，用每1帧目标源快照数据渲染百度地图上层热力图，并动态刷新出来。同时设定相同目标源距离阈值，通过阈值找到了每1帧快照的同一个污染源迁移点，在前端动态刷新的过程中，把同一个迁移点用线连接起来，直观地展示确切的污染源点位置。

一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测与溯源方法，通过高斯过程回归对细颗粒物空间分布和短时间变化态势进行预测，然后通过萤火虫算法对预测数据进行多模函数的求解，进而完成对细颗粒物的溯源。

利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化，其步骤如下：

步骤1.创建模型的样本训练集并为高斯过程回归模型选择协方差函数；

步骤2.利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化；

步骤3.根据优化的超参数得出模型，输出预测。

对萤火虫算法在移动步长方面改进，针对每一个萤火虫个体进行自适应步长设计，让每一个萤火虫个体拥有属于自己的步长。

对萤火虫算法在自探索方面改进，当萤火虫个体在其感知范围内找不到优秀个体时，在此范围内没有找到优秀的个体，那就以现在的位置为起点，进行对数螺旋搜索，使得搜索面积变大，并且在移动过程中移动两圈，如果移动后，新位置的荧光素水平比原来低，而且周围没有优秀的个体，萤火虫会以下式跳出该位置：

x_i(t+1)＝x_i(t)+(1-2rand())·s_i(t)

其中x_i(t)表示时间t萤火虫i所处位置，s_i(t)表示萤火虫步长，rand()是(0,1)的随机数；如果更新后荧光素值仍然小于原来位置水平，则取消位置更新。

Claims (4)

1.一种基于高斯过程回归和萤火虫算法相结合的细颗粒物预测与溯源方法，其特征在于，包括以下步骤：

(1)对萤火虫算法进行优化，得到萤火虫优化算法；

(2)创建模型的样本训练集并为高斯过程回归模型选择协方差函数；

(3)利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化；

(4)根据优化的超参数得出模型，输出预测。

2.根据权利要求1所述的方法，其特征在于，所述对萤火虫算法进行优化，得到萤火虫优化算法，包括：

(1.1)荧光素更新阶段

在荧光素更新阶段期间，每个萤火虫向其之前的荧光素水平增加目标位置在目标函数中的适应性成正比的荧光素量，此外，荧光素值的一部分被减去以模拟荧光素随时间的衰减，荧光素更新规则由下式给出：

l_i(t+1)＝(1-ρ)l_i(t)+γJ(x_i(t+1))

其中，l_i(t)表示与萤火虫i在t时刻的荧光素大小，ρ是荧光素衰减常数(0＜ρ＜1)，γ是荧光素增强常数，J(x_i(t))表示t时i位置处目标函数的值；

(1.2)运动阶段

在运动阶段，每个萤火虫使用概率机制决定朝着一个荧光素值高于其自身的邻居移动，对于每个萤火虫i，向邻居j移动的概率由下式给出：

萤火虫根据上式给出的p_ij(t)选择一个j∈N_i(t)的萤火虫；

其中，j∈N_i(t)，表示时间t处的萤火虫i的邻居的集合，d_ij(t)表示时间t处的萤火虫i和j之间的欧几里得距离，并且表示在时间t与萤火虫i相关的邻域范围

萤火虫运动的离散时间模型为：

其中，x_i(t)∈R^m是在时间t时萤火虫i在m维实空间R^m中的位置，||.||表示欧几里得范数算子，s(s＞0)表示是步长；

(1.3)邻域范围更新阶段

令r₀为每个萤火虫的初始邻域范围，为了自适应地更新每个萤火虫的邻域范围，应用以下规则：

其中，β是常数参数，n_t是用于控制邻居数量的参数；

通过下式，计算得到邻域集合中个体与目标萤火虫个体的平均距离：

的范围：在初始化时令当时，令q＝0.96；

(1.4)萤火虫自探索方面

当萤火虫个体在其感知范围内找不到优秀个体时，在此范围内没有找到优秀的个体，那就以现在的位置为起点，进行对数螺旋搜索，使得搜索面积变大，并且在移动过程中移动两圈，如果移动后，新位置的荧光素水平比原来低，而且周围没有优秀的个体，萤火虫会以下式跳出该位置：

x_i(t+1)＝x_i(t)+(1-2rand())·s_i(t)

如果更新后荧光素值仍然小于原来位置水平，则取消位置更新

其中，rand()是(0,1)的随机数。

3.根据权利要求1所述的方法，其特征在于，所述创建模型的样本训练集并为高斯过程回归模型选择协方差函数，包括：

模型训练样本集为：

D＝{(x⁽ⁱ⁾,y⁽ⁱ⁾)|i＝1,2,3,L,n}＝(X,y)

其中，x⁽ⁱ⁾∈R^d为d维输入矢量，表示输入的特征向量，X＝[x⁽¹⁾,x⁽²⁾,L,x⁽ⁿ⁾]表示为d×n维输入矩阵，y⁽ⁱ⁾∈R表示第i个样本对应的输出标量，而y则表示由每个样本对应的输出标量组成的输出矢量，回归问题就要通过训练样本集学习到输入矩阵X与输出矢量y的映射关系：(f(·):R^d→R)，从而得到新的测试点x_*对应的输出f(x_*)

平方指数协方差函数形式如下：

其中，来控制先验方差的，最大允许协方差为当某个函数在y轴的数据范围很大时，其数值相应增大，M＝diag(l²)，其中l∈Rd被称为方差尺度，用来控制协方差衰减率的各向同性长度尺度参数，其中ω_i∈l表示输入的第i个特征的权重，其中集合为超参数。

4.根据权利要求1所述的方法，其特征在于，所述利用改进的粒子群优化算法对高斯过程回归模型的方程超参数优化，包括：

(4.1)设置PSO参数：首先选定粒子群大小N，设置最大迭代次数iter_max，精度误差为ε，c₁＝c₂＝2，ω_end＝0.9，ω_ini＝0.4，CR＝0.9,F＝0.9，并初始化每个粒子的x_i＝(x_i1,x_i2,x_i3,L,x_id,L,x_iD)和v_i＝(v_i1,v_i2,v_i3,L,v_id,L,v_iD)；

(4.2)通过下式估计每个粒子的适应度，并找出每个粒子的对应的pbest和gbest；

(4.3)在[0,1]内生成一个随机数r，在[1,D]生成一随机数k，并根据下式对粒子进行变异操作；

Trail_d＝gbest_d+F(pbest_1d+pbest_2d-pbest_3d-pbest_4d)(if r＜CR or d＝k)

其中，CR是交叉常数，缩放因子F是[0,2]内的常数，pbest_1d、pbest_2d、pbest_3d、pbest_4d是从pbest中随机选择的；

(4.4)根据下式更新x_i和v_i；

v_i(t+1)＝ωv_i(t)+c₁r₁[pbest_i-x_i(t)]+c₂r₂[gbest_i-x_i(t)]

x_i(t+1)＝x_i(t)+v_i(t+1)

其中，公式中i＝1,2,3,L,N，N表示粒子个数；r₁,r₂代表两个从U(0,1)中采样的随机序列；c₁,c₂是两个学习因子，c₁＝c₂＝2；ω称为关惯性因子，并且ω≥0

ω(t)＝(ω_ini-ω_end)(G_K-g)/G_k+ω_end

其中G_k为最大迭代次数，ω_ini权值为0.4，ω_end设定为0.9；

(4.5)重复以上过程，达到最大迭代次数或gbest的适应度满足精度误差的要求，该过程将结束。