CN104914227B - 基于多高斯核自优化相关向量机的污水水质软测量方法 - Google Patents

Abstract

本发明公开了一种基于多高斯核自优化相关向量机的污水水质软测量方法，包括以下步骤：1)剔除污水输入和输出的数据中的异常点，由于各输入变量量纲的不同，对其进行归一化处理，归一化到[0,1]区间中；2)多高斯核函数相关向量机软测量模型模块；3)多高斯核函数核参数自优化算法；4)遗传优化算法对初始参数寻优模块；5)多高斯核函数自优化相关向量机软测量模型建模。本发明通过自优化方法确定各尺度上的核参数，运用遗传优化算法对初始参数寻优，建立最优模型，在保证模型收敛性和稀疏性的情况下，有效提高污水中BOD输出精度。

Description

基于多高斯核自优化相关向量机的污水水质软测量方法

技术领域

本发明涉及软测量的技术领域，尤其是指一种基于多高斯核自优化相关向量机的污水水质软测量方法。

背景技术

随着全球经济的增长和人们生活水平的提高，城市生活污水和工业污水的排放量日益增大，有限水资源受到不同程度的污染，防止水污染保护人们赖以生存的水环境是人类面临的当务之急。保护水资源的一个重要方面就是要提高污水处理的技术水平和检测手段。污水处理过程相当复杂，具有参数时变，多变量耦合、强非线性，严重滞后等特点。使得出水水质不易在线测量，废水处理成本高。污水处理过程中一些自动化检测设备、仪表的功能还不完善，处理周期太长，远远达不到国家对环境保护的要求。同时国家对环境保护的投入加大，污水处理技术越来越受到更多的关注。国家发展规划中明确提出要研发并推广低能耗、有效的污水处理技术。解决污水生物处理过程中参数在线测量的方法有两种：一是改进测量仪表，二是深入研究软测量技术，这对于出水的优化排放有着重大的实际意义和应用价值。由于在线仪表设备投资大，维护困难以及分析周期长、准确性不高，而软仪表成本低廉，测量具有实时性，尽管IAWQ推出ASMs和ADM1等机理模型，但实际的污水处理现场条件太多，过程太复杂，很难建立精确的机理模型，而神经网络、支持向量机等非机理建模都存在局限性。相关向量机具有解更稀疏、核函数选择更自由，泛化能力更强，鲁棒性更好等优点，使其在污水领域的应用中越来越受重视。但是传统高斯核函数核参数均匀统一，限制了rvm模型的灵活性和预测精度。

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

现在的污水处理一般都采用稀释法、传感器等测量污水中BOD、COD的浓度，但由于分析测定这两个指标的周期较长，测量中时常出现误差，不能及时反应污水处理的现场情况，因而污水控制系统存在着较大的延时，不能发挥其最佳的性能。

软测量技术就是利用易测量过程变量，依据这些易测量的过程变量和难以直接测量的待测过程变量之间的数学模型，通过各种数学计算和估计方法，用计算机软件实现待测量过程变量的测量。是近年在过程控制和检测领域涌现的一种新技术。软测量作为现代复杂过程工业中较难甚至无法由硬件在线检测参量实时估计的有效手段，具有投资低、维护保养简单等优点。将软测量技术用于污水处理过程，能降低污水处理厂的能耗，节约成本，避免污水处理过程中严重滞后等不足。但是污水过程非常复杂，传统的机理建模无法适应，而神经网络、支持向量机等非机理建模都存在局限性，本发明提出多高斯核函数自优化相关向量机模型，打破了传统高斯核函数参数均匀统一对预测精度的限制，使模型更灵活。在保证模型收敛性和稀疏性的情况下，有效提高污水中BOD输出精度。该模型是对基本rvm的改进，不仅离线情况下获得良好效果，在此基础上展开基于多尺度高斯核自优化相关向量机的污水出水在线软测量研究也很有意义。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于多高斯核自优化相关向量机的污水水质软测量方法，通过自优化方法确定各尺度上的核参数，运用遗传优化算法对初始参数寻优，建立最优模型。在保证模型收敛性和稀疏性的情况下，有效提高污水中BOD输出精度。

为实现上述目的，本发明所提供的技术方案为：基于多高斯核自优化相关向量机的污水水质软测量方法，包括以下步骤：

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

2)多高斯核函数相关向量机软测量模型模块，其计算公式：

