CN106371610A - 一种基于脑电信号的驾驶疲劳的检测方法 - Google Patents

Abstract

本发明涉及一种基于脑电信号的驾驶疲劳的检测方法，具体包括以下方面：(1)选取与疲劳驾驶状态相关的特征量；(2)疲劳驾驶状态的脑电信号功率谱分析；(3)疲劳驾驶状态的脑电信号样本熵分析；(4)疲劳驾驶状态的脑电信号Kc复杂度分析；(5)支持向量机(SVM)(6)最小二乘支持向量机(LS‑SVM)；(7)粒子群优化(PSO)算法设置LS‑SVM的模型训练参数C和g。本发明通过对不同驾驶状态下提取到的脑电信号分别从功率谱的角度以及使用非线性动力学中的相关方法分别进行了研究，表明在准确性和可靠性上都取得较好的效果。

Description

一种基于脑电信号的驾驶疲劳的检测方法

技术领域

本发明属于交通驾驶技术领域，具体涉及一种基于脑电信号的驾驶疲劳的检测方法

背景技术

如何从与疲劳状态紧密相关的众多特征中挑选出合适的特征量，对于进行驾驶疲劳的准确识别具有非常重要的作用，同时为了实现对驾驶疲劳进行客观而准确的识别，往往在对提取出来的脑电信号进行相关处理后，得到选取的相关的特征向量后，再使用某些分类器根据这些特征向量进行判别。支持向量机是常用的方法之一，它在实际应用中虽然能够取得不错的效果，但是它需要借助二次规划的方法进行求解，在样本数量比较大的情况下需要耗费大量的内存以及运算时间。本发明通过对提取到的不同驾驶状态的脑电信号分别从功率谱的角度以及采用非线性动力学的相关方法进行了研究，探讨了不同的驾驶状态与这些特征值之间的存在的对应关系，进行相关分析比较后选取了脑电信号的Kc复杂度、样本熵以及相关节律波的平均功率的比值等多个特征值构成了特征向量空间，作为构成驾驶疲劳识别的主要特征空间，然后采用最小二乘支持向量机对提取的样本信号的特征向量进行驾驶疲劳识别模型的训练，从而避免了解二次规划的问题，提高了模型训练的速度，同时针对LS-SVM参数通常依靠经验进行选取的问题，使用了粒子群算法(PSO)进行了相应的优化处理，最后通过实验进行了分析，表明使用PSO进行参数优化的LS-SVM进行驾驶疲劳的识别，在准确性和可靠性上都取得较好的效果，并且在实际应用中具有较好的效果。

支持向量机(SupportVectorMachine，SVM)理论作为一种新型的机器学习方法，由Vapnik等人依据统计学理论中的VC维理论以及结构风险最小化原理提出来的，它是一种监督式的学习方法，被广泛用于对统计分类以及回归分析问题中，它在解决样本数较小、非线性以及高维模式的识别问题中有着特有的优势。它的主要思想如下：

(1)它专门是针对于样本有限的情况下使用的机器学习方法，在实现结构风险最小化的原则下，即对于给定的数据通过在逼近的精度以及逼近函数的复杂性之间寻找一个比较好的折衷，以期望获得最好的推广能力。

(2)它最终是通过将要解决的问题转化为一个凸二次规划的问题进行求解，最后得到它的全局最优解。

(3)它将实际问题通过非线性变换转换到高维特征空间，通过在高维空间中构造线性决策函数以实现原空间中的非线性决策函数，这能够巧妙的解决维数问题，从而能够保证具有较好的推广能力。

发明内容

本发明提供一种基于脑电信号的驾驶疲劳检测的方法，通过对不同状态下的脑电信号进行相应分析后，选取了Kc复杂度、样本熵以及相关节律波的平均功率的比值等多个特征值构成了相应的特征向量空间，同时针对支持向量机需要利用二次规划进行计算的问题，提出了使用LS-SVM的方法进行模型的训练，提高了模型训练的速度同时针对LS-SVM参数通常依靠经验进行选取的问题，使用了粒子群算法(PSO)进行了相应的优化处理，最后通过实验进行了分析，表明使用PSO进行参数优化的LS-SVM进行驾驶疲劳的识别，在准确性和可靠性上都取得较好的效果。

为了实现上述目的，本发明技术方案如下。

一种基于脑电信号的驾驶疲劳的检测方法，具体包括以下方面：

(1)选取与疲劳驾驶状态相关的特征量：特征选择是指通过对原始测量信号进行一系列相应处理后，得到一个由一组相关特征量组成的特征集，然后从该特征集中选取某些具有代表性、分类性能最佳的物理量作为区分不同行为的一组特征量，然后使用相关的分类识别方法，利用选取的特征量将不同的行为从特征空间中分离出来，因此对于学习分类而言都需要训练集样本：

T＝{(x₁,y₁),…,(x_l,y_l)}∈(X×Y)^l (1)

其中x_i∈X＝Rⁿ,y_i∈Y＝R,i＝1,…,l；

对于符合要求的训练集就是如何选择满足要求的特征量，具体是：

(a)需要选取适当数量的能够描述不同行为模式的样本作为训练集，但是又不能选取过多的训练集，因为采用过多的训练集会导致训练时间的增长；

(b)在输入样本中，代表每一样本的行为特征量需要尽量地能描述出该行为，因此选取的特征量不能太少；

(c)为了避免复杂运算，用来描述输入样本的特征量不能什么都包括，以致于训练后的模型过于敏感；

(2)疲劳驾驶状态的脑电信号功率谱分析：对处于不同状态下脑电信号的功率谱特性进行相关分析，将采集的脑电信号的长度选择为3秒，即对采集到的数据以每3秒截取一次并进行平均功率的计算，选取了从精力充沛状态过度到疲劳状态的脑电数据8分钟，将其截取为160段，并对截取到的这些时间段的脑电信号进行傅里叶变换：

$G=\underset{T\→\mathrm{\∞}}{\lim }\frac{1}{2T}\underset{-T}{\overset{T}{\∫}}{x}_T^2\left(t\right)dt=\frac{1}{2\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\underset{T\→\mathrm{\∞}}{\lim }\frac{1}{2T}{x}_T^2\left(f\right)df---\left(2\right)$

式中，f为连续频率变量；x_T(f)为信号的频谱密度；

其中：

$p_X\left(f\right)=\underset{T\→\mathrm{\∞}}{\lim }\frac{1}{2T}{x}_T^2\left(f\right)---\left(3\right)$

在式(3)中显示了处于不同频段上的信号功率的分布情况，同时也可以表示单位时间内和单位频段上的信号的能量，因此该积分函数又被称作为功率谱密度函数；经过傅里叶变换以后，得到x_n(e^jf)，然后对频谱的模进行平方，最后得到功率谱密度的计算公式：

$p_X\left(f\right)=\frac{1}{N}|X_N\left(f\right){|}^2---\left(4\right)$

由上面可知信号在某个频段h的平均功率密度为：

$G\left(h\right)=\frac{{\∫}_{f_l}^{f_h}{p}_X\left(f\right)df}{f_h-f_l}---\left(5\right)$

其中：f_h为频带h的上限，f_l为频带h的下限，p_X(f)为信号的功率谱密度；

则不同频段的平均功率比值则定义为：

$F=\frac{G\left(h_1\right)}{G\left(h_2\right)}---\left(6\right)$

式中h₁和h₂为脑电信号的不同频带，将F_(α+θ)/β,F_α/β这两个相关节律波的平均功率的比值选作进行驾驶疲劳判别的两个特征指标；

(3)疲劳驾驶状态的脑电信号样本熵分析：利用样本熵对脑电信号进行分析的具体步骤如下：计算选取长度为1000的数据作为一个数据窗口，以500个数据为移动步长，按以上步骤计算每次移动的数据的样本熵值，然后进行下1个数据窗口的样本熵值的计算，取2个数据窗口的样本熵的平均值作为时间序列的脑电数据样本熵值，通过对处于疲劳与清醒状态的脑电数据进行分析，得到不同状态下脑电信号的样本熵值；非疲劳状态时的样本熵值大于疲劳状态下的样本熵值，表明样本熵可以用于区分不同的驾驶状态，进一步可以得出疲劳状态下脑电信号计算的到样本熵值主要集中在0.2～0.4，而清醒状态的样本熵值集中在0.5～0.8之间；由此可知随着驾驶疲劳成都的加深脑电信号的样本熵值逐渐降低的，因此可看出脑电信号的样本熵值与人体精神状态之间存在相互对应关系，因此本发明将其选作为用于驾驶疲劳识别的相关特征量；

(4)疲劳驾驶状态的脑电信号Kc复杂度分析：首先采集得到的一段待分析的脑电时间序列，设为S(S₁,S₂,...,S_n)，第一步是对该序列进行粗粒化处理，具体就是先求取该序列的平均值，然后对原序列的每个元素进行判断并重新赋值，如果序列元素S_i大于平均值时则赋值S_i'＝1，如果序列元素S_i小于平均值时则赋值S_i'＝0，然后得到重构序列S'(S₁',S'₂,...,S'_n)，重新利用S_i'重新对S中的元素进行赋值，则得到新的序列S(S₁,S₂,...,S_n)，该新序列是一串由(0，1)字符组成的新序列；下一步则对重构序列进行子串划分与界定，将按照如下的规则进行：

