一种基于灰色理论及支持向量机的水质预测方法

摘要

本发明公开了一种基于灰色理论及支持向量机的水质预测方法，属于水质预测方法技术领域。该方法包括以下步骤：1）对采集的数据集进行小波分解；2）将步骤一中经过分解处理得到的细节子序列进行平移，并进行相空间重构；3）根据步骤二中经过相空间重构的细节子序列建立各细节子序列的预测模型；4）根据步骤一中经过分解处理得到的尺度子序列建立尺度子序列预测模型；5）输入数据进行预测，将得到的细节子序列的预测值进行反向平移，并将尺度序列模型和细节序列模型得到的预测值进行重构，得到最终预测结果。本方法与传统的水质预测方法相比具有更高的预测精度及更广的应用范围。同时，在预测稳定性和精度上，有更好地鲁棒性和复杂非线性映射的能力。

说明书

技术领域

本发明涉及一种水质预测方法，特别涉及一种基于灰色理论及支持向量机的水质预测方 法。

背景技术

水质预测是水资源管理和水污染控制的重要方法，是进行水质治理和水资源开发利用的 基础性工作，对及时掌握水质变化的发展趋势，为水质预测预警和水环境污染防治提供重要 的科学依据和技术支持。水质预测方法较多，典型的有基于数理统计、混沌理论和神经网络 方法。数理统计法存在对水质多因素预测难度大的问题，混沌理论的应用需要大量的水质信 息资料，神经网络预测方法能够刻画出水质时间序列复杂变化的非线性特征，但其本身存在 网络收敛速度慢、容易陷入局部极小值等问题，使得模型的泛化能力不强。因此需要一种先 进的预测算法，能够进一步提高水质预测模型的精度和效率。

灰色理论模型能够利用较少的系统信息，来描述系统的整体行为。GM(1,1)模型是一种呈 指数增长的模型，其预测精度受到原始监测序列光滑程度的限制，对近似服从、或符合指数 增长态势的动态序列具有良好的预测效果。支持向量机作为建立在统计学习理论VC维的基础 上，同时实现SRM(Structural Risk Minimization)准则的一种新的机器学习方法，在解决小样本、 非线性问题、高维模式识别和局部极小等方面具有优势，目前已被广泛的用于时间序列的预 测。但目前使用灰色理论和支持向量机对水质时间序列的预测仍然存在一些问题：其一，支 持向量机对水质时间序列的拟合预测精度受其参数选择的影响，已有的基于单因素参数寻优 的方法，受各因素之间没有交互作用条件的限制，不仅非常耗时，而且常常不能保证得到的 参数是最优的；其二，水质监测参数的时间序列变化受诸多因素的影响，某些监测值表现出 突发跳变等不规律变化的情况，导致在模型训练过程中常常得不到理想的模型；其三，采用 小波分解的方法虽然能够获得平稳的子序列，减小瞬变数据对预测模型的影响，但是大量输 入的训练数据仍然会影响模型的预测精度。

发明内容

有鉴于此，本发明的目的在于提供一种基于灰色理论及支持向量机的水质预测方法，该 方法结合灰色理论模型、支持向量机及小波变换等对水质进行高精度的预测。

为达到上述目的，本发明提供如下技术方案：

一种基于灰色理论及支持向量机的水质预测方法，包括以下步骤：步骤一：对采集的数 据集进行小波分解；步骤二：将步骤一中经过分解处理得到的细节子序列进行平移，并进行 相空间重构；步骤三：根据步骤二中经过相空间重构的细节子序列建立各细节子序列的预测 模型；步骤四：根据步骤一中经过分解处理得到的尺度子序列建立尺度子序列预测模型；步 骤五：输入数据进行预测，将得到的细节子序列的预测值进行反向平移，并将尺度序列模型 和细节序列模型得到的预测值进行重构，得到最终预测结果。

进一步，在步骤一中，选取DbN小波函数进行小波分解，并根据最小拟合误差确定分解 级数；将得到的分解序列进行单支重构，进而得到重构子序列，包括细节子序列和尺度子序 列。

进一步，在步骤二中，将步骤一中得到的细节子序列分别平移到正数区域；再根据相空 间重构理论，将其进行相空间重构，转化成矩阵的形式。