给定一组输入和目标值x_n∈R^M，且考虑目标函数只是一个标量，根据概率方程式，假设目标函数是模型的样本并且带有附加噪声：

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

式中ε_n为附加噪声，服从(0,σ²)的高斯分布，且彼此间相互独立，因此p(t_n|y(x_n),σ²)服从高斯分布，其分布由期望y(x_n)和方差σ²决定，其中y(x_n)可由核函数的加权模型表示：

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

w_i为模型权值；k是核函数，φ_i＝k(x,x_i)；t＝[t₁,t₂,...t_n]^T为N维列向量，假设t_n服从独立分布，数据集的似然估计概率为：

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

φ(x_n)＝[1,k(x_n,x₁),k(x_n,x₂),...k(x_n,x_n)]^T为(N+1)*1维矩阵，贝叶斯矩阵表示为：Φ＝[φ(x₁),φ(x₂),...,φ(x_n)]^T，其是N*(N+1)维；为了约束rvm模型中权值w，假设其服从(0,α_i ^-1)高斯分布，且设各权值间相互独立，其先验概率表示如下：

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

α为N+1维超参数，引入超参数导致算法的稀疏性；

根据贝叶斯准则，可得到后验概率公式：

给定新的测试样本x，预测相应的目标t_*，按照预测分布：

$p\left(t_{*}\right|t)=\∫p\left(t_{*}\right|w,\mathrm{\α},{\mathrm{\σ}}^2)p(w,\mathrm{\α},{\mathrm{\σ}}^2|t)d_w{d}_{\mathrm{\α}}{d}_{{\mathrm{\σ}}^2}$

无法计算后验概率，因为无法计算

所以将后验概率分解：p(w,α,σ²|t)＝p(w|t,α,σ²)p(α,σ²|t)

可得权重的后验概率为：

$\begin{array}{c}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}|\Σ{|}^{-1/2}\exp \{-\frac{1}{2}{(w-u)}^T{\Σ}^{-1}(w-u)\}\end{array}---\left(5\right)$

其后验协方差和均值分别为：

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

u＝σ^-2∑Φ^Tt (7)

因为p(α,σ²|t)∝p(t|α,σ²)p(α)p(σ²)，所以可以用后者等价前者：

$\begin{array}{c}p\left(t\right|a,{\mathrm{\σ}}^2)=\∫p\left(t\right|w,\mathrm{\β})p\left(w\right|\mathrm{\α})d_w\\ ={\left(2\mathrm{\π}\right)}^{-N/2}|{\mathrm{\σ}}^{2}I+{\mathrm{\ΦA}}^{-1}{\mathrm{\Φ}}^T{|}^{-1/2}\exp (-\frac{1}{2}{t}^{T}Ct)\end{array}$

其中矩阵C＝(σ²I+ΦA^-1Φ^T)^-1，对上述式子取对数得：

$L=\log p\left(t\right|\mathrm{\α},\mathrm{\β})=-\frac{1}{2}\left(Nlog\right(2\mathrm{\π})+log|C|+t^T{C}^{-1}t)---\left(8\right)$

为了找到超参数的优化值，我们对上述似然估计最大化得到超参数的更新公式：

$\begin{array}{cc}{a}_i=\frac{{\mathrm{\γ}}_i}{u_i}& {\mathrm{\σ}}^2=\frac{\left|\right|t-\mathrm{\Φ}u|{|}^2}{N-{\mathrm{\Σ}}_{i=1}^N{\mathrm{\γ}}_i}\end{array}---\left(9\right)$

式中γ_i＝1-α_i∑_ii

相较于线性核函数，多项式核函数，sigmoid核函数，高斯核函数使用得更多；传统高斯核函数核参数均匀统一，限制了rvm模型的灵活性和预测精度；基于提高预测精度的考虑，且不破坏其它性能，提出使用如下多高斯核函数：

$k(x_m,x_n)=\exp (-\underset{k=1}{\overset{d}{\Σ}}{\mathrm{\η}}_k{\left(x_{mk}-x_{nk}\right)}^2)$

其各个尺度上使用不同核参数，η_k是第k个核参数的平方的倒数，_d是输入向量的属性个数，则多尺度核参数可表示为：η＝(η₁,η₂,...η_d)；用φ_nm＝k(x_m,x_n)来表示贝叶斯矩阵的元素，则似然估计对第k个核参数的梯度为：

