CN105784556A - 一种基于自组织模糊神经网络的空气细颗粒物pm2.5软测量方法 - Google Patents

Abstract

一种基于自组织模糊神经网络的空气细颗粒物PM_2.5软测量方法既属于环境工程领域又属于检测技术领域。大气环境系统具有多变量、非线性、内部机理复杂、信息不完备等特性，难以通过机理分析建立其数学模型，而神经网络对于高度非线性和严重不确定性系统具有较好的处理能力。本发明针对PM_2.5难以预测的问题，基于偏最小二乘选取与PM_2.5相关的辅助变量，然后利用基于敏感度分析方法的自组织模糊神经网络建立相关变量和PM_2.5之间的软测量模型，对PM_2.5进行预测，取得了较好的效果，为环境管理决策者和群众提供及时准确的大气环境质量信息，有利于加强大气环境污染控制，防止严重污染的发生。

Description

一种基于自组织模糊神经网络的空气细颗粒物PM2.5软测量方法

技术领域

本发明涉及空气中细颗粒物PM_2.5的软测量方法。软测量是根据某种最优准则，选择一组既与主变量密切联系，又容易测量的辅助变量，通过构造某种数学模型，依靠事先学习和记忆实现对主变量的估计。PM_2.5的预测在空气污染防治中具有重要意义，将软测量方法应用于PM_2.5的预测中，既可节约空气污染监测成本，又能及时提供PM_2.5的浓度信息，有利于加强大气环境污染控制，是先进制造技术领域的重要分支，既属于环境工程领域，又属于检测技术领域。

背景技术

国务院印发的《大气污染防治行动计划》明确指出，到2017年全国地级及以上城市可吸入颗粒物浓度比2012年下降10％以上，其中北京市细颗粒物浓度控制在60μg/m³。但是，当前我国大气污染形势严峻，以PM_2.5等为特征污染物的区域性大气环境问题日益突出，损害公众身体健康，影响社会和经济的可持续发展。因此通过对PM_2.5的预测，及时防治空气污染，防止严重污染的发生已受到社会的广泛关注。然而，PM_2.5的浓度不仅和污染源、污染物有关，还和气象条件有关，其成分及比例，排放源组成及贡献率有明显的时空异质性，且大气环境系统具有多变量、非线性、内部机理复杂、信息不完备等特性，使PM_2.5的预测难度较大。

目前，各国环保部门广泛采用的PM_2.5测定方法有重量法、β射线吸收法和微量振荡天平法。重量法需要人工称重，操作繁琐费时。后两者属于自动监测方法，所需仪器仪表价格昂贵，维护困难，测量范围有限。常用的PM_2.5预测方法确定性化学建模方法所需的模型分辨率、气象初始条件、排放源的时空分布等参数难以确定，计算复杂，精度较低。线性回归建模方法不适合对本身非线性的大气环境系统进行建模，人工神经网络可解释性差，一般模糊神经网络结构固定。因此，研究新的PM_2.5预测方法解决PM_2.5难以预测的问题具有深远的现实意义。

本发明提出一种新的基于自组织模糊神经网络(SOFNN)的空气细颗粒物PM_2.5软测量方法。该自组织模糊神经网络利用频域上的模型输出敏感度分析方法确定模糊神经网络规则化层神经元输出对网络输出的贡献率，根据贡献率大小增删规则化层神经元，使模糊神经网络的结构得到自动调整，同时利用梯度下降算法训练模糊神经网络的中心、宽度和权值。然后，以通过偏最小二乘选取的辅助变量为输入，利用该自组织模糊神经网络建立PM_2.5软测量模型，实现对空气细颗粒物PM_2.5的预测。

发明内容

本发明获得了一种基于自组织模糊神经网络的空气细颗粒物PM_2.5软测量方法。该方法首先利用空气污染物浓度和气象实测数据基于偏最小二乘选取与PM_2.5相关的辅助变量，然后利用基于敏感度分析的自组织模糊神经网络建立相关变量和PM_2.5之间的软测量模型，对PM_2.5进行预测，解决了PM_2.5难以预测的问题。

一种基于自组织模糊神经网络的空气细颗粒物PM_2.5软测量方法设计包括以下步骤：

(1)确定辅助变量和预测变量。

①变量数据的标准化处理：将空气污染物浓度及气象变量的上一小时数据和PM_2.5浓度的下一小时数据一一对应，整理出L组数据，其中，L取值可为200左右，过小数据量不够，包含信息不完备，过大计算过于复杂，然后以空气污染物浓度和气象变量为自变量X，PM_2.5浓度为因变量Y：

X＝[x₁,x₂,…,x₇]^T,Y＝y(1)

其中，x₁,x₂,…,x₇分别表示温度、相对湿度、风速、CO、NO₂、O₃、SO₂浓度，y表示PM_2.5浓度，温度的单位是摄氏度，风速的单位是级，各空气污染物浓度的单位是μg/m³，相对湿度用百分数表示，没有单位。数据标准化处理的计算公式如下：

$g_1=\frac{D_1-D_{1,\min }}{D_{1,max}-D_{1,min}}---\left(2\right)$

其中，g₁是温度的上一小时实测数据数组标准化后得到的数组，D₁是标准化前温度的实测数据数组，D_1,min是数组D₁中的最小值，D_1,max是数组D₁中的最大值。相对湿度、风速、CO、NO₂、O₃、SO₂浓度的上一小时实测数据数组和PM_2.5浓度的下一小时实测数据数组分别是D₂，D₃，D₄，D₅，D₆，D₇和D₈，找出各数组内部的最大值和最小值，按照式(2)相同的处理方式得到标准化数组分别是：g₂，g₃，g₄，g₅，g₆，g₇和g₈。g₁，g₂，g₃，g₄，g₅，g₆和g₇构成了自变量X的标准化观测数据阵E₀，g₈是因变量Y的标准化观测数据阵F₀。