进一步，在步骤三中，将经过相空间重构的细节子序列分别输入支持向量机模型，结合 遗传算法进行参数优化，建立各细节子序列的预测模型。

进一步，在步骤四中，采用窗口滑动的形式，将尺度子序列x_t-m,x_t-m+1,x_t-m+2,…,x_t-1输入 灰色理论模型训练，预测t时刻的值，其中m为嵌入维数，建立尺度子序列预测模型。

进一步，在步骤五中，将细节子序列进行反向平移，将尺度序列模型及细节序列模型得 到预测的值根据以下公式进行重构，最终得到预测结果；即：

初始时间序列S(t)可以由尺度子序列A_N和细节子序列D₁,D₂,…,D_N表示：

$S\left(t\right)=A_N+\underset{i=1}{\overset{N}{\Sigma }}{D}_i.$

本发明的有益效果在于：本发明提出的预测模型与传统的支持向量机预测模型和神经网 络预测模型相比，具有更高的预测精度及更广的应用范围。同时，在预测稳定性和精度上， 有更好地鲁棒性和复杂非线性映射的能力。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为本发明所述方法的流程图；

图2为预测模型建立过程示意图。

具体实施方式

本发明提出的水质预测模型，通过使用小波变换等数据处理方法，将水质时间序列变换 到多个尺度上，有效的减小了瞬变数据对建模的影响；其次，结合遗传算法对支持向量回归 机的参数进行寻优，比传统的基于多次单因素实验寻找最优参数的方法更高效；最后，用对 平稳序列具有良好预测效果的灰色模型，对分解后的尺度序列进行预测。通过少量的数据建 立灰色预测模型，避免模型受到历史数据变化趋势的影响，从而提高预测精度。

下面将结合附图，对本发明的优选实施例进行详细的描述。

图1为本发明所述水质预测方法的流程图，如图所示，本方法包括5个步骤，1）对采集 的数据集进行小波分解；2）将步骤一中经过分解处理得到的细节子序列进行平移，并进行相 空间重构；3）根据步骤二中经过相空间重构的细节子序列建立各细节子序列的预测模型；4） 根据步骤一中经过分解处理得到的尺度子序列建立尺度子序列预测模型；5）输入数据进行预 测，将得到的细节子序列的预测值进行反向平移，并将尺度序列模型和细节序列模型得到的 预测值进行重构，得到最终预测结果。

小波多分辨率分析的主要思想是用不同的分辨率来逐级逼近待分析序列。任意一个初始 时间序列为S(t),t＝1,2,…,n。通过高通和低通滤波器对S(t)的一步分解可以将其变换为低频 部分cA₁和高频部分cD₁。可以继续对低频部分cA₁进行分解，当进行尺度1,2,…,N步分解后， 就得到一个低频部分cA_N和N个高频部分cD₁,cD₂,…,cD_N。为使得分解序列与初始时间序列在 时间上具有一一对应的关系，根据小波多分辨分析的分解和重构特性，将低频部分cA_N和高 频部分cD₁,cD₂,…,cD_N进行单支重构，可以得到低频部分的尺度子序列A_N和高频部分的细节 子序列D₁,D₂,…,D_N。分解后细节系数D₁为随机分量，呈现无规律的波动；系数D₂,…,D_N数 量级较小，有较强的规律性，一定程度上表示原始曲线的峰值；尺度系数A_N波形平滑，数值 较大，表示原始序列的趋势走向。

初始时间序列S(t)可以由分解序列A_N,D₁,D₂,…,D_N表示：

$S\left(t\right)=A_N+\underset{i=1}{\overset{N}{\Sigma }}{D}_i.$

由小波变换对原始的时间序列S(t)进行多分辨分解，得到N阶分解下的尺度子序列A_N， 细节子序列D₁,D₂,…,D_N。

在各子序列分别进行模型训练前，细节子序列D₁,D₂,…,D_N中数据由于存在多位小数，且 存在数量级的差别，为减小小数部分的拟合预测误差、降低数量级差，将序列中的数据统一 平移到正数区域。若原始序列为D_i＝{d₁,d₂,d₃,…,d_n},i＝1,2,…,N，对序列做平移处理后的序 列为：

D_i⁺＝{d₁+Q,d₂+Q,d₃+Q,…,d_n+Q}Q∈N+,i＝1,2,…,N

式中平移因子Q的值为使得序列D_i⁺中的所有数据为正数的最小正整数。