$\frac{{\∂}_L}{{\∂}_{{\mathrm{\η}}_k}}=\underset{n=1}{\overset{N}{\Σ}}\underset{m=1}{\overset{N}{\Σ}}\frac{{\∂}_L}{{\∂}_{{\mathrm{\φ}}_{nm}}}\frac{{\∂}_{{\mathrm{\φ}}_{nm}}}{{\∂}_{{\mathrm{\η}}_k}}$

该式的第一部分与贝叶斯核参数无关，用来表示，得到矩阵：

D＝(C^-1tt^TC^-1-C^-1)ΦA^-1

＝β[(t-y)u^T-Φ∑] (10)

似然估计对核参数求导为：

3)多高斯核函数核参数自优化算法

为了避免核参数过大或过小引起过平滑和过适应问题，增加限制条件从而避免出现上述问题：

$\underset{k=1}{\overset{d}{\Σ}}{\mathrm{log\η}}_k=d\log \frac{1}{h^2}$

h是核参数几何平均值，则似然估计模型变成：

该模型对核参数求导为：

$\frac{\∂L^{-}}{\∂{\mathrm{\η}}_k}=\frac{\∂L}{\∂{\mathrm{\η}}_k}+\frac{\mathrm{\λ}}{{\mathrm{\η}}_k},k=1,...,d$

$\frac{\∂L^{-}}{\∂{\mathrm{\η}}_k}=\underset{m=1}{\overset{N}{\Σ}}\underset{n=1}{\overset{N}{\Σ}}-D_{nm}{\mathrm{\Φ}}_{nm}{(x_{mk}-x_{nk})}^2+\frac{\mathrm{\λ}}{{\mathrm{\η}}_k},k=1,...,d$

当上式第一部分对核参数求导时，所得梯度同号，为了保证所得核参数均大于零且最大似然估计尽可能的大，可用如下公式更新核参数，

若不满足同号，为了保证核参数均大于零，且最大似然估计尽可能的大，梯度小于零的核参数其更新公式为：

η_inew＝p kη_i (13)

梯度大于零的核参数其更新公式为：

式中i表示梯度小于零的属性列，j表示梯度大于零的属性列，p表示大于零的梯度个数；

4)遗传优化算法对初始参数寻优模块

在建模过程中，多高斯核函数自优化相关向量机模型有两类参数，一类属于核函数参数，核参数的个数由输入属性个数决定，核参数可通过上述算法中介绍的自优化方法来确定；另一类是相关向量机模型的初始参数，多高斯核函数自优化相关向量机的初始参数包括核参数几何平均值h和pk，其中核参数几何平均值决定初始核参数的值，pk影响核参数的自优化过程；提出采用遗传优化算法，对多高斯核函数自优化相关向量机模型中的初始参数h,pk进行优化，具体优化算法的流程如下：

4.1)设定种群规模，遗传进化代数，染色体采用二进制编码，每个算子由两段编码组成，两段编码分别代表两个参数h,pk，并确定各参数取值范围，生成初始种群；

4.2)用初始种群训练多高斯核函数自优化相关向量机模型，目标函数取测试样本的均方根误差函数RMSE；

4.3)个体进行选择、重组、变异，终止条件若达到最大进化数，则停止进化，输出优化参数，否则继续寻优；

5)多高斯核函数自优化相关向量机软测量模型建模

污水处理中与BOD出水水质息息相关的输入属性为可降解固体浓度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；输入x_n∈R^M，是第n个输入值，t_n是BOD的输出值，则出水BOD的预测算法步骤如下：

5.1)对上述输入属性的污水数据进行归一化，去误差平滑处理，然后确定训练数据和测试数据；

5.2)根据训练样本中的污水数据建立初始化模型，初始化模型参数:α＝[1/N²,1/N²,...,1/N²]^T，σ²＝0.1var(t)，采用上面提到的遗传优化算法得到最佳初始核参数几何平均值h＝24.3372,pk＝0.2379,初始核参数η＝[1/h²,1/h²,...,1/h²]^T；

5.3)循环迭代

5.3.1)按上式(6),(7)更新u,∑；

5.3.2)按上式(9)更新超参数α,β；

5.3.3)如果超参数α_i大于给定的大值，则认为该超参数趋于无穷，从而将相应的权值行设为零，并且忽略相应的贝叶斯矩阵列；否则，训练数据的相应训练样本则为相关向量；

5.3.4)如果最大似然估计值变大，连续两代间的核参数最大变化值大于一定的小值，则按上面公式(12)，(13)，(14)更新核参数η，并更新贝叶斯矩阵Φ，否则停止更新；

