CN104680015A - 一种基于快速相关向量机的污水处理在线软测量方法 - Google Patents

Abstract

本发明公开了一种基于快速相关向量机的污水处理在线软测量方法，包含以下顺序的步骤：通过快速边际似然算法估计超参数，得到模型的权重值和样本偏差值；然后建立快速相关向量机在线预测模型，对模型参数寻优，实现了污水中BOD的精确快速测量。本发明的测量方法，能够满足实时性的要求，建立最优预测模型，预测精度得到了提高，效果显著，性能得到了改善，快速相关向量机建立的污水水质在线软测量模型预测精度高、泛化能力强、更新时间短，对于节省污水处理厂运营费用，实时反映污水水质状况，对污水处理自动控制***具有重要的意义。

Description

一种基于快速相关向量机的污水处理在线软测量方法

技术领域

本发明涉及污水处理领域，特别涉及一种基于快速相关向量机的污水处理在线软测量方法。

背景技术

污水处理是经济发展和水资源保护不可或缺的组成部分。随着国民经济的快速增长，污水排放量也大大增加，而污水处理厂太少，处理周期太长，远远达不到国家对环境保护的要求。同时国家对环境保护的投入加大，污水处理技术越来越受到更多的关注。国家发展规划中明确提出要研发并推广低能耗、有效的污水处理技术。

污水排放标准中，衡量是否达标的参数指标有：化学需氧量COD、生化需氧量BOD、氨氮、磷、固体悬浮物等。其中生化需氧量BOD和化学需氧量COD反映水被有机污染的程序，BOD/COD的比率反映出了污水的生物降解能力。这两个参数的测量对控制污水处理具有非常重要的价值。化学需氧量COD是指，水样在一定条件下，以氧化1升水样中还原性物质所消耗的氧化剂的量为指标，折算成每升水样全部被氧化后，需要的氧的毫克数，以mg/L表示。生化需氧量BOD是指微生物在一定的温度和时间条件下分解氧化有机物所消耗的溶解氧量，以mg/L表示。

现在的污水处理一般都采用稀释法、传感器等测量污水中BOD、COD的浓度，但由于分析测定这两个指标的周期较长，测量中时常出现误差，不能及时反应污水处理的现场情况，因而污水控制***存在着较大的延时，不能发挥其最佳的性能。本发明提出一种新的测量BOD的软测量方法，通过快速边际似然算法来对相关向量机的训练过程进行改进，能够更快地使超参数达到稳定值，从而得到权重值和偏差值，并构建在线的快速相关向量机软测量模型，实现对污水处理出水BOD的测量。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供一种基于快速相关向量机的污水处理在线软测量方法。

本发明的目的通过以下的技术方案实现：

一种基于快速相关向量机的污水处理在线软测量方法，包含以下顺序的步骤：

A、通过快速边际似然算法估计超参数，得到模型的权重值和样本偏差值；

B、然后建立快速相关向量机在线预测模型，对模型参数寻优，实现了污水中BOD的精确快速测量。

所述的基于快速相关向量机的污水处理在线软测量方法，具体包含以下步骤：

S1.剔除输入和输出的数据中的异常点，由于各输入变量量纲的不同，对其进行归一化处理，归一化到[0,1]区间中；

S2.给定污水数据集{(x_n,t_n),n＝1,2,...,N},x_n∈R^d,t_n∈R,N是样本数,为简便起见,只考虑标量目标函数,我们遵循标准的概率公式,假定:

t_n＝y(x_n；w)+ε_n    (1)

其中,y(·)是非线性函数，ε_n是均值为0,方差为σ²的高斯噪声，即因此有t_n～N(y(x_n),σ²)，这表明t_n满足均值y(x_n)为方差为σ²的高斯噪声分布，与支持向量机类似，函数y(x)定义为

$y(x;w)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i}K(x,x_i)+w_0---\left(2\right)$

其中，由公式φ_i(x)＝K(x,x_i)来确定基函数，其核被训练向量参数化，假定t_n是相互独立的，则整个训练集的似然函数可写为

$p\left(t\right|w,{\mathrm{\σ}}^2)={\left({2\mathrm{\π\σ}}^2\right)}^{-N/2}\exp (-\frac{{\left|\right|t-\mathrm{\Φw}\left|\right|}^2}{{2\mathrm{\σ}}^2})---\left(3\right)$

式中t＝[t₁,t₂,...,t_N]^T,w＝[w₀,w₁,...,w_M]^T,Φ是一个N×(N+1)的设计矩阵，Φ＝[φ₁,φ₂,...,φ_M]是组非线性基函数，φ(x_n)＝[1,K(x_n,x₁),K(x_n,x₂),...,K(x_n,x_N)]^T；

由于在模型中有和训练样本差不多的参数个数，从(3)式中得到的w和σ²的最大似然估计值有可能导致模型过拟合；为了避免过度拟合，通常的做法是给参数强加一些限制条件；在这里,我们从贝叶斯概率框架出发,通过定义一个先验概率分布来限制参数w和σ²；

这里选择一个比较平滑的函数,定义w的先验概率分布为零均值的高斯分布:

$p\left(w\right|\mathrm{\α})=\underset{j=0}{\overset{N}{\mathrm{\Π}}}N\left(w_j\right|0,{\mathrm{\α}}_j^{-1})---\left(4\right)$