细节子序列{x₁,x₂,x₃,…,x_n}进行模型训练前，需根据Takens理论进行相空间重构，转化 成矩阵的形式，获得数据之间的关联关系来得到尽可能多的信息量：

$Y=\left(\begin{array}{c}{x}_{m+1}\\ x_{m+2}\\ ·\\ ·\\ ·\\ x_n\end{array}\right)$

式中m为重构相空间的嵌入维数。重构相空间的嵌入维数m反映重构矩阵的信息量，不同的 m值对模型的预测效果有不同的影响，可以根据预测误差值的大小来优化选取相空间的嵌入 维数。

设为预测的目标值，将之前的目标值X_t＝{x_t-m,x_t-m+1,x_t-m+2,…,x_t-1}作为相关量， 建立自相关输入X_t与输出之间一一映射的关系：f:R^m→R，使得：

${\hat{Y}}_t=f\left(X_t\right)=f(x_{t-m},x_{t-m+1},x_{t-m+2},...,x_{t-1}).$

在本发明所述步骤中用到了灰色理论模型GM：

GM(1,1)模型可以用一个单变量的一阶微分方程表示：

$\frac{\mathrm{dx}}{\mathrm{dt}}+\mathrm{ax}=u$

已知拥有n个时间序列数据的原始序列为：x⁽⁰⁾＝(x⁽⁰⁾(1),x⁽⁰⁾(2),...,x⁽⁰⁾(n))

首先进行一阶累加生成新的序列：

x⁽¹⁾＝(x⁽¹⁾(1),x⁽¹⁾(2),...,x⁽¹⁾(n))

其中：$x^{\left(1\right)}\left(k\right)=\underset{i=1}{\overset{k}{\Sigma }}{x}^{\left(0\right)}\left(i\right),k=\mathrm{1,2},...,n.$根据一阶微分方程得到估计参数向量$\hat{\phi }=\left(\left(\begin{array}{c}a\\ u\end{array}\right)\right),$其中a是 发展灰数参数，u是控制灰数参数，可以得到：

$x^{\left(0\right)}(k+1)=a[-\frac{1}{2}(x^{\left(1\right)}\left(k\right)+x^{\left(1\right)}(k+1))]+u,k=\mathrm{1,2},...,n$

根据方程参数向量可用最小二乘法求取：

$\hat{\phi }={\left(B^{T}B\right)}^{-1}{B}^{T}Y$

其中$Y=\left(\left(\begin{array}{c}{x}^{\left(0\right)}\left(2\right)\\ x^{\left(0\right)}\left(3\right)\\ ·\\ ·\\ ·\\ x^{\left(0\right)}\left(n\right)\end{array}\right)\right),B=\left(\left(\begin{array}{cc}-\frac{1}{2}[x^{\left(1\right)}\left(1\right)+x^{\left(1\right)}\left(2\right)]& 1\\ -\frac{1}{2}[x^{\left(1\right)}\left(2\right)+x^{\left(1\right)}\left(3\right)]& 1\\ ·& ·\\ ·& ·\\ ·& ·\\ -\frac{1}{2}[x^{\left(1\right)}(n-1)+x^{\left(1\right)}\left(n\right)]& 1\end{array}\right)\right).$微分方程的解为：

${\hat{x}}^{\left(1\right)}(k+1)=[x^{\left(0\right)}\left(1\right)-\frac{u}{a}]e^{-\mathrm{ak}}+\frac{u}{a},k=\mathrm{1,2},...,n$

对进行逆累加生成还原，可得到预测值，即为GM(1，1)预测模型：

${\hat{x}}^{\left(0\right)}(k+1)={\hat{x}}^{\left(1\right)}(k+1)-{\hat{x}}^{\left(1\right)}\left(k\right),k=\mathrm{1,2},...,n.$

在本发明中还用到了支持向量回归算法：若给出一个训练样本集{(X_i,Y_i)|i＝1,2,…,l}， 其中X_i∈R^N为N维向量，Y_i∈R为相对应的输出变量。支持向量回归的基本思想是将非线性 的原始数据X_i,i＝1,2,…,l，映射到高维的特征空间H(Hibert空间)，即达到将低维空间的非 线性问题转化为高维空间的线性问题。当给出新的测试样本数据时，然后找出一个最优回归 超平面(OptimalHyper-plane)，能够通过回归超平面函数得到目标值。