5.3.5)如果连续两代间的超参数α最大变化值小于给定的小值，则迭代结束；

5.4)输出模型权值w，核参数η，相关向量机，噪声方差σ²。

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

1、本发明建立了一种基于多高斯核函数自优化相关向量机的软测量模型，模型通过自优化方法确定各尺度上的核参数，同时运用遗传优化算法对初始参数寻优，建立最优模型，在保证模型收敛性和稀疏性的情况下，有效提高污水中BOD输出精度，该模型是对基本rvm的改进，不仅离线情况下获得良好效果，在此基础上展开基于多尺度高斯核自优化相关向量机的污水出水在线软测量研究也很有意义。

2、本发明的核函数由多个高斯函数组成，具有多尺度核参数，打破了传统高斯核函数参数均匀统一对预测精度的限制，使模型更灵活。

附图说明

图1为本发明模型对出水BOD结果拟合图。

图2为本发明模型与基本rvm模型预测误差对比图。

图3为本发明模型与基本rvm模型收敛性对比图。

具体实施方式

下面结合具体实施例对本发明作进一步说明。

本实施例所述的基于多高斯核自优化相关向量机的污水水质软测量方法，通过自优化方法确定各尺度上的核参数，运用遗传优化算法对初始参数寻优，建立最优模型，在保证模型收敛性和稀疏性的情况下，有效提高污水中BOD输出精度。其具体情况如下：

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

2)多高斯核函数相关向量机软测量模型模块，其计算公式：

给定一组输入和目标值x_n∈R^M，且考虑目标函数只是一个标量，根据概率方程式，假设目标函数是模型的样本并且带有附加噪声：

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

式中ε_n为附加噪声，服从(0,σ²)的高斯分布，且彼此间相互独立，因此p(t_n|y(x_n),σ²)服从高斯分布，其分布由期望y(x_n)和方差σ²决定，其中y(x_n)可由核函数的加权模型表示：

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

w_i为模型权值；k是核函数，φ_i＝k(x,x_i)；t＝[t₁,t₂,...t_n]^T为N维列向量，假设t_n服从独立分布，数据集的似然估计概率为：

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

φ(x_n)＝[1,k(x_n,x₁),k(x_n,x₂),...k(x_n,x_n)]^T为(N+1)*1维矩阵，贝叶斯矩阵表示为：Φ＝[φ(x₁),φ(x₂),...,φ(x_n)]^T，其是N*(N+1)维；为了约束rvm模型中权值w，假设其服从(0,α_i ^-1)高斯分布，且设各权值间相互独立，其先验概率表示如下：

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

α为N+1维超参数，引入超参数导致算法的稀疏性；

根据贝叶斯准则，可得到后验概率公式：

给定新的测试样本x，预测相应的目标t_*，按照预测分布：

$p\left(t_{*}\right|t)=\∫p\left(t_{*}\right|w,\mathrm{\α},{\mathrm{\σ}}^2)p(w,\mathrm{\α},{\mathrm{\σ}}^2|t)d_w{d}_{\mathrm{\α}}{d}_{{\mathrm{\σ}}^2}$

无法计算后验概率，因为无法计算

所以将后验概率分解：p(w,α,σ²|t)＝p(w|t,α,σ²)p(α,σ²|t)

可得权重的后验概率为：

$\begin{array}{c}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}|\Σ{|}^{-1/2}\exp \{-\frac{1}{2}{(w-u)}^T{\Σ}^{-1}(w-u)\}\end{array}---\left(5\right)$

其后验协方差和均值分别为：

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

u＝σ^-2∑Φ^Tt (7)

因为p(α,σ²|t)∝p(t|α,σ²)p(α)p(σ²)，所以可以用后者等价前者：

$\begin{array}{c}p\left(t\right|a,{\mathrm{\σ}}^2)=\∫p\left(t\right|w,\mathrm{\β})p\left(w\right|\mathrm{\α})d_w\\ ={\left(2\mathrm{\π}\right)}^{-N/2}|{\mathrm{\σ}}^{2}I+{\mathrm{\ΦA}}^{-1}{\mathrm{\Φ}}^T{|}^{-1/2}\exp (-\frac{1}{2}{t}^{T}Ct)\end{array}$

其中矩阵C＝(σ²I+ΦA^-1Φ^T)^-1，对上述式子取对数得：