在一个字符串S(S₁,S₂,...,S_n)后再加一个字符串Q(q₁,q₂,...,q_n)，这样就得到组合后新的字符串序列SQ，令SQ序列删去末尾的一个字符所得字符串命名为SQv，接下来需要对Q是否属于SQv的子串进行判断，如果属于，则把下一个字符加到Q后面，增长一位再进行重复判断，直到增长后的Q不再属于SQv的一个子串，则用标记符号“.”进行序列的前后界定；下一步将“.”前的所有字符定义为新的S，重新以其后的字符开始构造Q，重复以上的判断和界定过程，直到整个字符串的序列判断结束；统计最终由“.”界定划分出的S的子串数目就是序列的Kc复杂度；

由于Kc复杂度的计算对数据的长度没有特定的要求，但是数据长度过长将影响计算的效率，而数据长度过短将不具有代表性，通过多次实验，最后本发明选择数据长度为1000，得到的不同状态下脑电信号的的Kc复杂度值结果；处于不同驾驶状态的Kc复杂度值存在着一定的差异性，其中非疲劳驾驶状态的Kc复杂度值主要集中在0.15～0.23，而处于疲劳驾驶状态得到的Kc复杂度值主要集中在0.08～0.14左右，说明受试者在非疲劳状态下进行驾驶操作的大脑复杂性高于疲劳状态下进行驾驶操作的大脑复杂性，因此表明Kc复杂度可以作为驾驶疲劳检测的另一个指标。

(5)支持向量机：考虑样本其中x_i代表第i个输入样本，y_i则是对应输入的输出；首先假设由子集y_i＝+1代表的类和y_i＝-1代表的类是能够进行线性可分的，那么用于进行分离的超平面形式的决策曲面方程可以表示为如下：

w^Tx+b＝0 (7)

x——输入变量；

w——可调的权值向量；

b——偏值；

那么式(1)则可以写成：

w^Tx+b≥0,当y_i＝+1 (8)

w^Tx+b≤0,当y_i＝-1 (9)

由式(8)和(9)则可合并为：

y_i(w^Tx+b)≥0 (10)

对于一个给定的权值向量w和偏置b，数据点到(1)定义的超平面的最近间隔称为分离边缘，用ρ表示，而支持向量机的目标则是寻找一个使得分离边缘ρ最大的超平面，而这个决策曲面则被称为最优超平面；

假设w₀和b₀分别为得到的权值向量和偏置的最优值；那么此时，输入空间里表示的多维线性决策面的最优超平面则可以表示如下：

$w_0^{T}x+b_0=0---\left(11\right)$

假如x_k对应于y_k，当满足时，那么则称(x_k,y_k)是支持向量的，那么此时支持向量到最优超平面的距离为：