②计算的最大特征值对应的单位特征向量W₁，求取自变量组的第一个成分其得分向量d₁＝E₀W₁，残差矩阵 其中α₁和β₁是模型效应负荷量，且

③计算的最大特征值对应的单位特征向量W₂，求取自变量组的第二个成分其得分向量d₂＝E₁W₂，残差矩阵 其中α₂和β₂是模型效应负荷量，且将和代入Y＝t₁β₁+t₂β₂中得到提取2个成分时因变量Y关于自变量X的偏最小二乘回归方程式：

y＝a₁x₁+…+a₇x₇(3)

其中，a₁,…a₇分别是提取2个成分时自变量x₁,…,x₇对应的回归系数。

④利用交叉有效性检验提取2个成分时模型精度是否满足要求：首先每次从E₀和F₀中舍去第q组数据(q＝1,2,…,L)，利用剩下的L-1组数据建立提取2个成分时的偏最小二乘回归式，并求出在舍去的数据处的预测值y_(-q)(2)，真实值是y_q。对q＝1,2,…,L重复以上操作，即得提取2个成分时因变量Y的预测误差平方和PRESS(2)：

$PRESS\left(2\right)={\Σ}_{q=1}^L{(y_q-y_{(-q)}(2\left)\right)}^2---\left(4\right)$

然后利用E₀和F₀拟合含1个成分的偏最小二乘回归方程式，记第q组数据的预测值为y_(q)(1)，则提取1个成分时Y的误差平方和SS(1)为：

$SS\left(1\right)={\Σ}_{q=1}^L{(y_q-y_{\left(q\right)}(1\left)\right)}^2---\left(5\right)$

定义提取2个成分时的交叉有效性为Q²(2)＝1-PRESS(2)/SS(1)，若Q²(2)＜1-0.95²，则模型精度满足要求，可停止提取成分，否则继续提取成分。

⑤通过上述步骤得到因变量Y关于自变量X的偏最小二乘回归方程式后，根据回归方程式中回归系数的绝对值，除去绝对值最小的自变量，将剩下的自变量选作辅助变量，记为r₁,r₂,…,r₆，PM_2.5浓度作为预测变量，是自组织模糊神经网络的期望输出，记为o_d，辅助变量和预测变量的准化观测数据阵分别记为E和F，选取E和F的前I组数据作为训练样本，后I′组数据作为测试样本(I′＝L-I，I′<＝I)。

(2)设计用于PM_2.5预测的自组织模糊神经网络的初始拓扑结构。本发明用于PM_2.5预测的自组织模糊神经网络共四层：输入层、RBF层、规则化层和输出层。输入是选取的辅助变量，输出是PM_2.5浓度，确定其6-P-P-1的初始连接方式，即输入层神经元数和辅助变量的个数同是6，RBF层神经元数是P，规则化层神经元数是P，其中P是正整数，在[1-10]之间取值，且避免P值过大使自组织模糊神经网络的计算负荷过大，输出层神经元数是1。自组织模糊神经网络的中心、宽度和权值的初始值设定在(0,1)之间，随机设定即可，不影响自组织模糊神经网络的预测效果，只影响其收敛速度。自组织模糊神经网络的输入表示为r＝[r₁,r₂,…,r₆]^T，[r₁,r₂,…,r₆]^T是[r₁,r₂,…,r₆]的转置，期望输出是o_d。第k组训练样本表示为r(k)＝[r₁(k),r₂(k),…,r₆(k)]^T，k＝1,2,…,I，用第k组训练样本训练自组织模糊神经网络时，自组织模糊神经网络的各层输出依次为：

①输入层：该层有6个神经元：

u_i(k)＝r_i(k)(6)

其中，u_i(k)是输入层第i个神经元的输出，该层的输入向量为r(k)＝[r₁(k),r₂(k),…,r₆(k)]^T，r₁(k),r₂(k),…,r₆(k)分别表示对应的辅助变量。

②RBF层：该层选取RBF神经元的高斯函数作为隶属函数对输入变量进行模糊化处理，RBF层有P个神经元，该层第j个神经元的输出为：

其中，c_ij(k)，σ_ij(k)分别是模糊神经网络的中心和宽度。

③规则化层：该层神经元数目与RBF层相同，该层第l个神经元的输出v_l(k)为：

其中，是RBF层第l个神经元的输出。

④输出层：该层有1个神经元，该层的输出表示PM_2.5浓度，如下式所示：

$o\left(k\right)={\&Sigma;}_{l=1}^P{w}_l\left(k\right)v_l\left(k\right)---\left(9\right)$

其中，w_l(k)是规则化层第l个神经元与输出层神经元之间的连接权值。式(6)-式(9)中，i＝1,2,…,6，j＝1,2,…,P，l＝1,2,…,P。自组织模糊神经网络的训练均方根误差(RMSE)为：

$RMSE={\mathrm{\&Sigma;}}_{k=1}^I\frac{{\left(o_d\right(k)-o(k\left)\right)}^2}{I}---\left(10\right)$

其中，o_d(k)和o(k)分别是第k组训练样本的期望输出和用第k组训练样本训练模糊神经网络时的网络输出，训练模糊神经网络的目的是使式(10)定义的训练RMSE达到期望值。

(3)用训练样本训练自组织模糊神经网络。在训练过程中，利用频域上的模型输出敏感度分析方法确定模糊神经网络规则化层神经元输出对网络输出的贡献率，根据贡献率大小增删规则化层神经元，以此来自动调整模糊神经网络的结构，并利用梯度下降算法训练模糊神经网络的中心、宽度和权值，直到网络精度满足信息处理的需求，具体为：