$L=\log p\left(t\right|\mathrm{\α},\mathrm{\β})=-\frac{1}{2}\left(Nlog\right(2\mathrm{\π})+log|C|+t^T{C}^{-1}t)---\left(8\right)$

为了找到超参数的优化值，对上述似然估计最大化得到超参数的更新公式：

$\begin{array}{cc}{a}_i=\frac{{\mathrm{\γ}}_i}{u_i}& {\mathrm{\σ}}^2=\frac{\left|\right|t-\mathrm{\Φ}u|{|}^2}{N-{\mathrm{\Σ}}_{i=1}^N{\mathrm{\γ}}_i}\end{array}---\left(9\right)$

式中γ_i＝1-α_i∑_ii

相较于线性核函数，多项式核函数，sigmoid核函数，高斯核函数使用得更多；传统高斯核函数核参数均匀统一，限制了rvm模型的灵活性和预测精度；基于提高预测精度的考虑，且不破坏其它性能，本文提出使用如下多高斯核函数：

$k(x_m,x_n)=\exp (-\underset{k=1}{\overset{d}{\Σ}}{\mathrm{\η}}_k{\left(x_{mk}-x_{nk}\right)}^2)$

其各个尺度上使用不同核参数，η_k是第k个核参数的平方的倒数，_d是输入向量的属性个数，则多尺度核参数可表示为：η＝(η₁,η₂,...η_d)；用φ_nm＝k(x_m,x_n)来表示贝叶斯矩阵的元素，则似然估计对第k个核参数的梯度为：

$\frac{{\∂}_L}{{\∂}_{{\mathrm{\η}}_k}}=\underset{n=1}{\overset{N}{\Σ}}\underset{m=1}{\overset{N}{\Σ}}\frac{{\∂}_L}{{\∂}_{{\mathrm{\φ}}_{nm}}}\frac{{\∂}_{{\mathrm{\φ}}_{nm}}}{{\∂}_{{\mathrm{\η}}_k}}$

该式的第一部分与贝叶斯核参数无关，用来表示，得到矩阵：

D＝(C^-1tt^TC^-1-C^-1)ΦA^-1

＝β[(t-y)u^T-Φ∑] (10)

似然估计对核参数求导为： $\frac{{\∂}_L}{{\∂}_{{\mathrm{\η}}_k}}=\underset{m=1}{\overset{N}{\Σ}}\underset{n=1}{\overset{N}{\Σ}}-D_{nm}{\mathrm{\Φ}}_{nm}{(x_{mk}-x_{nk})}^2$

3)多高斯核函数核参数自优化算法

为了避免核参数过大或过小引起过平滑和过适应问题，增加限制条件从而避免出现上述问题：

$\underset{k=1}{\overset{d}{\Σ}}{\mathrm{log\η}}_k=d\log \frac{1}{h^2}$

h是核参数几何平均值，则似然估计模型变成：

该模型对核参数求导为：

$\frac{\∂L^{-}}{\∂{\mathrm{\η}}_k}=\frac{\∂L}{\∂{\mathrm{\η}}_k}+\frac{\mathrm{\λ}}{{\mathrm{\η}}_k},k=1,...,d$

$\frac{\∂L^{-}}{\∂{\mathrm{\η}}_k}=\underset{m=1}{\overset{N}{\Σ}}\underset{n=1}{\overset{N}{\Σ}}-D_{nm}{\mathrm{\Φ}}_{nm}{(x_{mk}-x_{nk})}^2+\frac{\mathrm{\λ}}{{\mathrm{\η}}_k},k=1,...,d$

当上式第一部分对核参数求导时，所得梯度同号，为了保证所得核参数均大于零且最大似然估计尽可能的大，可用如下公式更新核参数，

若不满足同号，为了保证核参数均大于零，且最大似然估计尽可能的大，梯度小于零的核参数其更新公式为：

η_inew＝p kη_i (13)

梯度大于零的核参数其更新公式为：

式中i表示梯度小于零的属性列，j表示梯度大于零的属性列，p表示大于零的梯度个数；

4)遗传优化算法对初始参数寻优模块

在建模过程中，多高斯核函数自优化相关向量机模型有两类参数，一类属于核函数参数，核参数的个数由输入属性个数决定，核参数可通过上述算法中介绍的自优化方法来确定；另一类是相关向量机模型的初始参数，多高斯核函数自优化相关向量机的初始参数包括核参数几何平均值h和pk，其中核参数几何平均值决定初始核参数的值，pk影响核参数的自优化过程；本文采用遗传优化算法，对多高斯核函数自优化相关向量机模型中的初始参数h,pk进行优化，具体优化算法的流程如下：