其中ε-SVR支持向量机求解最优回归超平面的问题，为一个二次凸规划问题：

$\min \frac{1}{2}{\left|\right|\omega \left|\right|}^2+C\underset{i=1}{\overset{l}{\Sigma }}({\xi }_i+{\xi }_i^{*})$

上两式中ω,ω∈R^N为可调的权值函数；C为惩罚因子，为模型的复杂度与拟合精度之间的折 衷；ε为不敏感损失值，是拟合允许的误差；b为阈值；为X_i在高维特征 空间的映射值；ξ_i,是超出误差ε的溢出样本点偏差值。

由Lagrange优化方法将（1）式转化为对偶优化问题：

$s.t.\left(\begin{array}{c}\underset{i=1}{\overset{l}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})=0\\ {\alpha }_i,{\alpha }_i^{*}\in [0,C]\end{array}\right)$

式中α_i、为拉格朗日乘子。

高维空间中的向量点积可以由满足Mercer条件的对称核函数K(X_i,X_j)代 替，从而巧妙地降低了非线性映射的计算量。常用的核函数有多项式核函数、径向基核函数 和Sigmoid核函数等。由于径向基核函数在非线性系统识别方面更具有优势，本发明因此选 取径向基核函数：

$K(X_i,X_j)=\exp [-\frac{{\left|\right|X_i-X_j\left|\right|}^2}{2{\sigma }^2}]$

最后得到最优回归超平面的近似函数：

$f\left(X\right)=\underset{i=1}{\overset{l}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})K(X,X_i)+b$

并且再由KTT条件可以求得阈值：

$\left(\begin{array}{c}b=\frac{1}{N_{\mathrm{SV}}}\left\{\underset{0<a_i<C}{\Sigma }\right[Y_i-\underset{X_i,X_j\in \mathrm{SV}}{\Sigma }(\alpha -{\alpha }_i^{*})K(X_i,X_j)-ϵ]\\ +\underset{0<a_i^{*}<C}{\Sigma }[Y_i-\underset{X_i,X_j\in \mathrm{SV}}{\Sigma }({\alpha }_i-{\alpha }_i^{*})K(X_i,X_j)-ϵ]\}\end{array}\right)$

式中N_SV为支持向量总数，SV为支持向量集，说明最优回归超平面的构建是基于支持向量。

在本发明中，结合遗传算法对支持向量回归机的参数进行寻优，比传统的基于多次单因 素实验寻找最优参数的方法更高效。

遗传算法是以自然选择和遗传理论为基础，模拟了生物进化过程中适者生存规则与群体 内部染色体的随机信息交换机制，在模式定理和积木块假设保证能够找到全局最优解的条件 下，能自适应的控制搜索过程的高效全局搜索算法。本发明基于遗传算法的思想，建立适应 度函数，来解决支持向量机参数寻优这类多元单峰值的优化问题。

K-折交叉检验及遗传算法的适应度函数：

在K-折交叉检验中，首先将训练集分为K个数量相同的子集，然后取其中的一个子集作 为模型测试集，剩下的K-1个子集作为模型训练集。这样重复K次，使得训练集中的所有样 本都能够被模型预测一次。通过交叉验证分割训练集进行训练测试的方法，充分验证了训练 集中各数据的特性，使得预测结果更好地反映了模型实际效果，因此能够防止被选择的模型 出现过拟合的情况。目前在进行支持向量回归模型训练的时候，也采用交叉验证的方法。

通过遗传算法优化支持向量回归模型的参数时，需要比较种群个体的适应度大小，这里 适应度函数定义为负的均方误差（MSE）：

$f\left(x\right)=-\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{(Y_i-{\hat{Y}}_i)}^2$

式中N为训练样本的数目；Y_i为实际值；为使用交叉验证方法获得的预测值。

遗传算法参数优化步骤为：

惩罚因子用C表示、损失值用E表示、径向基核函数的参数用G表示。

步骤1：设置初始值：最大进化代数、种群规模和交叉验证折数K等。随机初始化一个 种群Init_Pop,进行模型训练，求得种群中各个体适应度值Valu_F；