①利用训练样本和梯度下降算法训练给定初始模糊神经网络，网络每产生一组输出参数更新一次，对I组训练样本进行反复训练直到满足结构修改步数(取值为20-30步之间，过小对模糊神经网络规则化层神经元输出信息采集不够，过大会降低网络结构调整效率)后转向②，设置最大迭代步数是1000，期望训练RMSE设定为0.01。梯度下降算法参数调整公式如下：

$\left\{\begin{array}{c}{c}_{ij}(k+1)=c_{ij}\left(k\right)-\mathrm{\&eta;}\&part;e\left(k\right)/\&part;c_{ij}\left(k\right)\\ {\mathrm{\&sigma;}}_{ij}(k+1)={\mathrm{\&sigma;}}_{ij}\left(k\right)-\mathrm{\&eta;}\&part;e\left(k\right)/\&part;{\mathrm{\&sigma;}}_{ij}\left(k\right)\\ w_l(k+1)=w_l\left(k\right)-\mathrm{\&eta;}\&part;e\left(k\right)/\&part;w_l\left(k\right)\end{array}\right.---\left(11\right)$

其中，c_ij(k)、σ_ij(k)和w_l(k)分别是第k组训练样本训练自组织模糊神经网络时网络的中心、宽度和权值，c_ij(k+1)、σ_ij(k+1)和w_l(k+1)分别是第k+1组训练样本训练自组织模糊神经网络时网络的中心、宽度和权值，目标函数e(k)＝(o_d(k)-o(k))²/2，o_d(k)和o(k)分别是第k组训练样本的期望输出和用第k组训练样本训练自组织模糊神经网络时的网络输出，η是学习率，文中取值是0.005，η过大会使自组织模糊神经网络的误差振荡，过小会降低自组织模糊神经网络的收敛速度。

②在频域中计算规则化层第h个神经元的输出v_h对网络输出o的贡献率ST_h(h＝1,…,P)，如下式所示：

${\mathrm{ST}}_h=\frac{S_h}{{\&Sigma;}_{j=1}^P{S}_j}---\left(12\right)$

$S_h=\frac{A_{{\mathrm{\&omega;}}_h}^2+B_{{\mathrm{\&omega;}}_h}^2}{(A_{{\mathrm{\&omega;}}_h}^2+B_{{\mathrm{\&omega;}}_h}^2)+{\&Sigma;}_{\mathrm{\&omega;}=1}^{\max \left({\mathrm{\&omega;}}_{~h}\right)}(A_{\mathrm{\&omega;}}^2+B_{\mathrm{\&omega;}}^2)}---\left(13\right)$

其中，S_h是v_h对o的总效应，A_ω和B_ω及和分别是o的傅里叶展开式f(s)在频率ω和ω_h处的傅里叶系数，且 ω_h是v_h的基频，max(ω_～h)是除v_h之外的其他输入参数基频的最大值，取ω_h＝2Mmax(ω_～h)，M是干扰因子，文中取值是4，f(s)和v_h分别用式(14)和式(15)表示：

$f\left(s\right)={\&Sigma;}_{\mathrm{\&omega;}=-\mathrm{\&infin;}}^{\mathrm{\&infin;}}\left(A_{\mathrm{\&omega;}}cos\right(\mathrm{\&omega;}s)+B_{\mathrm{\&omega;}}sin(\mathrm{\&omega;}s\left)\right)---\left(14\right)$

$v_h\left(s\right)=\frac{p_h+b_h}{2}+\frac{p_h-b_h}{\mathrm{\&pi;}}arcsin\left(sin\right({\mathrm{\&omega;}}_{h}s)---\left(15\right)$

其中，p_h和b_h分别是v_h的最大值和最小值。

③如果ST_h>＝λ₁(文中λ₁＝0.3)，则对规则化层第h个神经元进行分裂。为了降低网络结构调整对网络误差的影响，采用下式设定分裂得到的新的神经元的初始参数：

$\left\{\begin{array}{c}{c}_{\&CenterDot;new1}=c_{\&CenterDot;new2}=c_{\&CenterDot;h}\left(k\right),{\mathrm{\&sigma;}}_{\&CenterDot;new1}={\mathrm{\&sigma;}}_{\&CenterDot;new2}={\mathrm{\&sigma;}}_{\&CenterDot;h}\left(k\right)\\ w_{new1}={\mathrm{\&mu;w}}_h\left(k\right),w_{new2}=(1-\mathrm{\&mu;})w_h\left(k\right)\end{array}\right.---\left(16\right)$

其中，new1和new2是两个新的神经元，c._new1、σ._new1和w_new1分别是神经元new1的中心向量、宽度向量和其与输出层神经元间的连接权值。c._new2、σ._new2和w_new2分别是神经元new2的中心向量、宽度向量和权值。c._h(k)、σ._h(k)和w_h(k)分别是网络结构调整前神经元h的中心向量、宽度向量和权值，μ服从标准正态分布。

④如果ST_h<λ₂(文中λ₂＝0.05)，删除该神经元，同时对网络参数进行调整，如下式所示：

$\{\begin{array}{c}{c}_{\&CenterDot;h-m}=c_{\&CenterDot;h-m}\left(k\right),{\mathrm{\&sigma;}}_{\&CenterDot;h-m}={\mathrm{\&sigma;}}_{\&CenterDot;h-m}\left(k\right)\\ w_{h-m}=w_{h-m}\left(k\right)+w_h\left(k\right)v_h\left(k\right)/v_{h-m}\left(k\right)\end{array}---\left(17\right)$

其中，神经元h-m是与神经元h欧氏距离最小的规则化层神经元，且ST_h-m>＝λ₂，c._h-m、σ._h-m和w_h-m分别是网络修剪后神经元h-m的中心向量、宽度向量和权值，c._h-m(k)、σ._h-m(k)和w_h-m(k)分别是网络修剪前神经元h-m的中心向量、宽度向量和权值，w_h(k)是网络修剪前神经元h的权值，v_h(k)和v_h-m(k)分别是网络修剪前神经元h和神经元h-m的输出。

⑤利用梯度下降算法更新网络参数，若在某步训练过程中训练RMSE<＝0.01或算法迭代了1000步时则停止计算，否则达到结构修改步数后转向②。

(4)对测试样本进行检测。利用测试样本对训练好的自组织模糊神经网络进行测试，自组织模糊神经网络的输出即为PM_2.5的预测结果。

附图说明

图1是本发明的偏最小二乘回归方程式回归系数绝对值的直方图。

图2是本发明的PM_2.5软测量模型拓扑结构图。

图3是本发明的PM_2.5软测量模型针对标准化数据的训练RMSE变化图。

图4是本发明的PM_2.5软测量模型训练过程结构变化图。

图5是本发明的PM_2.5软测量模型训练结果敏感度图。

图6是本发明的PM_2.5软测量模型训练结果图。

图7是本发明的PM_2.5软测量模型训练结果误差图。

图8是本发明的PM_2.5软测量模型测试结果图。

图9是本发明的PM_2.5软测量模型测试结果误差图。

具体实施方式

本发明获得了一种基于自组织模糊神经网络的空气细颗粒物PM_2.5软测量方法。该软测量方法以PM_2.5为输出，以通过偏最小二乘选取的辅助变量为输入，利用基于敏感度分析的自组织模糊神经网络建立PM_2.5软测量模型，对PM_2.5进行预测。

实验数据来自于石家庄世纪公园2014年10月1日至10日空气污染物浓度和气象小时数据，剔除异常及缺失数据后将温度、相对湿度、风速、CO、NO₂、O₃、SO₂浓度的上一小时数据和PM_2.5浓度的下一小时数据一一对应，共整理出150组数据。

一种基于自组织模糊神经网络的空气细颗粒物PM_2.5软测量方法设计包括以下步骤：

(1)确定预测变量和辅助变量。文中利用空气污染物浓度和气象变量实测数据，基于偏最小二乘选取与PM_2.5相关的辅助变量，具体如下：

①变量数据的标准化处理：以温度、相对湿度、风速、CO、NO₂、O₃、SO₂浓度为自变量X，PM_2.5浓度为因变量Y。温度、相对湿度、风速、CO、NO₂、O₃、SO₂浓度的上一小时实测数据数组和PM_2.5浓度的下一小时实测数据数组分别表示为D₁，D₂，D₃，D₄，D₅，D₆，D₇和D₈，找出各数组内部的最大值和最小值，按照式(2)的处理方式得到标准化后的数组分别是：g₁，g₂，g₃，g₄，g₅，g₆，g₇和g₈。g₁，g₂，g₃，g₄，g₅，g₆和g₇构成了自变量X的标准化观测数据阵E₀，g₈是因变量Y的标准化观测数据阵F₀。

②偏最小二乘选取辅助变量：利用E₀和F₀基于偏最小二乘求得因变量Y关于自变量X的回归方程式如下，图1给出了X对Y的回归系数绝对值的直方图。

PM_2.5＝0.2018×T+0.1365×RH-0.0829×WS+0.4009

×CO+0.3336×NO₂-0.0463×O₃+0.1005×SO₂(18)

其中，T是温度，RH是相对湿度，WS是风速。根据式(18)中回归系数的绝对值，除去绝对值最小即和PM_2.5浓度相关性最小的自变量：O₃浓度，最终选取温度、相对湿度、风速、CO、NO₂、SO₂浓度作为辅助变量，表示为r＝[r₁,r₂,…,r₆]^T，PM_2.5浓度作为预测变量，表示自组织模糊神经网络的期望输出，记为o_d。g₁，g₂，g₃，g₄，g₅和g₇构成了辅助变量的标准化观测数据阵E，g₈是预测变量的标准化观测数据阵F。取E和F的前100组数据用作训练样本，后50组数据用作测试样本。

(2)初始化用于PM_2.5预测的自组织模糊神经网络。本发明用于PM_2.5预测的自组织模糊神经网络共四层：输入层、RBF层、规则化层和输出层，其拓扑结构如图2所示。根据辅助变量的选取结果，以温度、相对湿度、风速、CO、NO₂、SO₂浓度为输入，PM_2.5浓度为输出，确定自组织模糊神经网络的输入层神经元数是6，输出层神经元数是1，RBF层和规则化层初始神经元数P设定为1，即自组织模糊神经网络的初始连接方式是6-1-1-1。自组织模糊神经网络的中心、宽度和权值赋初值为(0,1)上的伪随机数，用第k组训练样本训练自组织模糊神经网络时，自组织模糊神经网络的各层输出按照式(6)-式(9)计算，式(10)定义的训练RMSE的期望值设定为0.01。

(3)用训练样本训练自组织模糊神经网络。在训练过程中，利用频域上的模型输出敏感度分析方法确定模糊神经网络规则化层神经元输出对网络输出的贡献率，分裂贡献率大的规则化层神经元，删除贡献率小的规则化层神经元，以此来自动调整模糊神经网络的结构，并利用梯度下降算法训练模糊神经网络的中心、宽度和权值，直到网络精度满足信息处理的需求。具体为：

①利用训练样本和梯度下降算法训练给定初始模糊神经网络，模糊神经网络的中心、宽度和权值按照式(11)进行更新，网络每产生一组输出参数更新一次，对100组训练样本进行反复训练直到满足结构修改步数(文中设定为30步)后转向②，设置最大迭代步数为1000。

②找出每一个规则化层神经元输出v_h的最大值p_h和最小值b_h，取干扰因子M＝4，利用式(12)计算v_h对网络输出o的贡献率ST_h。

③如果ST_h>＝λ₁＝0.3，则对规则化层第h个神经元进行分裂。为了降低网络结构调整对网络误差的影响，采用式(16)设定分裂得到的两个新的神经元的中心、宽度和权值。

④如果ST_h<λ₂＝0.05，删除该神经元，同时利用式(17)调整网络参数。

⑤利用梯度下降算法更新网络参数，若在某步训练过程中训练RMSE<＝0.01或算法迭代了1000步时则停止计算，否则达到结构修改步数后转向②。

(4)对测试样本进行检测。将测试样本作为训练好的自组织模糊神经网络的输入，自组织模糊神经网络的输出即为PM_2.5的预测结果。

图3是PM_2.5软测量模型训练RMSE变化图。图4是PM_2.5软测量模型训练过程结构变化图。图5是PM_2.5软测量模型训练结果敏感度图，图5中每个规则化层神经元的敏感度都在0.05-0.3之间，说明利用敏感度分析调整模糊神经网络结构的有效性。图6是PM_2.5软测量模型训练结果图，X轴：训练样本点，Y轴：PM_2.5浓度(μg/m³)，黑色实线是PM_2.5训练期望值，虚线是自组织模糊神经网络训练输出值。图7是PM_2.5软测量模型训练结果误差图。图8是PM_2.5软测量模型测试结果图，X轴：测试样本点，Y轴：PM_2.5浓度(μg/m³)，黑色实线是PM_2.5测试期望值，虚线是自组织模糊神经网络测试输出值。图9是PM_2.5软测量模型测试结果误差图。

表1-17是本发明的实验数据，表1是自变量O₃的实测值，表2-7是训练样本，表8是PM_2.5训练期望值，表9是自组织模糊神经网络训练输出值，表10-15是测试样本，表16是PM_2.5测试期望值，表17是自组织模糊神经网络测试输出值。

表1.自变量O₃的实测值(μg/m³)：

11	32	0	55	3	14	7	66	32	1
										2	40	0	32	0	3	28	63	31	2
0	47	0	29	0	1	62	72	36	10
										1	47	0	27	0	0	92	72	29	30
2	45	0	6	0	1	77	65	24	51
										1	43	0	4	0	0	37	54	56	56
0	41	0	3	0	0	44	24	47	63
										0	36	1	0	0	0	41	35	31	68
0	14	4	0	1	1	12	38	15	54
										0	8	9	0	4	0	1	35	5	44
0	2	23	0	10	0	3	39	6	27
										0	0	33	0	47	1	5	34	2	22
1	0	46	0	39	1	6	42	2	7
										10	0	54	4	26	2	24	45	1	4
24	0	58	8	24	2	59	45	1	2

训练样本：

表2辅助变量温度的训练输入值(摄氏度)

10	11	18	17	8	22	16	22	16	25
										9	13	18	16	8	19	16	23	15	27
10	15	17	16	8	18	16	22	19	28
										9	18	16	16	12	17	16	20	19	27
9	19	13	17	14	17	16	16	18	25
										10	20	12	11	21	16	17	14	19	23
7	21	16	11	20	16	18	14	18	22
										7	21	17	12	19	16	19	14	19	19
7	21	18	12	23	16	21	13	21	17
										7	19	17	11	22	16	22	12	23	19

表3.辅助变量相对湿度的训练输入值(％)

表4.辅助变量风速的训练输入值(级)

0	1	2	2	0	1	1	1	1	2
										0	1	2	2	0	1	1	1	1	2
1	1	1	2	0	0	1	1	1	1
										1	1	1	2	0	1	1	1	2	1
1	1	1	2	1	1	1	0	2	1
										1	2	1	0	2	1	1	0	2	1
0	1	1	1	2	1	1	1	2	1
										1	2	1	1	1	1	1	1	2	1
1	1	2	1	1	1	1	1	2	1
										1	2	2	0	1	1	1	1	2	1

表5.辅助变量CO的训练输入值(μg/m³)

0.3740	0.5050	0.4540	1.0740	0.7850	1.2330	0.7440	0.4890	1.4840	1.7340
										0.4050	0.5620	0.4650	1.0930	0.6410	1.1900	0.9900	1.0490	1.6180	1.1800
0.4920	0.6120	0.4980	1.1870	0.5800	1.1660	0.9200	0.7220	1.5070	0.8030
										0.4750	0.8330	0.5920	1.1600	0.6000	1.1750	0.6120	0.6500	1.4680	0.5270
0.5290	0.4960	0.6310	1.1670	0.6150	1.1640	0.5560	0.7030	1.4510	0.6360
										0.4350	0.2340	0.7530	1.3500	0.7500	1.1710	0.6060	0.9350	1.4330	0.8790
0.4670	0.2460	0.8490	1.4980	0.9700	1.1130	0.5360	1.0990	1.5270	1.0370
										0.5780	0.3070	0.9770	1.5830	1.0750	0.9640	0.6040	1.1680	1.5870	0.8780
0.6280	0.3760	1.0790	1.2240	1.1020	0.6350	0.5190	1.2570	1.9840	0.9900
										0.5080	0.4240	1.1090	1.0300	1.1860	0.6380	0.4540	1.3270	2.1900	1.2780

表6.辅助变量NO₂的训练输入值(μg/m³)

表7.辅助变量SO₂的训练输入值(μg/m³)

6	8	18	5	42	43	30	58	21	48
										11	9	14	5	39	41	7	31	20	50
10	11	12	8	31	37	9	26	16	38
										15	12	13	7	28	23	7	25	23	13
13	5	13	5	23	15	5	23	20	15
										5	4	18	6	22	11	12	26	15	17
3	4	15	18	24	9	20	34	22	23
										6	3	12	27	25	10	24	43	31	25
8	6	10	42	26	34	23	39	39	33
										7	13	7	47	32	39	23	31	51	32

表8.预测变量PM_2.5的训练期望值(μg/m³)

55	90	42	70	43	148	224	233	85	152
										57	105	62	60	46	142	225	243	98	125
58	90	98	46	45	167	219	220	96	41
										60	80	112	37	45	173	218	217	108	45
59	71	104	31	52	183	201	41	109	49
										59	59	94	27	22	190	211	52	94	57
75	48	79	30	22	190	237	66	94	65
										68	41	85	31	28	183	254	85	119	68
65	35	76	29	150	206	250	83	126	84
										67	34	74	31	148	227	236	87	124	93

表9.预测变量PM_2.5的网络训练输出值(μg/m³)

测试样本：

表10.辅助变量温度的测试输入值(摄氏度)

15	20	18	15	14	20	16	14	21	16
										15	20	18	15	15	19	14	14	22	14
17	21	17	15	15	18	14	16	22	13
										18	21	17	15	15	17	14	18	20	12
19	19	15	15	15	16	14	21	16	11

表11.辅助变量相对湿度的测试输入值(％)

93	63	76	97	98	74	93	90	29	49
										92	59	78	98	98	79	95	94	29	62
85	56	79	98	98	83	97	64	28	68
										80	57	86	97	96	89	97	49	32	78
71	68	95	97	90	92	90	28	55	82

表12.辅助变量风速的测试输入值(级)

1	1	2	1	1	2	1	1	2	1
										1	1	1	0	1	1	0	1	2	1
1	1	2	1	1	1	1	1	2	1
										1	1	1	1	0	1	1	2	1	1
1	1	2	1	0	1	1	2	1	0

表13.辅助变量CO的测试输入值(μg/m³)

2.4240	0.9170	1.0240	1.1750	1.0240	1.5530	1.4770	1.2280	0.4540	0.9930
										2.3220	0.8980	1.4030	1.2830	1.1050	1.2880	1.1910	1.4660	0.6420	0.6870
2.0450	0.8940	1.4250	1.0920	0.9130	1.3530	1.2260	1.2290	0.4350	0.8650
										1.6020	0.9750	1.3440	1.0620	0.9130	1.4290	1.3190	0.7060	1.1900	0.9330
1.1200	0.9880	1.5770	0.9200	1.3680	1.5550	1.1220	0.4510	1.2740	0.9630

表14.辅助变量NO₂的测试输入值(μg/m³)

表15.辅助变量SO₂的测试输入值(μg/m³)

42	46	33	13	9	29	15	6	11	18
										76	85	28	10	10	38	8	9	12	23
47	97	27	9	13	29	7	12	9	29
										33	71	22	10	14	17	6	14	12	17
32	41	16	9	18	15	6	9	16	11

表16.预测变量PM_2.5的测试期望值(μg/m³)

123	82	76	127	87	121	146	94	12	27
										153	83	115	103	94	147	108	86	8	49
193	81	143	90	92	142	95	81	5	78
										179	89	147	87	83	149	102	33	15	79
121	81	137	91	89	136	94	19	21	69

表17.预测变量PM_2.5的网络测试输出值(μg/m³)

131.8512	97.2187	85.7133	127.7970	73.6662	90.4250	150.1511	81.9682	18.3922	35.9242
										139.3385	74.3293	70.5680	99.3380	76.5607	132.6773	88.7555	87.0407	16.2408	45.7539
178.0600	74.3482	114.5047	92.9579	82.8539	151.1960	96.7332	87.2388	23.3373	64.1599
										182.4573	73.7601	141.4654	85.5874	74.9998	150.6066	80.3553	71.6608	26.9252	72.1336
145.9618	85.0494	124.2878	83.6534	69.3975	161.4584	87.3236	19.0184	43.3264	64.6049

Claims (1)

1.一种基于自组织模糊神经网络的空气细颗粒物PM_2.5软测量方法，其特征包括以下步骤：

(1)确定辅助变量和预测变量；

①变量数据的标准化处理：将空气污染物浓度及气象变量的上一小时数据和PM_2.5浓度的下一小时数据一一对应，整理出L组数据，然后以空气污染物浓度和气象变量为自变量X，PM_2.5浓度为因变量Y：

X＝[x₁,x₂,L,x₇]^T,Y＝y(1)

其中，x₁,x₂,L,x₇分别表示温度、相对湿度、风速、CO、NO₂、O₃、SO₂浓度，y表示PM_2.5浓度，温度的单位是摄氏度，风速的单位是级，各空气污染物浓度的单位是μg/m³，相对湿度用百分数表示，没有单位；数据标准化处理的计算公式如下：

$g_1=\frac{D_1-D_{1,min}}{D_{1,max}-D_{1,min}}---\left(2\right)$

其中，g₁是温度的上一小时实测数据数组标准化后得到的数组，D₁是标准化前温度的实测数据数组，D_1,min是数组D₁中的最小值，D_1,max是数组D₁中的最大值；相对湿度、风速、CO、NO₂、O₃、SO₂浓度的上一小时实测数据数组和PM_2.5浓度的下一小时实测数据数组分别是D₂，D₃，D₄，D₅，D₆，D₇和D₈，找出各数组内部的最大值和最小值，按照式(2)相同的处理方式得到标准化数组分别是：g₂，g₃，g₄，g₅，g₆，g₇和g₈；g₁，g₂，g₃，g₄，g₅，g₆和g₇构成了自变量X的标准化观测数据阵E₀，g₈是因变量Y的标准化观测数据阵F₀；

②计算的最大特征值对应的单位特征向量W₁，求取自变量组的第一个成分其得分向量d₁＝E₀W₁，残差矩阵其中α₁和β₁是模型效应负荷量，且

③计算的最大特征值对应的单位特征向量W₂，求取自变量组的第二个成分其得分向量d₂＝E₁W₂，残差矩阵其中α₂和β₂是模型效应负荷量，且将和代入Y＝t₁β₁+t₂β₂中得到提取2个成分时因变量Y关于自变量X的偏最小二乘回归方程式：

y＝a₁x₁+L+a₇x₇(3)

其中，a₁,La₇分别是提取2个成分时自变量x₁,L,x₇对应的回归系数；

④利用交叉有效性检验提取2个成分时模型精度是否满足要求：首先每次从E₀和F₀中舍去第q组数据，利用剩下的L-1组数据建立提取2个成分时的偏最小二乘回归式，并求出在舍去的数据处的预测值y_(-q)(2)，真实值是y_q；对q＝1,2,…,L重复以上操作，即得提取2个成分时因变量Y的预测误差平方和PRESS(2)：

$PRESS\left(2\right)={\mathrm{\&Sigma;}}_{q=1}^L{(y_q-y_{(-q)}(2\left)\right)}^2---\left(4\right)$

然后利用E₀和F₀拟合含1个成分的偏最小二乘回归方程式，记第q组数据的预测值为y_(q)(1)，则提取1个成分时Y的误差平方和SS(1)为：

$SS\left(1\right)={\mathrm{\&Sigma;}}_{q=1}^L{(y_q-y_{\left(q\right)}(1\left)\right)}^2---\left(5\right)$

定义提取2个成分时的交叉有效性为Q²(2)＝1-PRESS(2)/SS(1)，若Q²(2)＜1-0.95²，则模型精度满足要求，停止提取成分，否则继续提取成分；

⑤通过上述步骤得到因变量Y关于自变量X的偏最小二乘回归方程式后，根据回归方程式中回归系数的绝对值，除去绝对值最小的自变量，将剩下的自变量选作辅助变量，记为r₁,r₂,L,r₆，PM_2.5浓度作为预测变量，是自组织模糊神经网络的期望输出，记为o_d，辅助变量和预测变量的准化观测数据阵分别记为E和F，选取E和F的前I组数据作为训练样本，后I′组数据作为测试样本；

(2)设计用于PM_2.5预测的自组织模糊神经网络的初始拓扑结构；本发明用于PM_2.5预测的自组织模糊神经网络共四层：输入层、RBF层、规则化层和输出层；输入是选取的辅助变量，输出是PM_2.5浓度，确定其6-P-P-1的初始连接方式，即输入层神经元数和辅助变量的个数同是6，RBF层神经元数是P，规则化层神经元数是P，其中P是正整数，在[1-10]之间取值，输出层神经元数是1；自组织模糊神经网络的中心、宽度和权值的初始值设定在(0,1)之间；自组织模糊神经网络的输入表示为r＝[r₁,r₂,L,r₆]^T，[r₁,r₂,L,r₆]^T是[r₁,r₂,L,r₆]的转置，期望输出是o_d；第k组训练样本表示为r(k)＝[r₁(k),r₂(k),L,r₆(k)]^T，k＝1,2,…,I，用第k组训练样本训练自组织模糊神经网络时，自组织模糊神经网络的各层输出依次为：

①输入层：该层有6个神经元：

u_i(k)＝r_i(k)(6)

其中，u_i(k)是输入层第i个神经元的输出，该层的输入向量为r(k)＝[r₁(k),r₂(k),L,r₆(k)]^T，r₁(k),r₂(k),L,r₆(k)分别表示对应的辅助变量；

②RBF层：该层选取RBF神经元的高斯函数作为隶属函数对输入变量进行模糊化处理，RBF层有P个神经元，该层第j个神经元的输出为：

其中，c_ij(k)，σ_ij(k)分别是模糊神经网络的中心和宽度；

③规则化层：该层神经元数目与RBF层相同，该层第l个神经元的输出v_l(k)为：

其中，是RBF层第l个神经元的输出；

④输出层：该层有1个神经元，该层的输出表示PM_2.5浓度，如下式所示：

$o\left(k\right)={\mathrm{\&Sigma;}}_{l=1}^P{w}_l\left(k\right)v_l\left(k\right)---\left(9\right)$

其中，w_l(k)是规则化层第l个神经元与输出层神经元之间的连接权值；式(6)-式(9)中，i＝1,2,…,6，j＝1,2,…,P，l＝1,2,…,P；自组织模糊神经网络的训练均方根误差RMSE为：

$RMSE={\mathrm{\&Sigma;}}_{k=1}^I\frac{{\left(o_d\right(k)-o(k\left)\right)}^2}{I}---\left(10\right)$

其中，o_d(k)和o(k)分别是第k组训练样本的期望输出和用第k组训练样本训练模糊神经网络时的网络输出，训练模糊神经网络的目的是使式(10)定义的训练RMSE达到期望值；

(3)用训练样本训练自组织模糊神经网络；在训练过程中，利用频域上的模型输出敏感度分析方法确定模糊神经网络规则化层神经元输出对网络输出的贡献率，根据贡献率大小增删规则化层神经元，以此来自动调整模糊神经网络的结构，并利用梯度下降算法训练模糊神经网络的中心、宽度和权值，直到网络精度满足信息处理的需求，具体为：

①利用训练样本和梯度下降算法训练给定初始模糊神经网络，网络每产生一组输出参数更新一次，对I组训练样本进行反复训练直到满足结构修改步数，结构修改步数取值为20-30后转向②，设置最大迭代步数是1000，期望训练RMSE设定为0.01；梯度下降算法参数调整公式如下：

$\left\{\begin{array}{c}{c}_{ij}(k+1)=c_{ij}\left(k\right)-\mathrm{\&eta;}\&part;e\left(k\right)/\&part;c_{ij}\left(k\right)\\ {\mathrm{\&sigma;}}_{ij}(k+1)={\mathrm{\&sigma;}}_{ij}\left(k\right)-\mathrm{\&eta;}\&part;e\left(k\right)/\&part;{\mathrm{\&sigma;}}_{ij}\left(k\right)\\ w_l(k+1)=w_l\left(k\right)-\mathrm{\&eta;}\&part;e\left(k\right)/\&part;w_l\left(k\right)\end{array}\right.---\left(11\right)$

其中，c_ij(k)、σ_ij(k)和w_l(k)分别是第k组训练样本训练自组织模糊神经网络时网络的中心、宽度和权值，c_ij(k+1)、σ_ij(k+1)和w_l(k+1)分别是第k+1组训练样本训练自组织模糊神经网络时网络的中心、宽度和权值，目标函数e(k)＝(o_d(k)-o(k))²/2，o_d(k)和o(k)分别是第k组训练样本的期望输出和用第k组训练样本训练自组织模糊神经网络时的网络输出，η是学习率，取值是0.005，②在频域中计算规则化层第h个神经元的输出v_h对网络输出o的贡献率ST_h，h＝1,…,P如下式所示：

${\mathrm{ST}}_h=\frac{S_h}{{\mathrm{\&Sigma;}}_{j=1}^P{S}_j}---\left(12\right)$

$S_h=\frac{A_{{\mathrm{\&omega;}}_h}^2+B_{{\mathrm{\&omega;}}_h}^2}{(A_{{\mathrm{\&omega;}}_h}^2+B_{{\mathrm{\&omega;}}_h}^2)+{\mathrm{\&Sigma;}}_{\mathrm{\&omega;}=1}^{\max \left({\mathrm{\&omega;}}_{~h}\right)}(A_{\mathrm{\&omega;}}^2+B_{\mathrm{\&omega;}}^2)}---\left(13\right)$

其中，S_h是v_h对o的总效应，A_ω和B_ω及和分别是o的傅里叶展开式f(s)在频率ω和ω_h处的傅里叶系数，且 ω_h是v_h的基频，max(ω_～h)是除v_h之外的其他输入参数基频的最大值，取ω_h＝2Mmax(ω_～h)，M是干扰因子，取值是4，f(s)和v_h分别用式(14)和式(15)表示：

$f\left(s\right)={\mathrm{\&Sigma;}}_{\mathrm{\&omega;}=-\mathrm{\&infin;}}^{\mathrm{\&infin;}}\left(A_{\mathrm{\&omega;}}cos\right(\mathrm{\&omega;}s)+B_{\mathrm{\&omega;}}\sin (\mathrm{\&omega;}s\left)\right)---\left(14\right)$

$v_h\left(s\right)=\frac{p_h+b_h}{2}+\frac{p_h-b_h}{\mathrm{\&pi;}}\arcsin \left(\sin \right({\mathrm{\&omega;}}_{h}s)---\left(15\right)$

其中，p_h和b_h分别是v_h的最大值和最小值；

③如果ST_h>＝λ₁，，λ₁＝0.3，则对规则化层第h个神经元进行分裂；为了降低网络结构调整对网络误差的影响，采用下式设定分裂得到的新的神经元的初始参数：

$\left\{\begin{array}{c}{c}_{\&CenterDot;new1}=c_{\&CenterDot;new2}=c_{\&CenterDot;h}\left(k\right),{\mathrm{\&sigma;}}_{\&CenterDot;new1}={\mathrm{\&sigma;}}_{\&CenterDot;new2}={\mathrm{\&sigma;}}_{\&CenterDot;h}\left(k\right)\\ w_{new1}={\mathrm{\&mu;w}}_h\left(k\right),w_{new2}=(1-\mathrm{\&mu;})w_h\left(k\right)\end{array}\right.---\left(16\right)$

其中，new1和new2是两个新的神经元，c_·new1、σ_·new1和w_new1分别是神经元new1的中心向量、宽度向量和其与输出层神经元间的连接权值；c_·new2、σ_·new2和w_new2分别是神经元new2的中心向量、宽度向量和权值；c_·h(k)、σ_·h(k)和w_h(k)分别是网络结构调整前神经元h的中心向量、宽度向量和权值，μ服从标准正态分布；

④如果ST_h<λ₂，λ₂＝0.05，删除该神经元，同时对网络参数进行调整，如下式所示：

$\left\{\begin{array}{c}{c}_{\&CenterDot;h-m}=c_{\&CenterDot;h-m}\left(k\right),{\mathrm{\&sigma;}}_{\&CenterDot;h-m}={\mathrm{\&sigma;}}_{\&CenterDot;h-m}\left(k\right)\\ w_{h-m}=w_{h-m}\left(k\right)+w_h\left(k\right)v_h\left(k\right)/v_{h-m}\left(k\right)\end{array}\right.---\left(17\right)$

其中，神经元h-m是与神经元h欧氏距离最小的规则化层神经元，且ST_h-m>＝λ₂，c_.h-m、σ_.h-m和w_h-m分别是网络修剪后神经元h-m的中心向量、宽度向量和权值，c_.h-m(k)、σ_.h-m(k)和w_h-m(k)分别是网络修剪前神经元h-m的中心向量、宽度向量和权值，w_h(k)是网络修剪前神经元h的权值，v_h(k)和v_h-m(k)分别是网络修剪前神经元h和神经元h-m的输出；

⑤利用梯度下降算法更新网络参数，若在某步训练过程中训练RMSE<＝0.01或算法迭代了1000步时则停止计算，否则达到结构修改步数后转向②；

(4)对测试样本进行检测；利用测试样本对训练好的自组织模糊神经网络进行测试，自组织模糊神经网络的输出即为PM_2.5的预测结果。