4.1)设定种群规模，遗传进化代数，染色体采用二进制编码，每个算子由两段编码组成，两段编码分别代表两个参数h,pk，并确定各参数取值范围，生成初始种群；

4.2)用初始种群训练多高斯核函数自优化相关向量机模型，目标函数取测试样本的均方根误差函数RMSE；

4.3)个体进行选择、重组、变异，终止条件若达到最大进化数，则停止进化，输出优化参数，否则继续寻优；

5)多高斯核函数自优化相关向量机软测量模型建模

污水处理中与BOD出水水质息息相关的输入属性为可降解固体浓度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；输入x_n∈R^M，是第n个输入值，t_n是BOD的输出值，则出水BOD的预测算法步骤如下：

5.1)对上述输入属性的污水数据进行归一化，去误差平滑处理，然后确定训练数据和测试数据；

5.2)根据训练样本中的污水数据建立初始化模型，初始化模型参数:α＝[1/N²,1/N²,...,1/N²]^T，σ²＝0.1var(t)，采用上面提到的遗传优化算法得到最佳初始核参数几何平均值h＝24.3372,pk＝0.2379,初始核参数η＝[1/h²,1/h²,...,1/h²]^T；

5.3)循环迭代

5.3.1)按上式(6),(7)更新u,∑；

5.3.2)按上式(9)更新超参数α,β；

5.3.3)如果超参数α_i大于给定的大值，则认为该超参数趋于无穷，从而将相应的权值行设为零，并且忽略相应的贝叶斯矩阵列；否则，训练数据的相应训练样本则为相关向量；

5.3.4)如果最大似然估计值变大，连续两代间的核参数最大变化值大于一定的小值，则按上面公式(12)，(13)，(14)更新核参数η，并更新贝叶斯矩阵Φ，否则停止更新；

5.3.5)如果连续两代间的超参数α最大变化值小于给定的小值，则迭代结束；

5.4)输出模型权值w，核参数η，相关向量机，噪声方差σ²。

下面我们以加州大学数据库(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组用作新数据测试模型精度。本发明模型对出水BOD结果拟合如图1所示，本发明模型与基本rvm模型预测误差对比如图2所示，本发明模型与基本rvm模型收敛性对比如图3所示。

算法步骤：

1.对上述输入属性的污水数据进行归一化，去误差平滑处理，然后确定训练数据和测试数据。

2.根据训练样本中的污水数据建立初始化模型，初始化模型参数:α＝[1/N²,1/N²,...,1/N²]^T，σ²＝0.1var(t)，采用上节提到的遗传优化算法得到最佳初始核参数几何平均值h＝24.3372,pk＝0.2379,初始核参数η＝[1/h²,1/h²,...,1/h²]^T。

3.循环迭代

1)按式(6),(7)更新u,∑。

2)按式(9)更新超参数α,β。

3)如果超参数αi大于给定的大值，则认为该超参数趋于无穷，从而将相应的权值行设为零，并且忽略相应的贝叶斯矩阵列。否则，训练数据的相应训练样本则为相关向量。

4)如果最大似然估计值变大，连续两代间的核参数最大变化值大于一定的小值，则按公式(12)，(13)，(14)更新核参数η，并更新贝叶斯矩阵Φ，否则停止更新。

5)如果连续两代间的超参数α最大变化值小于某一小值，则迭代结束。

4.输出模型权值w，核参数η，相关向量机，噪声方差σ²。

以上所述之实施例子只为本发明之较佳实施例，并非以此限制本发明的实施范围，故凡依本发明之形状、原理所作的变化，均应涵盖在本发明的保护范围内。

Claims (1)

1.基于多高斯核自优化相关向量机的污水水质软测量方法，其特征在于，包括以下步骤：

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

2)多高斯核函数相关向量机软测量模型模块，其计算公式：

给定一组输入和目标值x_n∈R^M，且考虑目标函数只是一个标量，根据概率方程式，假设目标函数是模型的样本并且带有附加噪声：

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

式中ε_n为附加噪声，服从(0,σ²)的高斯分布，且彼此间相互独立，因此p(t_n|y(x_n),σ²)服从高斯分布，其分布由期望y(x_n)和方差σ²决定，其中y(x_n)可由核函数的加权模型表示：

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