步骤2：根据父代种群Init_Pop中个体，计算最高适应度，并选取适应度最高时的C、 E、G值，然后将它们分别作为当前最佳适应度值Best_Fit、最佳Best_C、Best_E、Best_G 值；

步骤3：选择算子：从种群Init_Pop中选出一部分适应度高的优良个体构成一个新的种 群；

步骤4：交叉算子：基于步骤3得到的种群，完成种群中的个体染色体的交叉重组，形 成一个新的种群；

步骤5：变异算子：基于步骤4得到的种群，完成种群中个体染色体的基因突变，而增 加群体多样性，形成一个新的种群，然后进行模型训练，求得到该种群中各个体适应度值 Valu_M；

步骤6：形成子代种群：比较种群Init_Pop和步骤5得到的种群个体适应度值Valu_F 和Valu_M，将步骤5得到种群中适应度值高的个体替代种群Init_Pop中适应度低的个体， 得到适应度水平更高的新种群Init_Pop。此时Init_Pop更新作为子代，并求的该子代种群 此时的最高适应度Max_Fit和取最高适应度的C、E、G值。

步骤7：比较步骤6中的最高适应度Max_Fit与当前最佳适应度值Best_Fit，若 Max_Fit大于Best_Fit，则更新最佳适应度值Best_Fit和最佳Best_C、Best_E、Best_G 值；否则不做任何操作，转向下一步继续执行。

步骤8：未达到最大进化代数，转向步骤3；否则算法停止，求得最优Best_C、Best_E、 Best_G值。

下面通过具体实施例对本发明所述预测方法的步骤进行详细说明，图2为预测模型建立 过程示意图，其过程如图所示：

步骤1在构建水质预测模型时，首先对数据集进行小波分解，首先选择合适的小波基 函数。研究表明，Daubechies小波基能较好地分析时间序列问题，故本发明选取DbN小波函 数。进行小波分解时，不同的分解级数会影响模型构建的消耗与预测的精度，需根据实验以 最小拟合误差为目标给出。若初始的水质时间序列为S(t),t＝1,2,…,n，首先将S(t)进行N阶 小波分解，得到一个低频序列cA_N和N个高频序列cD₁,cD₂,…,cD_N。再将分解序列进行单支 重构，得到重构子序列A_N,D₁,D₂,…,D_N。

步骤2将步骤1得到的细节子序列D₁,D₂,…,D_N按公式：

D_i⁺＝{d₁+Q,d₂+Q,d₃+Q,…,d_n+Q}Q∈N+,i＝1,2,…,N，分别平移到正数区域。再根据相空间重构理论，将其按$Y=\left(\begin{array}{c}{x}_{m+1}\\ x_{m+2}\\ ·\\ ·\\ ·\\ x_n\end{array}\right)$

的形式进行相空间重构，即转化成矩阵的形式。

步骤3将经过相空间重构的细节子序列分别输入支持向量机模型，结合遗传算法进行 参数优化，建立各细节子序列的预测模型。采用窗口滑动的形式，将尺度子序列 x_t-m,x_t-m+1,x_t-m+2,…,x_t-1输入灰色理论模型训练，预测t时刻的值，其中m为嵌入维数，建立尺 度子序列预测模型。

步骤4在进行预测分析时，根据步骤1和步骤2对数据进行处理，然后使用步骤3得 到的子序列预测模型进行预测。如输入的向量为X_t＝{x_t-m,x_t-m+1,x_t-m+2,…,x_t-1}时，

根据${\hat{x}}^{\left(0\right)}(k+1)={\hat{x}}^{\left(1\right)}(k+1)-{\hat{x}}^{\left(1\right)}\left(k\right),k=\mathrm{1,2},...,n,$此时可以得到第t时刻的尺度子序列的预测值：

${\hat{x}}^{\left(0\right)}\left(t\right)={\hat{x}}^{\left(1\right)}\left(t\right)-{\hat{x}}^{\left(1\right)}(t-1),t=\mathrm{5,6}...,n,$

根据式的最优回归函数，此时可以得到第t时刻的细节子序列 的预测值：

${\hat{Y}}_t=\underset{i=1}{\overset{l}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})K(X_t,X_i)+b,t=m+1,m+2,...,m+n$

步骤5最后，将细节子序列进行反向平移，将尺度序列模型及细节序列模型得到预测 的值根据公式进行重构，最终得到预测结果。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述 优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和 细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。