$r=\left\{\begin{array}{cc}\frac{1}{\left|\right|w_0\left|\right|}& y_k=+1\\ -\frac{1}{\left|\right|w_0\left|\right|}& y_k=-1\end{array}\right.---\left(12\right)$

其中的加号表示x_k在最优超平面的正面，负号则表示最优超平面的负面，式(11)可以得到分类间隔r：

$\mathrm{\ρ}=2r=\frac{2}{\left|\right|w_0\left|\right|}---\left(13\right)$

由上式得到的最优解则对应于分类间隔最大的情况，其等价于凸二次最优化问题：

$\left\{\begin{array}{c}\min \frac{1}{2}{w}^T\·w\\ s.t.:y_i(w^T{x}_i+b)\≥1,i=1,...,m\end{array}\right.---\left(14\right)$

由上式得到的最优解则是由拉格朗日(Lagrange)函数的鞍点决定：

$L(w,b,a)=\frac{1}{2}\left|\right|w|{|}^2-\underset{i=1}{\overset{N}{\Σ}}{a}_i\[y_i(w^T{x}_i+b)-1\]---\left(15\right)$

式中，

a_i,a_j——拉格朗日乘子；

由最优解条件(KKT条件)可得：

$\left.\begin{array}{c}\frac{\∂L(w,b,a)}{\∂b}=0\⇒\underset{i=1}{\overset{m}{\Σ}}{a}_i{y}_i=0\\ \frac{\∂L(w,b,a)}{\∂w}=0\⇒w=\underset{i=1}{\overset{m}{\Σ}}{a}_i{x}_i{y}_i\end{array}\right\}---\left(16\right)$

根据KKT(Karush-Kuhn-Tucker)条件，相应的拉格朗日乘子a_i应满足如下条件：

a_i[y_i(w^Tx_i+b)-1]＝0 (17)

从上式可以看出只有支持向量的a_i不为零；

根据核函数变换，则式(14)和(16)可转化为对偶问题求解，即：

$\left\{\begin{array}{c}W\left(a\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{a}_i{a}_j{y}_i{y}_j{x}_i^T{x}_j\\ st.\underset{i=1}{\overset{N}{\Σ}}{y}_i{a}_i=0,a_i\≥0,i=1,\mathrm{...},N\end{array}\right.---\left(18\right)$

然后由(16)可以求出最优权值向量w₀，以及利用式子(11)可以得到b₀，于是关于待测样本x，通过只算sgn[(w^Tx)+b]就可以判断其所属的类别；

但是实际应用当中的大部分问题都是属于线性不可分的情况，而造成这样的原因也是各种各样的，可能在样本的产生时造成一定误差，也可能在进行统计记录，标记的错误导致样本数据的误差等，这样也可能使的样本出现线性不可分的情况；为了处理当存在少量样本无法用线性方法分开的情况，因此在分类超平面的定义中引入松弛变量ξ_i≥0，即允许少量误差的存在，那么(14)的约束条件变为如下形式：

y_i(w^Tx+b)≥1-ξ_i,i＝1,2,...N (19)

此时的目标函数则变为：

$\left\{\begin{array}{c}\min \frac{1}{2}{w}^T\·w+C\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i\\ y_i(w^{T}x+b)\≥1-{\mathrm{\ξ}}_i,i=1,2,...N\end{array}\right.---\left(20\right)$

其中C为惩罚因子且C≥0，当选取的参数C比较大时，表示支持向量机模型对训练样本的质量具有较高的信心，相反，当选取的参数C比较小时，则认为训练样本中存在噪声；由(20)的到相应的式子如下所示：

$L(w,b,a)=\frac{1}{2}\left|\right|w|{|}^2+C\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i-\underset{i=1}{\overset{N}{\Σ}}{a}_i\[y_i(w^T{x}_i+b)+{\mathrm{\ξ}}_i-1\]-\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\β}}_i{\mathrm{\ξ}}_i---\left(21\right)$

那么根据KKT条件转化为其对偶问题如下:

$W\left(a\right)=\underset{i=1}{\overset{m}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{m}{\Σ}}{a}_i{a}_j{y}_i{y}_j{x}_i^T{x}_j---\left(22\right)$

$st.\underset{i=1}{\overset{m}{\Σ}}{y}_i{a}_i=0,0\≤a_i\≤C,i=1,...,m---\left(23\right)$

将输入空间Rⁿ中的点利用非线性映射Φ:Rⁿ→H映射到高维的Hilbert空间当中去，然后再使用线性可分的方法进行最优分类面的求解，而在这个求解的过程中，仅仅涉及到内积的运算，并且计算量也只依赖于训练样本数的个数，因此这个过程的就是如何找到一个核函数K(x_i,x_j),x_i,x_j∈Rⁿ，使得该核函数满足K(x_i,x_j)＝Φ^T(x_i)Φ(x_j)，这样就可以保证即使在不知道变换函数Φ形式的情况下也能够求得最终的结果，同时能够很好地避免维数灾难的发生，因此在最大间隔情况下，非线性支持向量机的目标函数为：

$W\left(a\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{a}_i{a}_j{y}_i{y}_{j}K(x_i,x_j)---\left(24\right)$

同样可得知，非线性支持向量机的最优化问题如下所示：

$\left\{\begin{array}{c}\min \frac{1}{2}{w}^T\·w+C\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i\\ st:\underset{i=1}{\overset{m}{\Σ}}{y}_i{a}_i=0,0\≤a_i\≤C,i=1,...,m\end{array}\right.---\left(25\right)$

它的对偶问题则为：

$\left\{\begin{array}{c}\max W\left(a\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{a}_i{a}_j{y}_i{y}_{j}K(x_i,x_j)\\ st:\underset{i=1}{\overset{N}{\Σ}}{y}_i{a}_i=0,0\≤a_i\≤C,i=1,\mathrm{...},N\end{array}\right.---\left(26\right)$

相应得到的分类函数为：

$f\left(x\right)=\mathrm{sgn}\[\left(w^T\mathrm{\Φ}\right(x\left)\right)+b\]=\mathrm{sgn}\[\underset{i=1}{\overset{N}{\Σ}}{y}_i{a}_{i}K(x_i,x_j)+b\]---\left(27\right)$

(6)最小二乘支持向量机：最小二乘支持向量机算法的描述如下：首先选择一个非线性的变换把n维输入和1维输出的样本向量(x₁,y₁)，(x₂,y₂)，…，(x_N,y_N)，x_i∈Rⁿ,y_i∈R，i＝1,2,...,N，l从原空间映射到高维特征空间，则分类LSSVM的优化问题则变为：

即将支持向量机的不等式约束改为等式约束，构造拉格朗日函数进行求解：

$L(w,b,a)=\frac{1}{2}\left|\right|w|{|}^2+\frac{C}{2}\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i^2-\underset{i=1}{\overset{N}{\Σ}}{a}_i\[y_i(w^T{x}_i+b)-1+{\mathrm{\ξ}}_i\]---\left(29\right)$

其中a_i——拉格朗日乘子；

最优解满足KTT条件则：进而可以得到下面的方程组：

对于i＝1,2,...,N，通过式子(30)消去w和ξ得到如下线性系统：

$\left[\begin{array}{cc}0& Y^T\\ Y& {\mathrm{ZZ}}^T+C^{-1}I\end{array}\right]\left[\begin{array}{c}b\\ a\end{array}\right]=\left[\begin{array}{c}0\\ 1\end{array}\right]---\left(31\right)$

式中：

令Ω＝ZZ^T+C^-1I，求解5-13可得到：

$\left\{\begin{array}{c}a={\mathrm{\Ω}}^{-1}(Y-b1)\\ b=\frac{1^T{\mathrm{\Ω}}^{-1}Y}{1^T{\mathrm{\Ω}}^{-1}1}\end{array}\right.---\left(33\right)$

最终得到的决策函数如下：

$f\left(x\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_{i}K(x,x_i)+b;---\left(34\right)$

(7)粒子群优化(PSO)算法设置LS-SVM的模型训练参数C和g：采用最小二乘法支持向量机进行模型建立的过程来需要确定两个参数，即惩罚参数C和核函数参数g，不同参数对模型性能以及泛化能力有不同的影响；因此，在建模的时候，选取适当的参数能够极大地改善识别的准确度，接下来具体研究LS-SVM的参数影响：惩罚参数C可以理解为对超出样本数据的一个管制制度，他直接影响模型的复杂程度和稳定性，惩罚参数C能够使模型在复杂度和训练误差之间做出一个平衡，使模型的推广能力得到保证，若惩罚参数C过小，则对超出的样本数据惩罚就小，训练的误差就会变大，模型的推广能力就会变差，当使用模型对新的数据进行学习时，错分率就会很高，而当惩罚参数C过大时，训练学习精度相应提高，但会对训练样本产生过度的依赖，对新数据的认识度不够，导致错分率增加，也会影响模型的推广能力，此外，惩罚参数C值也会对算法处理样本中的离群点(数据集合中一些元素和其他元素很不相似或不一致)产生影响，选取合适的值就能在一定程度上抗干扰，使模型稳定；

核函数参数g对支持向量的性能有很大的影响，对于高斯核中的核宽度δ而言，它反映了支持向量之间的相关程度，δ过小，支持向量之间的联系比较松弛，机器学习相对复杂，推广能力得不到保证，δ过大，支持向量之间的影响过强，预测模型很难达到足够的精度，所以在选取参数时，应综合考虑他们之间的相互影响；

因此，只有恰当的选择LS-SVM的参数才能得到最小的经验误差，同时提高模型的训练速度以及模型的泛化能力，使得预测结果更为准确；其中传统的模型参数选取主要通过交叉验证法以及经验选择方法等，但是使用经验选择发通常需要获知该重逢的经验知识，在通过判断进行参数的设置，要求比较高，因此一般不采用，而交叉验证法较为普遍，是大多数所采用的方法，便于实现，步骤简单，但是通过这种方法获得的参数精确度一般，并且需要花费较长的时间，本发明采用智能优化算法中的粒子群算法进行对模型参数的寻优，使其获得较好的性能，粒子群算法具体如下：首先在可解空间内初始化一群粒子，而每一个粒子都同时拥有位置、速度和适应度值三个特征，而其中的适应度值的大小则决定该粒子的优劣，同时粒子的移动主要是通过跟踪个体极值Pbest、群体极值Gbest以及过去速度的惯性来更新个体位置，其中个体极值Pbest是指单个粒子所经历的适应度值最优的位置，群体极值Gbest是指种群中所有粒子所经历适应度值最优的位置，粒子每迭代一次就移动一次位置，其适应度值就相应的需要重新计算一次，然后通过对新粒子个体极值的适应度值、群体极值的适应度值与迭代之前的极值比较来确定是否更新粒子的个体极值和群体极值的位置；

假设在一个D维的空间中，其中X＝(X₁,X₂,...,X_n)为由n个粒子组成的种群，其中向量X_i＝(x_i1,x_i2,...,x_iD)^T表示第i个粒子，对于第i个粒子，它的适应度值一般就是根据粒子X_i在目标函数中的值，速度为V_i＝(v_i1,v_i2,...,v_iD)^T，其个体极值为P_i＝(p_i1,p_i2,...,p_iD)^T，种群的全局极值为P_g＝(p_g1,p_g2,...,p_gD)^T

在每一次的迭代过程中，每个粒子都是通过个体极值、全局极值以及惯性速度来更新它的速度和位置信息，其更新公式(35)和(36)如下所示：

$v_{id}^{k+1}={\mathrm{wv}}_{id}^k+c_1{r}_1(p_{id}^k-x_{id}^k)+c_2{r}_2(p_{gd}^k-x_{id}^k)---\left(35\right)$

$x_{id}^{k+1}=x_{id}^k+v_{id}^{k+1}---\left(36\right)$

式(24)中，w为惯性权重，当w较大时有利于群体在较大的范围内进行搜索，从而避免陷入寻找局部最优，而较小的w能够保证算法能够快速收敛到最优值；d＝1,2,...,D，i＝1,2,...,n，k为当前迭代次数；v_id为粒子的速度；c₁和c₂为非负常数，称为加速度因子，是每一个粒子在P_i和P_g位置的随机加速权重，r₁和r₂是[0，1]之间的随机数；为了防止粒子群的盲目搜索，一般把粒子群的位置限定在[-X_max,X_max]和速度限定在[-V_max,V_max]的范围内；

根据上面的分析，本发明使用粒子群算法优化参数C和g，具体步骤如下：

(a)初始化粒子群并随机产生所有粒子的位置和速度并确定粒子的P_i和P_g；

(b)将每个粒子当前最优位置P_i及种群中所有粒子所经历的最好位置进P_g行比较，如果这个粒子的位置最优，则将其设置为当前的最好位置P_g，否则，P_g保持不变；

(c)更新粒子的速度、位置；

(d)若没有达到结束条件，如：预设的运算精度、迭代次数，便返回步骤(b)；

(e)开始下一次的迭代计算，否则取当前的P_g为最优解。

该发明的有益效果在于：本发明通过对脑电信号从功率谱角度，非线性动力学的相关方法等多个角度进行了其与疲劳状态的相关性，最后从中选取了多个特征值构成了采用LS-SVM进行驾驶疲劳识别的向量空间，解决了仅仅从单一角度进行分析带来的识别准确度的不高的问题，同时使用粒子群算法对最小二乘向量机参数进行了优化，解决了依靠经验选取得到参数导致的准确度精度偏低的情况，最后通过实验进行了分析，表明使用PSO进行参数优化的LS-SVM进行驾驶疲劳的识别，在准确性和可靠性上都取得较好的效果。

附图说明

图1为本发明实施例中所采用的PSO优化算法流程图。

具体实施例

上面结合附图和具体实施例进一步阐述本发明。

(1)选取与疲劳驾驶状态相关的特征量：特征选择是指通过对原始测量信号进行一系列相应处理后，得到一个由一组相关特征量组成的特征集，然后从该特征集中选取某些具有代表性、分类性能最佳的物理量作为区分不同行为的一组特征量，然后使用相关的分类识别方法，利用选取的特征量将不同的行为从特征空间中分离出来，因此对于学习分类而言都需要训练集样本：

T＝{(x₁,y₁),…,(x_l,y_l)}∈(X×Y)^l (1)

其中x_i∈X＝Rⁿ,y_i∈Y＝R,i＝1,…,l；

对于符合要求的训练集就是如何选择满足要求的特征量，具体是：

(a)需要选取适当数量的能够描述不同行为模式的样本作为训练集，但是又不能选取过多的训练集，因为采用过多的训练集会导致训练时间的增长；

(b)在输入样本中，代表每一样本的行为特征量需要尽量地能描述出该行为，因此选取的特征量不能太少；

(c)为了避免复杂运算，用来描述输入样本的特征量不能什么都包括，以致于训练后的模型过于敏感；

(2)疲劳驾驶状态的脑电信号功率谱分析：对处于不同状态下脑电信号的功率谱特性进行相关分析，将采集的脑电信号的长度选择为3秒，即对采集到的数据以每3秒截取一次并进行平均功率的计算，选取了从精力充沛状态过度到疲劳状态的脑电数据8分钟，将其截取为160段，并对截取到的这些时间段的脑电信号进行傅里叶变换：

$G=\underset{T\→\mathrm{\∞}}{\lim }\frac{1}{2T}\underset{-T}{\overset{T}{\∫}}{x}_T^2\left(t\right)dt=\frac{1}{2\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\underset{T\→\mathrm{\∞}}{\lim }\frac{1}{2T}{x}_T^2\left(f\right)df---\left(2\right)$

式中，f为连续频率变量；x_T(f)为信号的频谱密度；

其中：

$p_X\left(f\right)=\underset{T\→\mathrm{\∞}}{\lim }\frac{1}{2T}{x}_T^2\left(f\right)---\left(3\right)$

在式(3)中显示了处于不同频段上的信号功率的分布情况，同时也可以表示单位时间内和单位频段上的信号的能量，因此该积分函数又被称作为功率谱密度函数；经过傅里叶变换以后，得到x_n(e^jf)，然后对频谱的模进行平方，最后得到功率谱密度的计算公式：

$p_X\left(f\right)=\frac{1}{N}|X_N\left(f\right){|}^2---\left(4\right)$

由上面可知信号在某个频段h的平均功率密度为：

$G\left(h\right)=\frac{{\∫}_{f_l}^{f_h}{p}_X\left(f\right)df}{f_h-f_l}---\left(5\right)$

其中：f_h为频带h的上限，f_l为频带h的下限，p_X(f)为信号的功率谱密度；

则不同频段的平均功率比值则定义为：

$F=\frac{G\left(h_1\right)}{G\left(h_2\right)}---\left(6\right)$

式中h₁和h₂为脑电信号的不同频带，将F_(α+θ)/β,F_α/β这两个相关节律波的平均功率的比值选作进行驾驶疲劳判别的两个特征指标；

(3)疲劳驾驶状态的脑电信号样本熵分析：利用样本熵对脑电信号进行分析的具体步骤如下：计算选取长度为1000的数据作为一个数据窗口，以500个数据为移动步长，按以上步骤计算每次移动的数据的样本熵值，然后进行下1个数据窗口的样本熵值的计算，取2个数据窗口的样本熵的平均值作为时间序列的脑电数据样本熵值，通过对处于疲劳与清醒状态的脑电数据进行分析，得到不同状态下脑电信号的样本熵值；非疲劳状态时的样本熵值大于疲劳状态下的样本熵值，表明样本熵可以用于区分不同的驾驶状态，进一步可以得出疲劳状态下脑电信号计算的到样本熵值主要集中在0.2～0.4，而清醒状态的样本熵值集中在0.5～0.8之间；由此可知随着驾驶疲劳成都的加深脑电信号的样本熵值逐渐降低的，因此可看出脑电信号的样本熵值与人体精神状态之间存在相互对应关系，因此本发明将其选作为用于驾驶疲劳识别的相关特征量；

(4)疲劳驾驶状态的脑电信号Kc复杂度分析；：首先采集得到的一段待分析的脑电时间序列，设为S(S₁,S₂,...,S_n)，第一步是对该序列进行粗粒化处理，具体就是先求取该序列的平均值，然后对原序列的每个元素进行判断并重新赋值，如果序列元素S_i大于平均值时则赋值S_i'＝1，如果序列元素S_i小于平均值时则赋值S_i'＝0，然后得到重构序列S'(S₁',S'₂,...,S'_n)，重新利用S_i'重新对S中的元素进行赋值，则得到新的序列S(S₁,S₂,...,S_n)，该新序列是一串由(0，1)字符组成的新序列；下一步则对重构序列进行子串划分与界定，将按照如下的规则进行：

在一个字符串S(S₁,S₂,...,S_n)后再加一个字符串Q(q₁,q₂,...,q_n)，这样就得到组合后新的字符串序列SQ，令SQ序列删去末尾的一个字符所得字符串命名为SQv，接下来需要对Q是否属于SQv的子串进行判断，如果属于，则把下一个字符加到Q后面，增长一位再进行重复判断，直到增长后的Q不再属于SQv的一个子串，则用标记符号“.”进行序列的前后界定；下一步将“.”前的所有字符定义为新的S，重新以其后的字符开始构造Q，重复以上的判断和界定过程，直到整个字符串的序列判断结束；统计最终由“.”界定划分出的S的子串数目就是序列的Kc复杂度。

为了更清晰说明计算过程，举例进行说明，比如100101，则其复杂度可由下列步骤计算：

第一个字符→1.

2、S＝1，Q＝0，SQ＝10，SQv＝1，Q不属于SQv→1.0.

3、S＝10，Q＝0，SQ＝100，SQv＝10，Q属于SQv→1.0.0

4·S＝10，Q＝01，SQ＝1001，SQv＝100，Q不属于QSv→1.0.01.

5.S＝1001，Q＝0，SQ＝10010，SQv＝1001，Q属于SQv→1.0.01.0

6.S＝1001，Q＝01，SQ＝100101，SQv＝10010，Q属于SQv→1.0.01.01

在第六步完成了整个字符串判断，结果中“.”界定划分出S的子串数目是4个，所以这个序列的复杂度c(n)＝4，即c(100101)＝4。

由于结果并无明显的上界，不易进行比较和分析，所以Lempel和Ziv又对S_i的(0，1)字符串序列进行了进一步的研究，发现几乎所有的序列c(n)都趋向于一定值：

$\underset{n\→\mathrm{\∞}}{\lim }c\left(n\right)=b\left(n\right)=n/\log n---\left(7\right)$

所以c(n)体现的是随机序列的渐近行为，为了能够更加直观得到利用复杂度来对序列的结果进行分析，则利用趋向的定值b(n)对c(n)进行归一化处理，得到新的复杂度结果即如下所示：

C(n)＝c(n)/b(n) (8)

由于Kc复杂度的计算对数据的长度没有特定的要求，但是数据长度过长将影响计算的效率，而数据长度过短将不具有代表性，通过多次实验，最后本发明选择数据长度为1000，得到的不同状态下脑电信号的的Kc复杂度值结果；处于不同驾驶状态的Kc复杂度值存在着一定的差异性，其中非疲劳驾驶状态的Kc复杂度值主要集中在0.15～0.23，而处于疲劳驾驶状态得到的Kc复杂度值主要集中在0.08～0.14左右，说明受试者在非疲劳状态下进行驾驶操作的大脑复杂性高于疲劳状态下进行驾驶操作的大脑复杂性，因此表明Kc复杂度可以作为驾驶疲劳检测的另一个指标。

(5)支持向量机：考虑样本其中x_i代表第i个输入样本，y_i则是对应输入的输出；首先假设由子集y_i＝+1代表的类和y_i＝-1代表的类是能够进行线性可分的，那么用于进行分离的超平面形式的决策曲面方程可以表示为如下：

w^Tx+b＝0 (9)

x——输入变量；

w——可调的权值向量；

b——偏值；

那么式(1)则可以写成：

w^Tx+b≥0,当y_i＝+1 (10)

w^Tx+b≤0,当y_i＝-1 (11)

由式(10)和(11)则可合并为：

y_i(w^Tx+b)≥0 (12)

对于一个给定的权值向量w和偏置b，数据点到(1)定义的超平面的最近间隔称为分离边缘，用ρ表示，而支持向量机的目标则是寻找一个使得分离边缘ρ最大的超平面，而这个决策曲面则被称为最优超平面；

假设w₀和b₀分别为得到的权值向量和偏置的最优值；那么此时，输入空间里表示的多维线性决策面的最优超平面则可以表示如下：

$w_0^{T}x+b_0=0---\left(13\right)$

假如x_k对应于y_k，当满足时，那么则称(x_k,y_k)是支持向量的，那么此时支持向量到最优超平面的距离为：

$r=\left\{\begin{array}{cc}\frac{1}{\left|\right|w_0\left|\right|}& y_k=+1\\ -\frac{1}{\left|\right|w_0\left|\right|}& y_k=-1\end{array}\right.---\left(14\right)$

其中的加号表示x_k在最优超平面的正面，负号则表示最优超平面的负面，式(13)可以得到分类间隔r：

$\mathrm{\ρ}=2r=\frac{2}{\left|\right|w_0\left|\right|}---\left(15\right)$

由上式得到的最优解则对应于分类间隔最大的情况，其等价于凸二次最优化问题：

$\left\{\begin{array}{c}\min \frac{1}{2}{w}^T\·w\\ s.t.:y_i(w^T{x}_i+b)\≥1,i=1,...,m\end{array}\right.---\left(16\right)$

由上式得到的最优解则是由拉格朗日(Lagrange)函数的鞍点决定：

$L(w,b,a)=\frac{1}{2}\left|\right|w|{|}^2-\underset{i=1}{\overset{N}{\Σ}}{a}_i\[y_i(w^T{x}_i+b)-1\]---\left(17\right)$

式中，

a_i,a_j——拉格朗日乘子；

由最优解条件(KKT条件)可得：

$\left.\begin{array}{c}\frac{\∂L(w,b,a)}{\∂b}=0\⇒\underset{i=1}{\overset{m}{\Σ}}{a}_i{y}_i=0\\ \frac{\∂L(w,b,a)}{\∂w}=0\⇒w=\underset{i=1}{\overset{m}{\Σ}}{a}_i{x}_i{y}_i\end{array}\right\}---\left(18\right)$

根据KKT(Karush-Kuhn-Tucker)条件，相应的拉格朗日乘子a_i应满足如下条件：

a_i[y_i(w^Tx_i+b)-1]＝0 (19)

从上式可以看出只有支持向量的a_i不为零；

根据核函数变换，则式(16)和(18)可转化为对偶问题求解，即：

$\left\{\begin{array}{c}W\left(a\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{a}_i{a}_j{y}_i{y}_j{x}_i^T,x_j\\ st.\underset{i=1}{\overset{N}{\Σ}}{y}_i{a}_i=0,a_i\≥0,i=1,\mathrm{...},N\end{array}\right.---\left(20\right)$

然后由(18)可以求出最优权值向量w₀，以及利用式子(13)可以得到b₀，于是关于待测样本x，通过只算sgn[(w^Tx)+b]就可以判断其所属的类别；

但是实际应用当中的大部分问题都是属于线性不可分的情况，而造成这样的原因也是各种各样的，可能在样本的产生时造成一定误差，也可能在进行统计记录，标记的错误导致样本数据的误差等，这样也可能使的样本出现线性不可分的情况；为了处理当存在少量样本无法用线性方法分开的情况，因此在分类超平面的定义中引入松弛变量ξ_i≥0，即允许少量误差的存在，那么(16)的约束条件变为如下形式：

y_i(w^Tx+b)≥1-ξ_i,i＝1,2,...N (21)

此时的目标函数则变为：

$\left\{\begin{array}{c}\min \frac{1}{2}{w}^T\·w+C\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i\\ y_i(w^{T}x+b)\≥1-{\mathrm{\ξ}}_i,i=1,2,...N\end{array}\right.---\left(22\right)$

其中C为惩罚因子且C≥0，当选取的参数C比较大时，表示支持向量机模型对训练样本的质量具有较高的信心，相反，当选取的参数C比较小时，则认为训练样本中存在噪声；由(22)的到相应的式子如下所示：

$L(w,b,a)=\frac{1}{2}\left|\right|w|{|}^2+C\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i-\underset{i=1}{\overset{N}{\Σ}}{a}_i\[y_i(w^T{x}_i+b)+{\mathrm{\ξ}}_i-1\]-\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\β}}_i{\mathrm{\ξ}}_i---\left(23\right)$

那么根据KKT条件转化为其对偶问题如下:

$W\left(a\right)=\underset{i=1}{\overset{m}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{m}{\Σ}}{a}_i{a}_j{y}_i{y}_j{x}_i^T{x}_j---\left(24\right)$

$st.\underset{i=1}{\overset{m}{\Σ}}{y}_i{a}_i=0,0\≤a_i\≤C,i=1,...,m---\left(25\right)$

将输入空间Rⁿ中的点利用非线性映射Φ:Rⁿ→H映射到高维的Hilbert空间当中去，然后再使用线性可分的方法进行最优分类面的求解，而在这个求解的过程中，仅仅涉及到内积的运算，并且计算量也只依赖于训练样本数的个数，因此这个过程的就是如何找到一个核函数K(x_i,x_j),x_i,x_j∈Rⁿ，使得该核函数满足K(x_i,x_j)＝Φ^T(x_i)Φ(x_j)，这样就可以保证即使在不知道变换函数Φ形式的情况下也能够求得最终的结果，同时能够很好地避免维数灾难的发生，因此在最大间隔情况下，非线性支持向量机的目标函数为：

$W\left(a\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{a}_i{a}_j{y}_i{y}_{j}K(x_i,x_j)---\left(26\right)$

同样可得知，非线性支持向量机的最优化问题如下所示：

$\left\{\begin{array}{c}\min \frac{1}{2}{w}^T\·w+C\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i\\ st:\underset{i=1}{\overset{m}{\Σ}}{y}_i{a}_i=0,0\≤a_i\≤C,i=1,...,m\end{array}\right.---\left(27\right)$

它的对偶问题则为：

$\left\{\begin{array}{c}\max W\left(a\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{a}_i{a}_j{y}_i{y}_{j}K(x_i,x_j)\\ st:\underset{i=1}{\overset{N}{\Σ}}{y}_i{a}_i=0,0\≤a_i\≤C,i=1,\mathrm{...},N\end{array}\right.---\left(28\right)$

相应得到的分类函数为：

$f\left(x\right)=\mathrm{sgn}\[\left(w^T\mathrm{\Φ}\right(x\left)\right)+b\]=\mathrm{sgn}\[\underset{i=1}{\overset{N}{\Σ}}{y}_i{a}_{i}K(x_i,x_j)+b\]---\left(29\right)$

虽然支持向量机的方法在很多应用中都能取得不错的效果，但是在SVM的求解过程中实际上是通过使用二次规划的方法来进行求解的，而在二次规划的求解过程中，当样本数目很大时，计算将耗费大量的机器内存以及运算时间。因此，在SVM的研究的基础上，提出了最小二乘支持向量机，它采用训练误差的平方代替松弛变量，同时将不等式约束变成等式约束，从而问题的求解则变成了对线性方程组的求解，从而避免了解二次规划的问题，并且能够提高求解的速度

(6)最小二乘支持向量机：最小二乘支持向量机算法的描述如下：首先选择一个非线性的变换把n维输入和1维输出的样本向量(x₁,y₁)，(x₂,y₂)，…，(x_N,y_N)，x_i∈Rⁿ,y_i∈R，i＝1,2,...,N，l从原空间映射到高维特征空间，则分类LSSVM的优化问题则变为：

即将支持向量机的不等式约束改为等式约束，构造拉格朗日函数进行求解：

$L(w,b,a)=\frac{1}{2}\left|\right|w|{|}^2+\frac{C}{2}\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i^2-\underset{i=1}{\overset{N}{\Σ}}{a}_i\[y_i(w^T{x}_i+b)-1+{\mathrm{\ξ}}_i\]---\left(31\right)$

其中a_i——拉格朗日乘子；

最优解满足KTT条件则：进而可以得到下面的方程组：

对于i＝1,2,...,N，通过式子(32)消去w和ξ得到如下线性系统：

$\left[\begin{array}{cc}0& Y^T\\ Y& {\mathrm{ZZ}}^T+C^{-1}I\end{array}\right]\left[\begin{array}{c}b\\ a\end{array}\right]=\left[\begin{array}{c}0\\ 1\end{array}\right]---\left(33\right)$

式中：

令Ω＝ZZ^T+C^-1I，求解可得到：

$\left\{\begin{array}{c}a={\mathrm{\Ω}}^{-1}(Y-b1)\\ b=\frac{1^T{\mathrm{\Ω}}^{-1}Y}{1^T{\mathrm{\Ω}}^{-1}1}\end{array}\right.---\left(35\right)$

最终得到的决策函数如下：

$f\left(x\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_{i}K(x,x_i)+b;---\left(36\right)$

由上面可以知道，最小二乘支持向量机是支持向量机的一种改进，其中最小二乘支持向量机把复杂的非线性二次规划问题化解成线性问题来进行求解，是支持向量机在二次损失函数下的一种形式，它既发扬了支持向量机本身的优点，也发挥出自己一些独特的优点，例如计算相对简单，计算速度提高，减少计算机内存等优点，虽然LSSVM算法相比其他方式有着更强的泛化能力，但它对模型本身的内部参数有着很强的依赖性，其中参数值的选取很大程度的关系着模型的性能以及它的识别精度。

(7)粒子群优化(PSO)算法设置LS-SVM的模型训练参数C和g：采用最小二乘法支持向量机进行模型建立的过程来需要确定两个参数，即惩罚参数C和核函数参数g，不同参数对模型性能以及泛化能力有不同的影响；因此，在建模的时候，选取适当的参数能够极大地改善识别的准确度，接下来具体研究LS-SVM的参数影响：惩罚参数C可以理解为对超出样本数据的一个管制制度，他直接影响模型的复杂程度和稳定性，惩罚参数C能够使模型在复杂度和训练误差之间做出一个平衡，使模型的推广能力得到保证，若惩罚参数C过小，则对超出的样本数据惩罚就小，训练的误差就会变大，模型的推广能力就会变差，当使用模型对新的数据进行学习时，错分率就会很高，而当惩罚参数C过大时，训练学习精度相应提高，但会对训练样本产生过度的依赖，对新数据的认识度不够，导致错分率增加，也会影响模型的推广能力，此外，惩罚参数C值也会对算法处理样本中的离群点(数据集合中一些元素和其他元素很不相似或不一致)产生影响，选取合适的值就能在一定程度上抗干扰，使模型稳定；

核函数参数g对支持向量的性能有很大的影响，对于高斯核中的核宽度δ而言，它反映了支持向量之间的相关程度，δ过小，支持向量之间的联系比较松弛，机器学习相对复杂，推广能力得不到保证，δ过大，支持向量之间的影响过强，预测模型很难达到足够的精度，所以在选取参数时，应综合考虑他们之间的相互影响；

因此，只有恰当的选择LS-SVM的参数才能得到最小的经验误差，同时提高模型的训练速度以及模型的泛化能力，使得预测结果更为准确；其中传统的模型参数选取主要通过交叉验证法以及经验选择方法等，但是使用经验选择发通常需要获知该重逢的经验知识，在通过判断进行参数的设置，要求比较高，因此一般不采用，而交叉验证法较为普遍，是大多数所采用的方法，便于实现，步骤简单，但是通过这种方法获得的参数精确度一般，并且需要花费较长的时间，本发明采用智能优化算法中的粒子群算法进行对模型参数的寻优，使其获得较好的性能，粒子群算法具体如下：首先在可解空间内初始化一群粒子，而每一个粒子都同时拥有位置、速度和适应度值三个特征，而其中的适应度值的大小则决定该粒子的优劣，同时粒子的移动主要是通过跟踪个体极值Pbest、群体极值Gbest以及过去速度的惯性来更新个体位置，其中个体极值Pbest是指单个粒子所经历的适应度值最优的位置，群体极值Gbest是指种群中所有粒子所经历适应度值最优的位置，粒子每迭代一次就移动一次位置，其适应度值就相应的需要重新计算一次，然后通过对新粒子个体极值的适应度值、群体极值的适应度值与迭代之前的极值比较来确定是否更新粒子的个体极值和群体极值的位置；

假设在一个D维的空间中，其中X＝(X₁,X₂,...,X_n)为由n个粒子组成的种群，其中向量X_i＝(x_i1,x_i2,...,x_iD)^T表示第i个粒子，对于第i个粒子，它的适应度值一般就是根据粒子X_i在目标函数中的值，速度为V_i＝(v_i1,v_i2,...,v_iD)^T，其个体极值为P_i＝(p_i1,p_i2,...,p_iD)^T，种群的全局极值为P_g＝(p_g1,p_g2,...,p_gD)^T

在每一次的迭代过程中，每个粒子都是通过个体极值、全局极值以及惯性速度来更新它的速度和位置信息，其更新公式(37)和(38)如下所示：

$v_{id}^{k+1}={\mathrm{wv}}_{id}^k+c_1{r}_1(p_{id}^k-x_{id}^k)+c_2{r}_2(p_{gd}^k-x_{id}^k)---\left(37\right)$

$x_{id}^{k+1}=x_{id}^k+v_{id}^{k+1}---\left(38\right)$

式(24)中，w为惯性权重，当w较大时有利于群体在较大的范围内进行搜索，从而避免陷入寻找局部最优，而较小的w能够保证算法能够快速收敛到最优值；d＝1,2,...,D，i＝1,2,...,n，k为当前迭代次数；v_id为粒子的速度；c₁和c₂为非负常数，称为加速度因子，是每一个粒子在P_i和P_g位置的随机加速权重，r₁和r₂是[0，1]之间的随机数；为了防止粒子群的盲目搜索，一般把粒子群的位置限定在[-X_max,X_max]和速度限定在[-V_max,V_max]的范围内；

根据上面的分析，本发明使用粒子群算法优化参数C和g，具体步骤如下：

(a)初始化粒子群并随机产生所有粒子的位置和速度并确定粒子的P_i和P_g；

(b)将每个粒子当前最优位置P_i及种群中所有粒子所经历的最好位置进P_g行比较，如果这个粒子的位置最优，则将其设置为当前的最好位置P_g，否则，P_g保持不变；

(c)更新粒子的速度、位置；

(d)若没有达到结束条件，如：预设的运算精度、迭代次数，便返回步骤(b)；

(e)开始下一次的迭代计算，否则取当前的P_g为最优解。

上面提出的最小二乘支持向量的模型训练的过程中，惩罚系数C和核函数参数g都需要由实验者依据其经验进行事先确定，然而这两个参数选择的合适与否与建模后的的识别准确度息息相关。本发明分别采用交叉验证法以及采用粒子群优化的方法分别进行了相关实验，寻找出相应的参数(C，g)，其中在使用PSO-LSSVM的算法中参数设置，初始化种群规模为30，最大迭代次数为200，c1＝c2＝2，惯性权重w0.9递减至0.4，初始化惩罚因子C的取值范围是[0.001，500]，核函数参数g的取值范围是[0.1，1000]，最大速度限制为10，由此得到的最优参数对比表如下所示：

选取由5名体型、特征各不相同，身体健康的志愿者进行模拟驾驶状态，分别采集处于不同状态下的脑电数据，并对采集到的脑电数据进行相应的后，采用第四章提供的特征提取方法提取不同状态下特征值，由此组成相关的特征向量X_i＝(F_α/β,F_(α+θ)/β,SampEn,Kc复杂度)，然后将实验中所取得的总体样本数据的75％作为训练样本，并对其进行样本分类，然后将剩余的25％作为测试样本。分别使用SVM、LSSVM以及PSO-LSSVM算法进行实验，得到对比分析结果。

在选取相同测试集下，使用SVM进行分类的识别准确率达到89.34％，而使用LSSVM进行分类的准确率为85.47％，使用P粒子优化算法对惩罚系数C和核函数参数g进行优化后，使用PSO-LSSVM的对驾驶疲劳进行识别的准确率有了进一步的提高达到92.67％，表明，使用粒子算法对参数进行优化能够提高识别的准确率。

Claims (1)

1.一种基于脑电信号的驾驶疲劳的检测方法，其特征在于：具体包括以下方面：

(1)选取与疲劳驾驶状态相关的特征量：特征选择是指通过对原始测量信号进行一系列相应处理后，得到一组由相关特征值组成的特征集，然后从该特征集中选取某些具有代表性、分类性能最佳的物理量作为区分不同行为的一组特征量，然后使用相关的分类识别方法，利用选取的特征量将不同的行为从特征空间中分离出来，因此对于学习分类而言都需要训练集样本：

T＝{(x₁,y₁),…,(x_l,y_l)}∈(X×Y)^l (1)

其中x_i∈X＝Rⁿ,y_i∈Y＝R,i＝1,…,l；

对于符合要求的训练集就是如何选择满足要求的特征量，具体是：

(a)需要选取适当数量的能够描述不同行为模式的样本作为训练集，但是又不能选取过多的训练集，因为采用过多的训练集会导致训练时间的增长；

(b)在输入样本中，代表每一样本的行为特征量需要尽量地能描述出该行为，因此选取的特征量不能太少；

(c)为了避免复杂运算，用来描述输入样本的特征量不能什么都包括，以致于训练后的模型过于敏感；

(2)疲劳驾驶状态的脑电信号功率谱分析：对处于不同状态下脑电信号的功率谱特性进行相关分析，将采集的脑电信号的长度选择为3秒，即对采集到的数据以每3秒截取一次并进行平均功率的计算，选取了从精力充沛状态过度到疲劳状态的脑电数据8分钟，将其截取为160段，并对截取到的这些时间段的脑电信号进行傅里叶变换：

$G=\underset{T\→\mathrm{\∞}}{\lim }\frac{1}{2T}\underset{-T}{\overset{T}{\∫}}{x}_T^2\left(t\right)dt=\frac{1}{2\mathrm{\π}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\underset{T\→\mathrm{\∞}}{\lim }\frac{1}{2T}{x}_T^2\left(f\right)df---\left(2\right)$

式中，f为连续频率变量；x_T(f)为信号的频谱密度；

其中：

$p_X\left(f\right)=\underset{T\→\mathrm{\∞}}{\lim }\frac{1}{2T}{x}_T^2\left(f\right)---\left(3\right)$

在式(3)中显示了处于不同频段上的信号功率的分布情况，同时也可以表示单位时间内和单位频段上的信号的能量，因此该积分函数又被称作为功率谱密度函数；经过傅里叶变换以后，得到x_n(e^jf)，然后对频谱的模进行平方，最后得到功率谱密度的计算公式：

$p_X\left(f\right)=\frac{1}{N}|X_N\left(f\right){|}^2---\left(4\right)$

由上面可知信号在某个频段h的平均功率密度为：

$G\left(h\right)=\frac{{\∫}_{f_l}^{f_h}{p}_X\left(f\right)df}{f_h-f_l}---\left(5\right)$

其中：f_h为频带h的上限，f_l为频带h的下限，p_X(f)为信号的功率谱密度；

则不同频段的平均功率比值则定义为：

$F=\frac{G\left(h_1\right)}{G\left(h_2\right)}---\left(6\right)$

式中h₁和h₂为脑电信号的不同频带，将F_(α+θ)/β,F_α/β这两个相关节律波的平均功率的比值选作进行驾驶疲劳判别的两个特征指标；

(3)疲劳驾驶状态的脑电信号样本熵分析：利用样本熵对脑电信号进行分析的具体步骤如下：计算选取长度为1000的数据作为一个数据窗口，以500个数据为移动步长，按以上步骤计算每次移动的数据的样本熵值，然后进行下1个数据窗口的样本熵值的计算，取2个数据窗口的样本熵的平均值作为时间序列的脑电数据样本熵值，通过对处于疲劳与清醒状态的脑电数据进行分析，得到不同状态下脑电信号的样本熵值；非疲劳状态时的样本熵值大于疲劳状态下的样本熵值，表明样本熵可以用于区分不同的驾驶状态，进一步可以得出疲劳状态下脑电信号计算的到样本熵值主要集中在0.2～0.4，而清醒状态的样本熵值集中在0.5～0.8之间；由此可知随着驾驶疲劳成都的加深脑电信号的样本熵值逐渐降低的，因此可看出脑电信号的样本熵值与人体精神状态之间存在相互对应关系，因此本发明将其选作为用于驾驶疲劳识别的相关特征量；

(4)疲劳驾驶状态的脑电信号Kc复杂度分析：首先采集得到的一段待分析的脑电时间序列，设为S(S₁,S₂,...,S_n)，第一步是对该序列进行粗粒化处理，具体就是先求取该序列的平均值，然后对原序列的每个元素进行判断并重新赋值，如果序列元素S_i大于平均值时则赋值S′₁＝1，如果序列元素S_i小于平均值时则赋值S_i'＝0，然后得到重构序列S'(S′₁,S'₂,...,S'_n)，重新利用S′_i重新对S中的元素进行赋值，则得到新的序列S(S₁,S₂,...,S_n)，该新序列是一串由(0，1)字符组成的新序列；下一步则对重构序列进行子串划分与界定，将按照如下的规则进行：

在一个字符串S(S₁,S₂,...,S_n)后再加一个字符串Q(q₁,q₂,...,q_n)，这样就得到组合后新的字符串序列SQ，令SQ序列删去末尾的一个字符所得字符串命名为SQv，接下来需要对Q是否属于SQv的子串进行判断，如果属于，则把下一个字符加到Q后面，增长一位再进行重复判断，直到增长后的Q不再属于SQv的一个子串，则用标记符号“.”进行序列的前后界定；下一步将“.”前的所有字符定义为新的S，重新以其后的字符开始构造Q，重复以上的判断和界定过程，直到整个字符串的序列判断结束；统计最终由“.”界定划分出的S的子串数目就是序列的Kc复杂度；

由于Kc复杂度的计算对数据的长度没有特定的要求，但是数据长度过长将影响计算的效率，而数据长度过短将不具有代表性，通过多次实验，最后本发明选择数据长度为1000，得到的不同状态下脑电信号的的Kc复杂度值结果；处于不同驾驶状态的Kc复杂度值存在着一定的差异性，其中非疲劳驾驶状态的Kc复杂度值主要集中在0.15～0.23，而处于疲劳驾驶状态得到的Kc复杂度值主要集中在0.08～0.14左右，说明受试者在非疲劳状态下进行驾驶操作的大脑复杂性高于疲劳状态下进行驾驶操作的大脑复杂性，因此表明Kc复杂度可以作为驾驶疲劳检测的另一个指标；

(5)支持向量机：考虑样本其中x_i代表第i个输入样本，y_i则是对应输入的输出；首先假设由子集y_i＝+1代表的类和y_i＝-1代表的类是能够进行线性可分的，那么用于进行分离的超平面形式的决策曲面方程可以表示为如下：

w^Tx+b＝0 (7)

x——输入变量；

w——可调的权值向量；

b——偏值；

那么式(7)则可以写成：

w^Tx+b≥0,当y_i＝+1 (8)

w^Tx+b≤0,当y_i＝-1 (9)

由式(8)和(9)则可合并为：

y_i(w^Tx+b)≥0 (10)

对于一个给定的权值向量w和偏置b，数据点到(1)定义的超平面的最近间隔称为分离边缘，用ρ表示，而支持向量机的目标则是寻找一个使得分离边缘ρ最大的超平面，而这个决策曲面则被称为最优超平面；

假设w₀和b₀分别为得到的权值向量和偏置的最优值；那么此时，输入空间里表示的多维线性决策面的最优超平面则可以表示如下：

$w_0^{T}x+b_0=0---\left(11\right)$

假如x_k对应于y_k，当满足时，那么则称(x_k,y_k)是支持向量的，那么此时支持向量到最优超平面的距离为：

$r=\left\{\begin{array}{cc}\frac{1}{\left|\right|w_0\left|\right|}& y_k=+1\\ -\frac{1}{\left|\right|w_0\left|\right|}& y_k=-1\end{array}\right.---\left(12\right)$

其中的加号表示x_k在最优超平面的正面，负号则表示最优超平面的负面，式(12)可以得到分类间隔r：

$\mathrm{\ρ}=2r=\frac{2}{\left|\right|w_0\left|\right|}---\left(13\right)$

由上式得到的最优解则对应于分类间隔最大的情况，其等价于凸二次最优化问题：

$\left\{\begin{array}{c}min\frac{1}{2}{w}^T\·w\\ s.t.:y_i(w^T{x}_i+b)\≥1,i=1,...,m\end{array}\right.---\left(14\right)$

由上式得到的最优解则是由拉格朗日(Lagrange)函数的鞍点决定：

$L(w,b,a)=\frac{1}{2}\left|\right|w|{|}^2-\underset{i=1}{\overset{N}{\Σ}}{a}_i\[y_i(w^T{x}_i+b)-1\]---\left(15\right)$

式中，

a_i,a_j——拉格朗日乘子；

由最优解条件(KKT条件)可得：

$\left.\begin{array}{c}\frac{\∂L(w,b,a)}{\∂b}=0\⇒\underset{i=1}{\overset{m}{\Σ}}{a}_i{y}_i=0\\ \frac{\∂L(w,b,a)}{\∂w}=0\⇒w=\underset{i=1}{\overset{m}{\Σ}}{a}_i{x}_i{y}_i\end{array}\right\}---\left(16\right)$

根据KKT(Karush-Kuhn-Tucker)条件，相应的拉格朗日乘子a_i应满足如下条件：

a_i[y_i(w^Tx_i+b)-1]＝0 (17)

从上式可以看出只有支持向量的a_i不为零；

根据核函数变换，则式(15)和(17)可转化为对偶问题求解，即：

$\left\{\begin{array}{c}W\left(a\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{a}_i{a}_j{y}_i{y}_j{x}_i^T{x}_j\\ \begin{array}{cc}st.\underset{i=1}{\overset{N}{\Σ}}{y}_i{a}_i=0& a_i\≥0,i=1,...,N\end{array}\end{array}\right.---\left(18\right)$

然后由(16)可以求出最优权值向量w₀，以及利用式子(17)可以得到b₀，于是关于待测样本x，通过只算sgn[(w^Tx)+b]就可以判断其所属的类别；

但是实际应用当中的大部分问题都是属于线性不可分的情况，而造成这样的原因也是各种各样的，可能在样本的产生时造成一定误差，也可能在进行统计记录，标记的错误导致样本数据的误差等，这样也可能使的样本出现线性不可分的情况；为了处理当存在少量样本无法用线性方法分开的情况，因此在分类超平面的定义中引入松弛变量ξ_i≥0，即允许少量误差的存在，那么(10)的约束条件变为如下形式：

y_i(w^Tx+b)≥1-ξ_i,i＝1,2,...N (19)

此时的目标函数则变为：

$\left\{\begin{array}{c}min\frac{1}{2}{w}^T\·w+C\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i\\ y_i(w^{T}x+b)\≥1-{\mathrm{\ξ}}_i,i=1,2,...N\end{array}\right.---\left(20\right)$

其中C为惩罚因子且C≥0，当选取的参数C比较大时，表示支持向量机模型对训练样本的质量具有较高的信心，相反，当选取的参数C比较小时，则认为训练样本中存在噪声；由(15)的到相应的式子如下所示：

$L(w,b,a)=\frac{1}{2}\left|\right|w|{|}^2+C\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i-\underset{i=1}{\overset{N}{\Σ}}{a}_i\[y_i(w^T{x}_i+b)+{\mathrm{\ξ}}_i-1\]-\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\β}}_i{\mathrm{\ξ}}_i---\left(16\right)$

那么根据KKT条件转化为其对偶问题如下:

$W\left(a\right)=\underset{i=1}{\overset{m}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{m}{\Σ}}{a}_i{a}_j{y}_i{y}_j{x}_i^T{x}_j---\left(17\right)$

$st.\underset{i=1}{\overset{m}{\Σ}}{y}_i{a}_i=0,0\≤a_i\≤C,i=1,...,m---\left(18\right)$

将输入空间Rⁿ中的点利用非线性映射Φ:Rⁿ→H映射到高维的Hilbert空间当中去，然后再使用线性可分的方法进行最优分类面的求解，而在这个求解的过程中，仅仅涉及到内积的运算，并且计算量也只依赖于训练样本数的个数，因此这个过程的就是如何找到一个核函数K(x_i,x_j),x_i,x_j∈Rⁿ，使得该核函数满足K(x_i,x_j)＝Φ^T(x_i)Φ(x_j)，这样就可以保证即使在不知道变换函数Φ形式的情况下也能够求得最终的结果，同时能够很好地避免维数灾难的发生，因此在最大间隔情况下，非线性支持向量机的目标函数为：

$W\left(a\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{a}_i{a}_j{y}_i{y}_{j}K(x_i,x_j)---\left(19\right)$

同样可得知，非线性支持向量机的最优化问题如下所示：

$\left\{\begin{array}{c}min\frac{1}{2}{w}^T\·w+C\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i\\ st:\underset{i=1}{\overset{m}{\Σ}}{y}_i{a}_i=0,0\≤a_i\≤C,i=1,...,m\end{array}\right.---\left(20\right)$

它的对偶问题则为：

$\left\{\begin{array}{c}\max W\left(a\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{a}_i{a}_j{y}_i{y}_{j}K(x_i,x_j)\\ st:\underset{i=1}{\overset{N}{\Σ}}{y}_i{a}_i=0,0\≤a_i\≤C,i=1,...,N\end{array}\right.---\left(21\right)$

相应得到的分类函数为：

$f\left(x\right)=\mathrm{sgn}\[\left(w^T\mathrm{\Φ}\right(x\left)\right)+b\]=\mathrm{sgn}\[\underset{i=1}{\overset{N}{\Σ}}{y}_i{a}_{i}K(x_i,x_j)+b\]---\left(22\right)$

(6)最小二乘支持向量机(LS-SVM)：最小二乘支持向量机算法的描述如下：首先选择一个非线性的变换把n维输入和1维输出的样本向量(x₁,y₁)，(x₂,y₂)，…，(x_N,y_N)，x_i∈Rⁿ,y_i∈R，i＝1,2,...,N，l从原空间映射到高维特征空间，则分类LS-SVM的优化问题则变为：

即将支持向量机的不等式约束改为等式约束，构造拉格朗日函数进行求解：

$L(w,b,a)=\frac{1}{2}\left|\right|w|{|}^2+\frac{C}{2}\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i^2-\underset{i=1}{\overset{N}{\Σ}}{a}_i\[y_i(w^T{x}_i+b)-1+{\mathrm{\ξ}}_i\]---\left(24\right)$

其中a_i——拉格朗日乘子；

最优解满足KTT条件则：进而可以得到下面的方程组：

对于i＝1,2,...,N，通过式子(25)消去w和ξ得到如下线性系统：

$\left[\begin{array}{cc}0& Y^T\\ Y& {\mathrm{ZZ}}^T+C^{-1}I\end{array}\right]\left[\begin{array}{c}b\\ a\end{array}\right]=\left[\begin{array}{c}0\\ 1\end{array}\right]$

式中：

令Ω＝ZZ^T+C^-1I，求解可得到：

$\left\{\begin{array}{c}a={\mathrm{\Ω}}^{-1}(Y-b1)\\ b=\frac{1^T{\mathrm{\Ω}}^{-1}Y}{1^T{\mathrm{\Ω}}^{-1}1}\end{array}\right.---\left(27\right)$

最终得到的决策函数如下：

$f\left(x\right)=\underset{i=1}{\overset{N}{\Σ}}{a}_{i}K(x,x_i)+b;---\left(28\right)$

(7)粒子群优化(PSO)算法设置LS-SVM的模型训练参数C和g：采用最小二乘法支持向量机进行模型建立的过程来需要确定两个参数，即惩罚参数C和核函数参数g，不同参数对模型性能以及泛化能力有不同的影响；因此，在建模的时候，选取适当的参数能够极大地改善识别的准确度，接下来具体研究LS-SVM的参数影响：惩罚参数C可以理解为对超出样本数据的一个管制制度，他直接影响模型的复杂程度和稳定性，惩罚参数C能够使模型在复杂度和训练误差之间做出一个平衡，使模型的推广能力得到保证，若惩罚参数C过小，则对超出的样本数据惩罚就小，训练的误差就会变大，模型的推广能力就会变差，当使用模型对新的数据进行学习时，错分率就会很高，而当惩罚参数C过大时，训练学习精度相应提高，但会对训练样本产生过度的依赖，对新数据的认识度不够，导致错分率增加，也会影响模型的推广能力，此外，惩罚参数C值也会对算法处理样本中的离群点(数据集合中一些元素和其他元素很不相似或不一致)产生影响，选取合适的值就能在一定程度上抗干扰，使模型稳定；

核函数参数g对支持向量的性能有很大的影响，对于高斯核中的核宽度δ而言，它反映了支持向量之间的相关程度，δ过小，支持向量之间的联系比较松弛，机器学习相对复杂，推广能力得不到保证，δ过大，支持向量之间的影响过强，预测模型很难达到足够的精度，所以在选取参数时，应综合考虑他们之间的相互影响；

因此，只有恰当的选择LS-SVM的参数才能得到最小的经验误差，同时提高模型的训练速度以及模型的泛化能力，使得预测结果更为准确；其中传统的模型参数选取主要通过交叉验证法以及经验选择方法等，但是使用经验选择发通常需要获知该重逢的经验知识，在通过判断进行参数的设置，要求比较高，因此一般不采用，而交叉验证法较为普遍，是大多数所采用的方法，便于实现，步骤简单，但是通过这种方法获得的参数精确度一般，并且需要花费较长的时间，本发明采用智能优化算法中的粒子群算法进行对模型参数的寻优，使其获得较好的性能，粒子群算法具体如下：首先在可解空间内初始化一群粒子，而每一个粒子都同时拥有位置、速度和适应度值三个特征，而其中的适应度值的大小则决定该粒子的优劣，同时粒子的移动主要是通过跟踪个体极值Pbest、群体极值Gbest以及过去速度的惯性来更新个体位置，其中个体极值Pbest是指单个粒子所经历的适应度值最优的位置，群体极值Gbest是指种群中所有粒子所经历适应度值最优的位置，粒子每迭代一次就移动一次位置，其适应度值就相应的需要重新计算一次，然后通过对新粒子个体极值的适应度值、群体极值的适应度值与迭代之前的极值比较来确定是否更新粒子的个体极值和群体极值的位置；

假设在一个D维的空间中，其中X＝(X₁,X₂,...,X_n)为由n个粒子组成的种群，其中向量X_i＝(x_i1,x_i2,...,x_iD)^T表示第i个粒子，对于第i个粒子，它的适应度值一般就是根据粒子X_i在目标函数中的值，速度为V_i＝(v_i1,v_i2,...,v_iD)^T，其个体极值为P_i＝(p_i1,p_i2,...,p_iD)^T，种群的全局极值为P_g＝(p_g1,p_g2,...,p_gD)^T

在每一次的迭代过程中，每个粒子都是通过个体极值、全局极值以及惯性速度来更新它的速度和位置信息，其更新公式(29)和(30)如下所示：

$v_{id}^{k+1}={\mathrm{wv}}_{id}^k+c_1{r}_1(p_{id}^k-x_{id}^k)+c_2{r}_2(p_{gd}^k-x_{id}^k)---\left(29\right)$

$x_{id}^{k+1}=x_{id}^k+v_{id}^{k+1}---\left(30\right)$

式(24)中，w为惯性权重，当w较大时有利于群体在较大的范围内进行搜索，从而避免陷入寻找局部最优，而较小的w能够保证算法能够快速收敛到最优值；d＝1,2,...,D，i＝1,2,...,n，k为当前迭代次数；v_id为粒子的速度；c₁和c₂为非负常数，称为加速度因子，是每一个粒子在P_i和P_g位置的随机加速权重，r₁和r₂是[0，1]之间的随机数；为了防止粒子群的盲目搜索，一般把粒子群的位置限定在[-X_max,X_max]和速度限定在[-V_max,V_max]的范围内；

根据上面的分析，本发明使用粒子群算法优化参数C和g，具体步骤如下：

(a)初始化粒子群并随机产生所有粒子的位置和速度并确定粒子的P_i和P_g；

(b)将每个粒子当前最优位置P_i及种群中所有粒子所经历的最好位置进P_g行比较，如果这个粒子的位置最优，则将其设置为当前的最好位置P_g，否则，P_g保持不变；

(c)更新粒子的速度、位置；

(d)若没有达到结束条件，如：预设的运算精度、迭代次数，便返回步骤(b)；

(e)开始下一次的迭代计算，否则取当前的P_g为最优解。