w_i为模型权值；k是核函数，φ_i＝k(x,x_i)；t＝[t₁,t₂,...t_n]^T为N维列向量，假设t_n服从独立分布，数据集的似然估计概率为：

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

φ(x_n)＝[1,k(x_n,x₁),k(x_n,x₂),...k(x_n,x_n)]^T为(N+1)*1维矩阵，贝叶斯矩阵表示为：Φ＝[φ(x₁),φ(x₂),...,φ(x_n)]^T，其是N*(N+1)维；为了约束rvm模型中权值w，假设其服从(0,α_i ^-1)高斯分布，且设各权值间相互独立，其先验概率表示如下：

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

α为N+1维超参数，引入超参数导致算法的稀疏性；

根据贝叶斯准则，可得到后验概率公式：

给定新的测试样本x，预测相应的目标t_*，按照预测分布：

$p\left(t_{*}\right|t)=\∫p\left(t_{*}\right|w,\mathrm{\α},{\mathrm{\σ}}^2)p(w,\mathrm{\α},{\mathrm{\σ}}^2|t)d_w{d}_{\mathrm{\α}}{d}_{{\mathrm{\σ}}^2}$

无法计算后验概率，因为无法计算

所以将后验概率分解：p(w,α,σ²|t)＝p(w|t,α,σ²)p(α,σ²|t)

可得权重的后验概率为：

$\begin{array}{c}p\left(w|t,\mathrm{\α},{\mathrm{\σ}}^2\right)=\frac{p\left(t|w,{\mathrm{\σ}}^2\right)p\left(w|\mathrm{\α}\right)}{p\left(t|\mathrm{\α},{\mathrm{\σ}}^2\right)}\\ ={\left(2\mathrm{\π}\right)}^{-\left(N+1\right)/2}|\Σ{|}^{-1/2}\exp \left\{-\frac{1}{2}{\left(w-u\right)}^T{\Σ}^{-1}\left(w-u\right)\right\}\end{array}---\left(5\right)$

其后验协方差和均值分别为：

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

u＝σ^-2∑Φ^Tt (7)

因为p(α,σ²|t)∝p(t|α,σ²)p(α)p(σ²)，所以可以用后者等价前者：

$\begin{array}{c}p\left(t\right|a,{\mathrm{\σ}}^2)=\∫p\left(t|w,\mathrm{\β}\right)p\left(w|\mathrm{\α}\right)d_w\\ ={\left(2\mathrm{\π}\right)}^{-N/2}|{\mathrm{\σ}}^{2}I+{\mathrm{\ΦA}}^{-1}{\mathrm{\Φ}}^T{|}^{-1/2}\exp \left(-\frac{1}{2}{t}^{T}Ct\right)\end{array}$

其中矩阵C＝(σ²I+ΦA^-1Φ^T)^-1，对上述式子取对数得：

$L=\log p\left(t\right|\mathrm{\α},\mathrm{\β})=-\frac{1}{2}\left(N\log \right(2\mathrm{\π})+log|C|+t^T{C}^{-1}t)---\left(8\right)$

为了找到超参数的优化值，对上述似然估计最大化得到超参数的更新公式：

$\begin{array}{cc}{a}_i=\frac{{\mathrm{\γ}}_i}{u_i}& {\mathrm{\σ}}^2=\frac{\left|\right|t-\mathrm{\Φ}u|{|}^2}{N-{\mathrm{\Σ}}_{i=1}^N{\mathrm{\γ}}_i}\end{array}---\left(9\right)$

式中γ_i＝1-α_i∑_ii

传统高斯核函数核参数均匀统一，限制了rvm模型的灵活性和预测精度；

基于提高预测精度的考虑，且不破坏其它性能，提出使用如下多高斯核函数：

$k(x_m,x_n)=\exp (-\underset{k=1}{\overset{d}{\Σ}}{\mathrm{\η}}_k{\left(x_{mk}-x_{nk}\right)}^2)$

其各个尺度上使用不同核参数，η_k是第k个核参数的平方的倒数，d是输入向量的属性个数，则多尺度核参数可表示为：η＝(η₁,η₂,...η_d)；用φ_nm＝k(x_m,x_n)来表示贝叶斯矩阵的元素，则似然估计对第k个核参数的梯度为：

$\frac{{\∂}_L}{{\∂}_{{\mathrm{\η}}_k}}=\underset{n=1}{\overset{N}{\Σ}}\underset{m=1}{\overset{N}{\Σ}}\frac{{\∂}_L}{{\∂}_{{\mathrm{\φ}}_{nm}}}\frac{{\∂}_{{\mathrm{\φ}}_{nm}}}{{\∂}_{{\mathrm{\η}}_k}}$