式(4)中:超参数α＝[α₀,α₁,...,α_N]^T，更重要的是，每个独立的超参数α_j只与其对应的权值w_j相关；通过这个限制条件，在经过大量的污水数据学习后,大部分超参数会趋近于无穷大,而与其对应的权值为0,从而使RVM具有较高的稀疏性；

现在已经定义了先验概率，从贝叶斯规则来看，对于给定的数据中未知数据，贝叶斯推理通过计算后验概率处理：

$p(w,\mathrm{\α},{\mathrm{\σ}}^2|t)=\frac{p\left(t\right|w,\mathrm{\α},{\mathrm{\σ}}^2)p(w,\mathrm{\α},{\mathrm{\σ}}^2)}{p\left(t\right)}---\left(5\right)$

给定一个测试点x_*,相应的污水出水水质预测值t_*的预测分布为

p(t_*|t)＝∫p(t_*|w,α,σ²)p(w,α,σ²|t)dwdαdσ²    (6)

根据贝叶斯公式,利用样本似然函数(4)和w先验分布(5)可得w的后验分布为

$p\left(w\right|\mathrm{t\α},\mathrm{\β})=\frac{p\left(w\right|\mathrm{\α})p\left(t\right|w,\mathrm{\β})}{p\left(t\right|\mathrm{\α},\mathrm{\β})}---\left(7\right)$

我们把后验概率分解为

p(w,α,σ²|t)＝p(|w|t,α,σ²)p(α,σ²|t)    (8)

因此对权重的后验概率分布为

$p\left(w\right|t,\mathrm{\α},{\mathrm{\σ}}^2)=\frac{p\left(t\right|w,{\mathrm{\σ}}^2)p\left(w\right|\mathrm{\α})}{p\left(t\right|\mathrm{\α},{\mathrm{\σ}}^2)}={\left(2\mathrm{\π}\right)}^{-(N+1)/2}{\left|\mathrm{\Σ}\right|}^{-1/2}\exp \{-\frac{1}{2}{(w-\mathrm{\μ})}^T{\mathrm{\Σ}}^{-1}(w-\mathrm{\μ})\}---\left(9\right)$

其协方差为∑＝(σ^-2Φ^TΦ+A)^-1    (10)

平均值为μ＝σ^-2∑Φ^Tt    (11)

其中矩阵A＝diag(α₀,α₁,...,α_N)

${\left({\mathrm{\α}}_i\right)}^{\mathrm{new}}=\frac{{\mathrm{\γ}}_i}{{\mathrm{\μ}}_i^2}---\left(12\right)$

${\left({\mathrm{\σ}}^2\right)}^{\mathrm{new}}=\frac{{\left|\right|t-\mathrm{\Φ\μ}\left|\right|}^2}{N-{\mathrm{\Σ}}_i{\mathrm{\γ}}_i}---\left(13\right)$

其中γ_i≡1-α_i∑_ii，∑_ii为协方差矩阵∑的第i个对角元素；最后通过公式(10)到公式(13)的迭代推理运算得到超参数α和方差σ²的估计值；输出的污水水质预测值为y_*＝μ^Tφ(x_*)，x_*是污水处理过程输入值；

S3.用快速边际似然算法估计超参数

针对相关向量机计算时间复杂度大、内存开销大的问题，采用了一种快速边际似然算法；它在训练过程中从空集开始动态地扩充基矩阵Φ，从而增大边际似然函数，或者去掉基矩阵Φ冗余的列来增大目标函数；

S4.对待预测的污水样本数据进行预测：将入水数据作为训练好的相关向量机软测量模型的输入，模型的输出即为出水BOD的预测结果。

所述的步骤S3，具体包含以下顺序的步骤：

相关向量机是通过最大化边际似然函数p(t|α,σ²)的方法确定超参数α和方差σ²的，等价于最大化为其对数；记L(α)＝log[p(t|α,σ²)]，整理有

$\begin{array}{c}L\left(\mathrm{\α}\right)=\log \left[p\left(t\right|\mathrm{\α},{\mathrm{\σ}}^2)\right]=\log [\∫p\left(t\right|w,{\mathrm{\σ}}^2)p\left(w\right|\mathrm{\α})]\\ =-\frac{1}{2}[N\log 2\mathrm{\π}+\log |C|+t^T{C}^{-1}t]\end{array}---\left(14\right)$

其中C＝σ²I+ΦA^-1Φ^T,t＝[t1,t2,…,t_N]^T；

为了便于最大化L(α)，对矩阵C进行等价变换，如下：

$\begin{array}{c}C={\mathrm{\σ}}^{2}I+{\mathrm{\ΦA}}^{-1}{\mathrm{\Φ}}^T={\mathrm{\σ}}^{2}I+\underset{m\≠i}{\mathrm{\Σ}}{\mathrm{\α}}_m^{-1}{\mathrm{\φ}}_m{\mathrm{\φ}}_m^T+{\mathrm{\α}}_i^{-1}{\mathrm{\φ}}_i{\mathrm{\φ}}_i^T\\ =C_{-i}+{\mathrm{\α}}_i^{-1}{\mathrm{\φ}}_i{\mathrm{\φ}}_i^T\end{array}---\left(15\right)$

其中此矩阵表示当α_i＝∞时，相应的基向量φ_i被移除后样本对应的协方差矩阵，根据矩阵相关性质整理可得