该式的第一部分与贝叶斯核参数无关，用来表示，得到矩阵：

D＝(C^-1tt^TC^-1-C^-1)ΦA^-1

＝β[(t-y)u^T-Φ∑] (10)

似然估计对核参数求导为：

3)多高斯核函数核参数自优化算法

为了避免核参数过大或过小引起过平滑和过适应问题，增加限制条件从而避免出现上述问题：

$\underset{k=1}{\overset{d}{\Σ}}{\mathrm{log\η}}_k=d\log \frac{1}{h^2}$

h是核参数几何平均值，则似然估计模型变成：

该模型对核参数求导为：

$\frac{\∂L^{-}}{\∂{\mathrm{\η}}_k}=\frac{\∂L}{\∂{\mathrm{\η}}_k}+\frac{\mathrm{\λ}}{{\mathrm{\η}}_k},k=1,\mathrm{...},d$

$\frac{\∂L^{-}}{\∂{\mathrm{\η}}_k}=\underset{m=1}{\overset{N}{\Σ}}\underset{n=1}{\overset{N}{\Σ}}-D_{nm}{\mathrm{\Φ}}_{nm}{(x_{mk}-x_{nk})}^2+\frac{\mathrm{\λ}}{{\mathrm{\η}}_k},k=1,\mathrm{...},d$

当上式第一部分对核参数求导时，所得梯度同号，为了保证所得核参数均大于零且最大似然估计尽可能的大，可用如下公式更新核参数，

若不满足同号，为了保证核参数均大于零，且最大似然估计尽可能的大，梯度小于零的核参数其更新公式为：

η_inew＝p kη_i (13)

梯度大于零的核参数其更新公式为：

式中i表示梯度小于零的属性列，j表示梯度大于零的属性列，p表示大于零的梯度个数；

4)遗传优化算法对初始参数寻优模块

在建模过程中，多高斯核函数自优化相关向量机模型有两类参数，一类属于核函数参数，核参数的个数由输入属性个数决定，核参数可通过上述算法中介绍的自优化方法来确定；另一类是相关向量机模型的初始参数，多高斯核函数自优化相关向量机的初始参数包括核参数几何平均值h和pk，其中核参数几何平均值决定初始核参数的值，pk影响核参数的自优化过程；提出采用遗传优化算法，对多高斯核函数自优化相关向量机模型中的初始参数h,pk进行优化，具体优化算法的流程如下：

4.1)设定种群规模，遗传进化代数，染色体采用二进制编码，每个算子由两段编码组成，两段编码分别代表两个参数h,pk，并确定各参数取值范围，生成初始种群；

4.2)用初始种群训练多高斯核函数自优化相关向量机模型，目标函数取测试样本的均方根误差函数RMSE；

4.3)个体进行选择、重组、变异，终止条件若达到最大进化数，则停止进化，输出优化参数，否则继续寻优；

5)多高斯核函数自优化相关向量机软测量模型建模

污水处理中与BOD出水水质息息相关的输入属性为可降解固体浓度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；输入x_n∈R^M，是第n个输入值，t_n是BOD的输出值，则出水BOD的预测算法步骤如下：

5.1)对上述输入属性的污水数据进行归一化，去误差平滑处理，然后确定训练数据和测试数据；

5.2)根据训练样本中的污水数据建立初始化模型，初始化模型参数:α＝[1/N²,1/N²,...,1/N²]^T，σ²＝0.1var(t)，初始核参数η＝[1/h²,1/h²,...,1/h²]^T；

5.3)循环迭代

5.3.1)按上式(6),(7)更新u,∑；

5.3.2)按上式(9)更新超参数α,β；

5.3.3)如果超参数α_i大于给定的大值，则认为该超参数趋于无穷，从而将相应的权值行设为零，并且忽略相应的贝叶斯矩阵列；否则，训练数据的相应训练样本则为相关向量；

5.3.4)如果最大似然估计值变大，连续两代间的核参数最大变化值大于一定的小值，则按上面公式(12)，(13)，(14)更新核参数η，并更新贝叶斯矩阵Φ，否则停止更新；

5.3.5)如果连续两代间的超参数α最大变化值小于给定的小值，则迭代结束；

5.4)输出模型权值w，核参数η，相关向量机，噪声方差σ²。