$\begin{array}{c}\left|C\right|=\left|C_{-i}\right||1+{\mathrm{\α}}_i^{-1}{\mathrm{\φ}}_i^T{C}_{-i}^{-1}{\mathrm{\φ}}_i|\\ C^{-1}=C_{-i}^{-1}-\frac{C_{-i}^{-1}{\mathrm{\φ}}_i{\mathrm{\φ}}_i^T{C}_{-i}^{-1}}{{\mathrm{\α}}_i+{\mathrm{\φ}}_i^T{C}_{-i}^{-1}{\mathrm{\φ}}_i}\end{array}---\left(16\right)$

因此公式(14)可以改写为

$\begin{array}{c}L\left(\mathrm{\α}\right)=-\frac{1}{2}[N\log 2\mathrm{\π}+\log |C|+t^T{C}_{-i}t-\log {\mathrm{\α}}_i+\log ({\mathrm{\α}}_i+{\mathrm{\φ}}_i^T{C}_{-i}^{-1}{\mathrm{\φ}}_i)-\frac{{\left({\mathrm{\φ}}_i^T{C}_{-i}^{-1}t\right)}^2}{{\mathrm{\α}}_i+{\mathrm{\φ}}_i^T{C}_{-i}^{-1}{\mathrm{\φ}}_i}]\\ =L\left({\mathrm{\α}}_{-i}\right)+\frac{1}{2}[\log {\mathrm{\α}}_i-\log ({\mathrm{\α}}_i+s_i)+\frac{{\left(q_i\right)}^2}{{\mathrm{\α}}_i+s_i}\\ =L\left({\mathrm{\α}}_{-i}\right)+l\left({\mathrm{\α}}_i\right)\end{array}---\left(17\right)$

注意L(α_-i)表示为当α_i＝∞时，相应的基本向量φ_i被移除后所对应的边界似然函数的对数，而l(α_i)表示边界似然的对数函数中只与α_i有关的独立部分，上式对α_i求偏导有

$\frac{\∂L\left(\mathrm{\α}\right)}{{\∂\mathrm{\α}}_i}=\frac{\∂l\left({\mathrm{\α}}_i\right)}{{\∂\mathrm{\α}}_i}=\frac{1}{2}[\frac{1}{{\mathrm{\α}}_i}-\frac{1}{{\mathrm{\α}}_i+{\mathrm{\φ}}_i^T{C}_{-i}^{-1}{\mathrm{\φ}}_i}-\frac{{\left({\mathrm{\φ}}_i^T{C}_{-i}^{-1}t\right)}^2}{{\mathrm{\α}}_i+{\mathrm{\φ}}_i^T{C}_{-i}^{-1}{\mathrm{\φ}}_i}]---\left(18\right)$

记 $S_i={\mathrm{\φ}}_i^T{C}_{-i}^{-1}{\mathrm{\φ}}_i,Q_i={\mathrm{\φ}}_i^T{C}_{-i}^{-1}t---\left(19\right)$

所以公式(18)可改写为

$\frac{\∂L\left(\mathrm{\α}\right)}{{\∂\mathrm{\α}}_i}=\frac{{\mathrm{\α}}_i^{-1}{S}_i^2-(Q_i^2-S_i)}{2{({\mathrm{\α}}_i+S_i)}^2}---\left(20\right)$

令公式(20)等于零，考虑到α_i是方差值必须为正，所以当时有

${\mathrm{\α}}_i=\frac{S_i^2}{Q_i^2-S_i}---\left(21\right)$

对L(α)关于α_i求二阶偏导有

$\begin{array}{c}\frac{{\∂}^{2}L\left(\mathrm{\α}\right)}{{\∂\mathrm{\α}}_i^2}=\frac{{-\mathrm{\α}}_i^{-2}{S}_i^2{({\mathrm{\α}}_i+S_i)}^2-2({\mathrm{\α}}_i+S_i)[{\mathrm{\α}}_i^{-1}{S}_i^2-(Q_i^2-S_i)]}{2{({\mathrm{\α}}_i+S_i)}^4}\\ =\frac{{\mathrm{\α}}_i^{-2}{S}_i^2}{2{({\mathrm{\α}}_i+S_i)}^2}-\frac{S_i^2}{{\mathrm{\α}}_i{({\mathrm{\α}}_i+S_i)}^3}-\frac{[{\mathrm{\α}}_i^{-1}{S}_i^2-(Q_i^2-S_i)]}{{({\mathrm{\α}}_i+S_i)}^3}\end{array}---\left(22\right)$

综合公式(20)和(21)进行分析可知

$\frac{{\∂}^{2}L\left(\mathrm{\α}\right)}{{\∂\mathrm{\α}}_i^2}{|}_{{\mathrm{\α}}_i=\frac{S_i^2}{Q_i^2-S_i}}=\frac{{-S}_i^2}{{2\mathrm{\α}}_i^2{({\mathrm{\α}}_i+S_i)}^2}---\left(23\right)$

所以当时，公式(23)左边的表达式是恒小于零的，并对以上推导公式分析可得，L(α)有唯一最大值点为

${\mathrm{\α}}_i=\left\{\begin{array}{cc}\frac{S_i^2}{Q_i^2-S_i}& Q_i^2>S_i\\ \∞& Q_i^2\≤S_i\end{array}\right.---\left(24\right)$

由此得到L(α)的最大值。

所述的贝叶斯L(α)通过以下方法最大化：

a、当基向量φ_i在模型中，即α_i<∞，但有则将φ_i从模型中删除，即令α_i＝∞，这样可以增大贝叶斯L(α)；

b、当基向量φ_i在模型中，即α_i＝∞，但有则将φ_i添加到模型中并利用公式(24)更新α_i，这样可以增大贝叶斯L(α)；

c、当基向量φ_i在模型中，即α_i<∞，但有则用公式(24)更新α_i，这样可以增大贝叶斯L(α)。

所述的快速相关向量机，其回归基本算法步骤如下：

Ⅰ、初始化σ²；

Ⅱ、用单个基向量φ_i初始化α_i，由公式(24)分析整理可得并设置其他的α_m(m≠i)为无穷大；

Ⅲ、计算∑、μ并对所有M个基函数φ_m初始化S_m和Q_m；

Ⅳ、从所有M个基函数φ_m集合中选择候选的基向量φ_i；

Ⅴ、计算 ${\mathrm{\&theta;}}_i\&equiv;Q_i^2-S_i;$

Ⅵ、若θ_i>0且α_i<∞(基向量φ_i在模型中)，重新估计α_i；

Ⅶ、若θ_i>0且α_i＝∞(基向量φ_i不在模型中)，添加φ_i到模型中并重新估计α_i；

Ⅷ、若θ_i≤0且α_i<∞，删除φ_i并设置α_i＝∞；

Ⅸ、估计噪声方差其中N为数据个数，M为基函数个数；

Ⅹ、重新计算协方差矩阵∑，权重矩阵μ以及相应迭代过程中的S_m和Q_m；

Ⅺ、若收敛或者达到最大迭代次数，则保存权重值及偏差值，此次训练结束；否则转步骤IV继续训练。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明采用一种基于快速相关向量机的污水出水水质在线软测量模型。首先建立快速相关向量机离线模型，然后根据工况实时更新模型，保留模型中重要的信息数据，并运用快速边际似然算法加快模型的学习速度，满足实时性的要求，建立最优预测模型，预测精度得到了提高，效果显著，性能得到了改善。快速相关向量机建立的污水水质在线软测量模型预测精度高、泛化能力强、更新时间短，对于节省污水处理厂运营费用，实时反映污水水质状况，对污水处理自动控制***具有重要的意义。

2、软测量的基本思想是把自动控制理论与生产过程知识有机结合起来，应用计算机技术，针对难于测量或暂时不能测量的重要变量(或称之为主导变量)，选择另外一些容易测量的变量(或称之为辅助变量)，通过构成某种数学关系来推断和估计，以软件来代替硬件(传感器)功能。这种方法响应迅速，能够连续给出主导变量信息，且具有投资低、维护保养简单等优点。将软测量技术用于污水处理过程，能降低污水处理厂的能耗，设备维护费等，并且对于传感器测量产生的误差有较正作用。但是传统的离线软测量模型是根据大量的数据一次训练建立的，之后不会再改变，这样就可能造成对新数据适应度不那么强。在线的软测量模型则是先建立离线模型，然后根据工况的变化模型也作相应的改变，始终与最新的数据保持近似同步更新，这样预测精度也得到了相应提高。

附图说明

图1是本发明模型在线模型出水BOD结果拟合图；

图2是本发明模型离线训练权重值分布图；

图3是本发明模型离线模型出水BOD结果拟合图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

本发明采用快速相关向量机建立离线模型；利用快速边际似然算法改进模型的训练过程，可以加快模型的在线训练速度，使超参数更快地达到稳定值。其算法步骤如下

算法步骤：

S1.初始化σ²；

S2.用单个基向量φ_i初始化α_i，由公式(24)分析整理可得并设置其他的α_m(m≠i)为无穷大；

S3.计算∑、μ并对所有M个基函数φ_m初始化S_m和Q_m；

S4.从所有M个基函数φ_m集合中选择候选的基向量φ_i；

S5.计算 ${\mathrm{\&theta;}}_i\&equiv;Q_i^2-S_i;$

S6.若θ_i>0且α_i<∞(基向量φ_i在模型中)，重新估计α_i；

S7.若θ_i>0且α_i＝∞(基向量φ_i不在模型中)，添加φ_i到模型中并重新估计α_i；

S8.若θ_i≤0且α_i<∞，删除φ_i并设置α_i＝∞；

S9.估计噪声方差其中N为数据个数，M为基函数个数；

S10.重新计算协方差矩阵∑，权重矩阵μ以及相应迭代过程中的S_m和Q_m；

S11.若收敛或者达到最大迭代次数，则保存权重值及偏差值，此次训练结束；否则转S4继续训练。

污水数据来源于加州大学数据库(UCI)中的污水数据。BOD是反映水体被有机物污染程度的综合指标。与悬浮固体浓度，进出水的化学需氧量，进水的生物需氧量，流量，PH值等变量密切相关温度。建模所需辅助变量为可降解固体浓度RD-SED-G，悬浮固体浓度RD-SS-G，生化需氧量RD-DBO-G，化学需氧量RD-DQO-G，初沉池中的生化需氧量RD-DBO-P，悬浮固体浓度RD-SS-P，二沉池中的生化需氧量RD-DBO-S，化学需氧量RD-DQO-S，入水中的生化需氧量DBO，化学需氧量DQO，二级处理中的化学需氧量DQO，生化需氧量DBO，悬浮固体浓度SS，PH值PH-S，可降解固体浓度SED，出水的化学需氧量DQO-S，可降解固体浓度SED-S，悬浮固体浓度SS-S，PH值PH-S。由上可知，输入属性19个，输出属性1个。选取处理后的400组数据，其中200组用于训练模型，200组用作新数据测试模型精度。

(2)RVM建模步骤如下：

污水数据集{(x_n,t_n),n＝1,2,...,N},x_n∈R^d,t_n∈R,N是样本数,假定：

t_n＝y(x_n；w)+ε_n    (1)

其中,y(·)是非线性函数，ε_n是均值为0，方差为σ²的高斯噪声，即因此有t_n～N(y(x_n),σ²)，函数y(x)定义为

$y(x;w)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{i}K(x,x_i)+w_0---\left(2\right)$

在公式中确定基函数,其核被训练向量参数化φ_i(x)＝K(x,x_i).假定t_n是相互独立的,则整个训练集的似然函数可写为

$p\left(t\right|w,{\mathrm{\&sigma;}}^2)={\left({2\mathrm{\&pi;\&sigma;}}^2\right)}^{-N/2}\exp (-\frac{{\left|\right|t-\mathrm{\&Phi;w}\left|\right|}^2}{{2\mathrm{\&sigma;}}^2})---\left(3\right)$

式中t＝[t₁,t₂,...,t_N]^T,w＝[w₀,w₁,...,w_M]^T,Φ是一个N×(N+1)的设计矩阵,Φ＝[φ₁,φ₂,...,φ_M]是组非线性基函数,φ(x_n)＝[1,K(x_n,x₁),K(x_n,x₂),...,K(x_n,x_N)]^T.

由于在模型中有和训练样本差不多的参数个数,从(3)式中得到的w和σ²的最大似然估计值有可能导致模型过拟合.为了避免过度拟合,通常的做法是给参数强加一些限制条件.在这里,我们从贝叶斯概率框架出发,通过定义一个先验概率分布来限制这些参数,

这里选择一个比较平滑的函数,定义w的先验概率分布为零均值的高斯分布:

$p\left(w\right|\mathrm{\&alpha;})=\underset{j=0}{\overset{N}{\mathrm{\&Pi;}}}N\left(w_j\right|0,{\mathrm{\&alpha;}}_j^{-1})---\left(4\right)$

式(4)中:超参数α＝[α₀,α₁,...,α_N]^T,更重要的是,每个独立的超参数α_j只与其对应的权值w_j相关.通过这个限制条件,在经过大量的污水数据学习后,大部分超参数会趋近于无穷大,而与其对应的权值为0,从而使RVM具有较高的稀疏性。

从贝叶斯规则来看，对于给定的数据中未知数据,贝叶斯推理通过计算后验概率处理，

$p(w,\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2|t)=\frac{p\left(t\right|w,\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2)p(w,\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2)}{p\left(t\right)}---\left(5\right)$

给定一个测试点x_*，相应的污水出水水质预测值t_*的预测分布为

p(t_*|t)＝∫p(t_*|w,α,σ²)p(w,α,σ²|t)dwdαdσ²    (6)

根据贝叶斯公式,利用样本似然函数(4)和w先验分布(5)可得w的后验分布为

$p\left(w\right|\mathrm{t\&alpha;},\mathrm{\&beta;})=\frac{p\left(w\right|\mathrm{\&alpha;})p\left(t\right|w,\mathrm{\&beta;})}{p\left(t\right|\mathrm{\&alpha;},\mathrm{\&beta;})}---\left(7\right)$

我们把后验概率分解为

p(w,α,σ²|t)＝p(|w|t,α,σ²)p(α,σ²|t)    (8)

因此对权重的后验概率分布为

$p\left(w\right|t,\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2)=\frac{p\left(t\right|w,{\mathrm{\&sigma;}}^2)p\left(w\right|\mathrm{\&alpha;})}{p\left(t\right|\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2)}={\left(2\mathrm{\&pi;}\right)}^{-(N+1)/2}{\left|\mathrm{\&Sigma;}\right|}^{-1/2}\exp \{-\frac{1}{2}{(w-\mathrm{\&mu;})}^T{\mathrm{\&Sigma;}}^{-1}(w-\mathrm{\&mu;})\}---\left(9\right)$

其协方差为

∑＝(σ^-2Φ^TΦ+A)^-1    (10)

平均值为

μ＝σ^-2∑Φ^Tt    (11)

其中矩阵A＝diag(α₀,α₁,...,α_N)

${\left({\mathrm{\&alpha;}}_i\right)}^{\mathrm{new}}=\frac{{\mathrm{\&gamma;}}_i}{{\mathrm{\&mu;}}_i^2}---\left(12\right)$

${\left({\mathrm{\&sigma;}}^2\right)}^{\mathrm{new}}=\frac{{\left|\right|t-\mathrm{\&Phi;\&mu;}\left|\right|}^2}{N-{\mathrm{\&Sigma;}}_i{\mathrm{\&gamma;}}_i}---\left(13\right)$

其中γ_i≡1-α_i∑_ii,∑_ii为协方差矩阵∑的第i个对角元素.最后通过(10)到(13)式的迭代推理运算得到超参数α和方差σ²的估计值。

输出的污水水质预测值为y_*＝μ^Tφ(x_*)，x_*是污水处理过程输入值。离线模型预测结果如图3所示。权值分布如图2所示。

(3)假设最新的污水输入属性为x_new，实际的出水BOD值为y_new，则出水BOD的在线更新算法如下步骤：

S1.根据历史数据使用上节所述快速相关向量机回归基本算法步骤建立初始模型；

S2.如果来了新数据，则对新来的污水数据x_new用公式(2)进行预测并计算偏差；否则，转到S8；

S3.添加新样本(x_new,y_new)到模型中并初始化权重w为0，α_new初始化为1；

S5.估计噪声方差其中N为数据个数，M为基函数个数；

S6.重新计算∑，μ以及相应迭代过程中的S_m和Q_m，计算公式分别为∑＝(σ^-2Φ^TΦ+A)^-1，μ＝σ^-2∑Φ^Tt，

S7.若收敛或者达到最大迭代次数，则保存权重值和偏差，转S2继续预测数据；否则转S4。

S8.程序结束。

在线模型实际运行预测结果如图1所示。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (5)

1.一种基于快速相关向量机的污水处理在线软测量方法，其特征在于，包含以下顺序的步骤：

A、通过快速边际似然算法估计超参数，得到模型的权重值和样本偏差值；

B、然后建立快速相关向量机在线预测模型，对模型参数寻优，实现了污水中BOD的精确快速测量。

2.根据权利要求1所述的基于快速相关向量机的污水处理在线软测量方法，其特征在于，具体包含以下步骤：

S1.剔除输入和输出的数据中的异常点，由于各输入变量量纲的不同，对其进行归一化处理，归一化到[0,1]区间中；

S2.给定污水数据集{(x_n,t_n),n＝1,2,...,N},x_n∈R^d,t_n∈R,N是样本数,为简便起见,只考虑标量目标函数,我们遵循标准的概率公式,假定:

t_n＝y(x_n；w)+ε_n                    (1)

其中,y(·)是非线性函数，ε_n是均值为0,方差为σ²的高斯噪声，即因此有t_n～N(y(x_n),σ²)，这表明t_n满足均值y(x_n)为方差为σ²的高斯噪声分布，与支持向量机类似，函数y(x)定义为

$y(x;w)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{i}K(x,x_i)+w_0---\left(2\right)$

其中，由公式φ_i(x)＝K(x,x_i)来确定基函数，其核被训练向量参数化，假定t_n是相互独立的，则整个训练集的似然函数可写为

$p\left(t\right|w,{\mathrm{\&sigma;}}^2)={\left(2\mathrm{\&pi;}{\mathrm{\&sigma;}}^2\right)}^{-N/2}\exp (-\frac{{\left|\right|t-\mathrm{\&Phi;w}\left|\right|}^2}{2{\mathrm{\&sigma;}}^2})---\left(3\right)$

式中t＝[t₁,t₂,...,t_N]^T,w＝[w₀,w₁,...,w_M]^T,Φ是一个N×(N+1)的设计矩阵，Φ＝[φ₁,φ₂,...,φ_M]是组非线性基函数，φ(x_n)＝[1,K(x_n,x₁),K(x_n,x₂),...,K(x_n,x_N)]^T；

由于在模型中有和训练样本差不多的参数个数，从(3)式中得到的w和σ²的最大似然估计值有可能导致模型过拟合；为了避免过度拟合，通常的做法是给参数强加一些限制条件；在这里,我们从贝叶斯概率框架出发,通过定义一个先验概率分布来限制参数w和σ²；

这里选择一个比较平滑的函数,定义w的先验概率分布为零均值的高斯分布:

$p\left(w\right|\mathrm{\&alpha;})=\underset{j=0}{\overset{N}{\mathrm{\&Pi;}}}N\left(w_j\right|0,{\mathrm{\&alpha;}}_j^{-1})---\left(4\right)$

式(4)中:超参数α＝[α₀,α₁,...,α_N]^T，更重要的是，每个独立的超参数α_j只与其对应的权值w_j相关；通过这个限制条件，在经过大量的污水数据学习后,大部分超参数会趋近于无穷大,而与其对应的权值为0,从而使RVM具有较高的稀疏性；

现在已经定义了先验概率，从贝叶斯规则来看，对于给定的数据中未知数据，贝叶斯推理通过计算后验概率处理：

$p(w,\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2|t)=\frac{p\left(t\right|w,\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2)p(w,\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2)}{p\left(t\right)}---\left(5\right)$

给定一个测试点x_*,相应的污水出水水质预测值t_*的预测分布为

p(t_*|t)＝∫p(t_*|w,α,σ²)p(w,α,σ²|t)dwdαdσ²          (6)

根据贝叶斯公式,利用样本似然函数(4)和w先验分布(5)可得w的后验分布为

$p\left(w\right|\mathrm{t\&alpha;},\mathrm{\&beta;})=\frac{p\left(w\right|\mathrm{\&alpha;})p\left(t\right|w,\mathrm{\&beta;})}{p\left(t\right|\mathrm{\&alpha;},\mathrm{\&beta;})}---\left(7\right)$

我们把后验概率分解为

p(w,α,σ²|t)＝p(|w|t,α,σ²)p(α,σ²|t)           (8)

因此对权重的后验概率分布为

$p\left(w\right|t,\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2)=\frac{p\left(t\right|w,{\mathrm{\&sigma;}}^2)p\left(w\right|\mathrm{\&alpha;})}{p\left(t\right|\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2)}={\left(2\mathrm{\&pi;}\right)}^{-(N+1)/2}{\left|\mathrm{\&Sigma;}\right|}^{-1/2}\exp \{-\frac{1}{2}{(w-\mathrm{\&mu;})}^T{\mathrm{\&Sigma;}}^{-1}(w-\mathrm{\&mu;})\}---\left(9\right)$

其协方差为Σ＝(σ^-2Φ^TΦ+A)^-1              (10)

平均值为μ＝σ^-2ΣΦ^Tt              (11)

其中矩阵A＝diag(α₀,α₁,...,α_N)

${\left({\mathrm{\&alpha;}}_i\right)}^{\mathrm{new}}=\frac{{\mathrm{\&gamma;}}_i}{{\mathrm{\&mu;}}_i^2}---\left(12\right)$

${\left({\mathrm{\&sigma;}}^2\right)}^{\mathrm{new}}=\frac{{\left|\right|t-\mathrm{\&Phi;\&mu;}\left|\right|}^2}{N-{\mathrm{\&Sigma;}}_i{\mathrm{\&gamma;}}_i}---\left(13\right)$

其中γ_i≡1-α_iΣ_ii，Σ_ii为协方差矩阵Σ的第i个对角元素；最后通过公式(10)到公式(13)的迭代推理运算得到超参数α和方差σ²的估计值；输出的污水水质预测值为y_*＝μ^Tφ(x_*)，x_*是污水处理过程输入值；

S3.用快速边际似然算法估计超参数

针对相关向量机计算时间复杂度大、内存开销大的问题，采用了一种快速边际似然算法；它在训练过程中从空集开始动态地扩充基矩阵Φ，从而增大边际似然函数，或者去掉基矩阵Φ冗余的列来增大目标函数；

S4.对待预测的污水样本数据进行预测：将入水数据作为训练好的相关向量机软测量模型的输入，模型的输出即为出水BOD的预测结果。

3.根据权利要求2所述的基于快速相关向量机的污水处理在线软测量方法，其特征在于，所述的步骤S3，具体包含以下顺序的步骤：

相关向量机是通过最大化边际似然函数p(t|α,σ²)的方法确定超参数α和方差σ²的，等价于最大化为其对数；记L(α)＝log[p(t|α,σ²)]，整理有

$\begin{array}{c}L\left(\mathrm{\&alpha;}\right)=\log \left[p\left(t\right|\mathrm{\&alpha;},{\mathrm{\&sigma;}}^2)\right]=\log [\&Integral;p\left(t\right|w,{\mathrm{\&sigma;}}^2)p\left(w\right|\mathrm{\&alpha;})]\\ =-\frac{1}{2}[N\log 2\mathrm{\&pi;}+\log |C|+t^T{C}^{-1}t]\end{array}---\left(14\right)$

其中C＝σ²I+ΦA^-1Φ^T,t＝[t1,t2,…,t_N]^T；

为了便于最大化L(α)，对矩阵C进行等价变换，如下：

$\begin{array}{c}C={\mathrm{\&sigma;}}^{2}I+\mathrm{\&Phi;}{A}^{-1}{\mathrm{\&Phi;}}^T={\mathrm{\&sigma;}}^{2}I+\underset{m\&NotEqual;i}{\mathrm{\&Sigma;}}{\mathrm{\&alpha;}}_m^{-1}{\mathrm{\&phi;}}_m{\mathrm{\&phi;}}_m^T+{\mathrm{\&alpha;}}_i^{-1}{\mathrm{\&phi;}}_i{\mathrm{\&phi;}}_i^T\\ =C_{-i}+{\mathrm{\&alpha;}}_i^{-1}{\mathrm{\&phi;}}_i{\mathrm{\&phi;}}_i^T\end{array}---\left(15\right)$

其中此矩阵表示当α_i＝∞时，相应的基向量φ_i被移除后样本对应的协方差矩阵，根据矩阵相关性质整理可得

$\left|C\right|=\left|C_i\right||1+{\mathrm{\&alpha;}}_i^{-1}{\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}{\mathrm{\&phi;}}_i|$

$C^{-1}=C_{-i}^{-1}-\frac{C_{-i}^{-1}{\mathrm{\&phi;}}_i{\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}}{{\mathrm{\&alpha;}}_i+{\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}{\mathrm{\&phi;}}_i}---\left(16\right)$

因此公式(14)可以改写为

$\begin{array}{c}L\left(\mathrm{\&alpha;}\right)=-\frac{1}{2}[N\log 2\mathrm{\&pi;}+\log |C|+t^T{C}_{-i}t-\log {\mathrm{\&alpha;}}_i+\log ({\mathrm{\&alpha;}}_i+{\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}{\mathrm{\&phi;}}_i)-\frac{{\left({\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}t\right)}^2}{{\mathrm{\&alpha;}}_i+{\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}{\mathrm{\&phi;}}_i}]\\ =L\left({\mathrm{\&alpha;}}_{-i}\right)+\frac{1}{2}[\log {\mathrm{\&alpha;}}_i-\log ({\mathrm{\&alpha;}}_i+s_i)+\frac{{\left(q_i\right)}^2}{{\mathrm{\&alpha;}}_i+s_i}]\\ =L\left({\mathrm{\&alpha;}}_{-i}\right)+l\left({\mathrm{\&alpha;}}_i\right)\end{array}---\left(17\right)$

注意L(α_-i)表示为当α_i＝∞时，相应的基本向量φ_i被移除后所对应的边界似然函数的对数，而l(α_i)表示边界似然的对数函数中只与α_i有关的独立部分，上式对α_i求偏导有

$\frac{\&PartialD;L\left(\mathrm{\&alpha;}\right)}{\&PartialD;{\mathrm{\&alpha;}}_i}=\frac{\&PartialD;l\left({\mathrm{\&alpha;}}_i\right)}{\&PartialD;{\mathrm{\&alpha;}}_i}=\frac{1}{2}[\frac{1}{{\mathrm{\&alpha;}}_i}-\frac{1}{{\mathrm{\&alpha;}}_i+{\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}{\mathrm{\&phi;}}_i}-\frac{{\left({\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}t\right)}^2}{{\mathrm{\&alpha;}}_i+{\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}{\mathrm{\&phi;}}_i}]---\left(18\right)$

记 $S_i={\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}{\mathrm{\&phi;}}_i,Q_i={\mathrm{\&phi;}}_i^T{C}_{-i}^{-1}t---\left(19\right)$

所以公式(18)可改写为

$\frac{\&PartialD;L\left(\mathrm{\&alpha;}\right)}{\&PartialD;{\mathrm{\&alpha;}}_i}=\frac{{\mathrm{\&alpha;}}_i^{-1}{S}_i^2-(Q_i^2-S_i)}{2{({\mathrm{\&alpha;}}_i+S_i)}^2}---\left(20\right)$

令公式(20)等于零，考虑到α_i是方差值必须为正，所以当时有

${\mathrm{\&alpha;}}_i=\frac{S_i^2}{Q_i^2-S_i}---\left(21\right)$

对L(α)关于α_i求二阶偏导有

$\begin{array}{c}\frac{{\&PartialD;}^{2}L\left(\mathrm{\&alpha;}\right)}{\&PartialD;{\mathrm{\&alpha;}}_i^2}=\frac{-{\mathrm{\&alpha;}}_i^{-2}{S}_i^2{({\mathrm{\&alpha;}}_i+S_i)}^2-2({\mathrm{\&alpha;}}_i+S_i)[{\mathrm{\&alpha;}}_i^{-1}{S}_i^2-(Q_i^2-S_i)]}{2{({\mathrm{\&alpha;}}_i+S_i)}^4}\\ =\frac{{\mathrm{\&alpha;}}_i^{-2}{S}_i^2}{2{({\mathrm{\&alpha;}}_i+S_i)}^2}-\frac{S_i^2}{{\mathrm{\&alpha;}}_i{({\mathrm{\&alpha;}}_i+S_i)}^3}-\frac{[{\mathrm{\&alpha;}}_i^{-1}{S}_i^2-(Q_i^2-S_i)]}{{({\mathrm{\&alpha;}}_i+S_i)}^3}\end{array}---\left(22\right)$

综合公式(20)和(21)进行分析可知

$\frac{{\&PartialD;}^{2}L\left(\mathrm{\&alpha;}\right)}{\&PartialD;{\mathrm{\&alpha;}}_i^2}{|}_{{\mathrm{\&alpha;}}_i=\frac{S_i^2}{Q_i^2-S_i}}=\frac{-S_i^2}{2{\mathrm{\&alpha;}}_i^2{({\mathrm{\&alpha;}}_i+S_i)}^2}---\left(23\right)$

所以当时，公式(23)左边的表达式是恒小于零的，并对以上推导公式分析可得，L(α)有唯一最大值点为

${\mathrm{\&alpha;}}_i=\left\{\begin{array}{cc}\frac{S_i^2}{Q_i^2-S_i}& Q_i^2>S_i\\ \&infin;& Q_i^2\&le;S_i\end{array}\right.---\left(24\right)$

由此得到L(α)的最大值。

4.根据权利要求3所述的基于快速相关向量机的污水处理在线软测量方法，其特征在于，所述的贝叶斯L(α)通过以下方法最大化：

a、当基向量φ_i在模型中，即α_i<∞，但有则将φ_i从模型中删除，即令α_i＝∞，这样可以增大贝叶斯L(α)；

b、当基向量φ_i在模型中，即α_i＝∞，但有则将φ_i添加到模型中并利用公式(24)更新α_i，这样可以增大贝叶斯L(α)；

c、当基向量φ_i在模型中，即α_i<∞，但有则用公式(24)更新α_i，这样可以增大贝叶斯L(α)。

5.根据权利要求3所述的基于快速相关向量机的污水处理在线软测量方法，所述的快速相关向量机，其回归基本算法步骤如下：

Ⅰ、初始化σ²；

Ⅱ、用单个基向量φ_i初始化α_i，由公式(24)分析整理可得并设置其他的α_m(m≠i)为无穷大；

Ⅲ、计算Σ、μ并对所有M个基函数φ_m初始化S_m和Q_m；

Ⅳ、从所有M个基函数φ_m集合中选择候选的基向量φ_i；

Ⅴ、计算 ${\mathrm{\&theta;}}_i=Q_i^2-S_i;$

Ⅵ、若θ_i>0且α_i<∞(基向量φ_i在模型中)，重新估计α_i；

Ⅶ、若θ_i>0且α_i＝∞(基向量φ_i不在模型中)，添加φ_i到模型中并重新估计α_i；

Ⅷ、若θ_i≤0且α_i<∞，删除φ_i并设置α_i＝∞；

Ⅸ、估计噪声方差其中N为数据个数，M为基函数个数；

Ⅹ、重新计算协方差矩阵Σ，权重矩阵μ以及相应迭代过程中的S_m和Q_m；

Ⅺ、若收敛或者达到最大迭代次数，则保存权重值及偏差值，此次训练结束；否则转步骤IV继续训练。