[发明专利]一种基于脑电信号的驾驶疲劳的检测方法

摘要

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

主权项

一种基于脑电信号的驾驶疲劳的检测方法，其特征在于：具体包括以下方面：(1)选取与疲劳驾驶状态相关的特征量：特征选择是指通过对原始测量信号进行一系列相应处理后，得到一组由相关特征值组成的特征集，然后从该特征集中选取某些具有代表性、分类性能最佳的物理量作为区分不同行为的一组特征量，然后使用相关的分类识别方法，利用选取的特征量将不同的行为从特征空间中分离出来，因此对于学习分类而言都需要训练集样本：T＝{(x1,y1),…,(xl,yl)}∈(X×Y)l   (1)其中xi∈X＝Rn,yi∈Y＝R,i＝1,…,l；对于符合要求的训练集就是如何选择满足要求的特征量，具体是：(a)需要选取适当数量的能够描述不同行为模式的样本作为训练集，但是又不能选取过多的训练集，因为采用过多的训练集会导致训练时间的增长；(b)在输入样本中，代表每一样本的行为特征量需要尽量地能描述出该行为，因此选取的特征量不能太少；(c)为了避免复杂运算，用来描述输入样本的特征量不能什么都包括，以致于训练后的模型过于敏感；(2)疲劳驾驶状态的脑电信号功率谱分析：对处于不同状态下脑电信号的功率谱特性进行相关分析，将采集的脑电信号的长度选择为3秒，即对采集到的数据以每3秒截取一次并进行平均功率的计算，选取了从精力充沛状态过度到疲劳状态的脑电数据8分钟，将其截取为160段，并对截取到的这些时间段的脑电信号进行傅里叶变换：G=limT→∞12T∫-TTxT2(t)dt=12π∫-∞∞limT→∞12TxT2(f)df---(2)]]>式中，f为连续频率变量；xT(f)为信号的频谱密度；其中：pX(f)=limT→∞12TxT2(f)---(3)]]>在式(3)中显示了处于不同频段上的信号功率的分布情况，同时也可以表示单位时间内和单位频段上的信号的能量，因此该积分函数又被称作为功率谱密度函数；经过傅里叶变换以后，得到xn(ejf)，然后对频谱的模进行平方，最后得到功率谱密度的计算公式：pX(f)=1N|XN(f)|2---(4)]]>由上面可知信号在某个频段h的平均功率密度为：G(h)=∫flfhpX(f)dffh-fl---(5)]]>其中：fh为频带h的上限，fl为频带h的下限，pX(f)为信号的功率谱密度；则不同频段的平均功率比值则定义为：F=G(h1)G(h2)---(6)]]>式中h1和h2为脑电信号的不同频带，将F(α+θ)/β,Fα/β这两个相关节律波的平均功率的比值选作进行驾驶疲劳判别的两个特征指标；(3)疲劳驾驶状态的脑电信号样本熵分析：利用样本熵对脑电信号进行分析的具体步骤如下：计算选取长度为1000的数据作为一个数据窗口，以500个数据为移动步长，按以上步骤计算每次移动的数据的样本熵值，然后进行下1个数据窗口的样本熵值的计算，取2个数据窗口的样本熵的平均值作为时间序列的脑电数据样本熵值，通过对处于疲劳与清醒状态的脑电数据进行分析，得到不同状态下脑电信号的样本熵值；非疲劳状态时的样本熵值大于疲劳状态下的样本熵值，表明样本熵可以用于区分不同的驾驶状态，进一步可以得出疲劳状态下脑电信号计算的到样本熵值主要集中在0.2～0.4，而清醒状态的样本熵值集中在0.5～0.8之间；由此可知随着驾驶疲劳成都的加深脑电信号的样本熵值逐渐降低的，因此可看出脑电信号的样本熵值与人体精神状态之间存在相互对应关系，因此本发明将其选作为用于驾驶疲劳识别的相关特征量；(4)疲劳驾驶状态的脑电信号Kc复杂度分析：首先采集得到的一段待分析的脑电时间序列，设为S(S1,S2,...,Sn)，第一步是对该序列进行粗粒化处理，具体就是先求取该序列的平均值，然后对原序列的每个元素进行判断并重新赋值，如果序列元素Si大于平均值时则赋值S′1＝1，如果序列元素Si小于平均值时则赋值Si'＝0，然后得到重构序列S'(S′1,S'2,...,S'n)，重新利用S′i重新对S中的元素进行赋值，则得到新的序列S(S1,S2,...,Sn)，该新序列是一串由(0，1)字符组成的新序列；下一步则对重构序列进行子串划分与界定，将按照如下的规则进行：在一个字符串S(S1,S2,...,Sn)后再加一个字符串Q(q1,q2,...,qn)，这样就得到组合后新的字符串序列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)支持向量机：考虑样本其中xi代表第i个输入样本，yi则是对应输入的输出；首先假设由子集yi＝+1代表的类和yi＝‑1代表的类是能够进行线性可分的，那么用于进行分离的超平面形式的决策曲面方程可以表示为如下：wTx+b＝0   (7)x——输入变量；w——可调的权值向量；b——偏值；那么式(7)则可以写成：wTx+b≥0,当yi＝+1   (8)wTx+b≤0,当yi＝‑1   (9)由式(8)和(9)则可合并为：yi(wTx+b)≥0   (10)对于一个给定的权值向量w和偏置b，数据点到(1)定义的超平面的最近间隔称为分离边缘，用ρ表示，而支持向量机的目标则是寻找一个使得分离边缘ρ最大的超平面，而这个决策曲面则被称为最优超平面；假设w0和b0分别为得到的权值向量和偏置的最优值；那么此时，输入空间里表示的多维线性决策面的最优超平面则可以表示如下：w0Tx+b0=0---(11)]]>假如xk对应于yk，当满足时，那么则称(xk,yk)是支持向量的，那么此时支持向量到最优超平面的距离为：r=1||w0||yk=+1-1||w0||yk=-1---(12)]]>其中的加号表示xk在最优超平面的正面，负号则表示最优超平面的负面，式(12)可以得到分类间隔r：ρ=2r=2||w0||---(13)]]>由上式得到的最优解则对应于分类间隔最大的情况，其等价于凸二次最优化问题：min12wT·ws.t.:yi(wTxi+b)≥1,i=1,...,m---(14)]]>由上式得到的最优解则是由拉格朗日(Lagrange)函数的鞍点决定：L(w,b,a)=12||w||2-Σi=1Nai[yi(wTxi+b)-1]---(15)]]>式中，ai,aj——拉格朗日乘子；由最优解条件(KKT条件)可得：∂L(w,b,a)∂b=0⇒Σi=1maiyi=0∂L(w,b,a)∂w=0⇒w=Σi=1maixiyi---(16)]]>根据KKT(Karush‑Kuhn‑Tucker)条件，相应的拉格朗日乘子ai应满足如下条件： ai[yi(wTxi+b)‑1]＝0   (17)从上式可以看出只有支持向量的ai不为零；根据核函数变换，则式(15)和(17)可转化为对偶问题求解，即：W(a)=Σi=1Nai-12Σi,j=1NaiajyiyjxiTxjst.Σi=1Nyiai=0ai≥0,i=1,...,N---(18)]]>然后由(16)可以求出最优权值向量w0，以及利用式子(17)可以得到b0，于是关于待测样本x，通过只算sgn[(wTx)+b]就可以判断其所属的类别；但是实际应用当中的大部分问题都是属于线性不可分的情况，而造成这样的原因也是各种各样的，可能在样本的产生时造成一定误差，也可能在进行统计记录，标记的错误导致样本数据的误差等，这样也可能使的样本出现线性不可分的情况；为了处理当存在少量样本无法用线性方法分开的情况，因此在分类超平面的定义中引入松弛变量ξi≥0，即允许少量误差的存在，那么(10)的约束条件变为如下形式：yi(wTx+b)≥1‑ξi,i＝1,2,...N   (19)此时的目标函数则变为：min12wT·w+CΣi=1Nξiyi(wTx+b)≥1-ξi,i=1,2,...N---(20)]]>其中C为惩罚因子且C≥0，当选取的参数C比较大时，表示支持向量机模型对训练样本的质量具有较高的信心，相反，当选取的参数C比较小时，则认为训练样本中存在噪声；由(15)的到相应的式子如下所示：L(w,b,a)=12||w||2+CΣi=1Nξi-Σi=1Nai[yi(wTxi+b)+ξi-1]-Σi=1Nβiξi---(16)]]>那么根据KKT条件转化为其对偶问题如下:W(a)=Σi=1mai-12Σi,j=1maiajyiyjxiTxj---(17)]]>st.Σi=1myiai=0,0≤ai≤C,i=1,...,m---(18)]]>将输入空间Rn中的点利用非线性映射Φ:Rn→H映射到高维的Hilbert空间当中去，然后再使用线性可分的方法进行最优分类面的求解，而在这个求解的过程中，仅仅涉及到内积的运算，并且计算量也只依赖于训练样本数的个数，因此这个过程的就是如何找到一个核函数K(xi,xj),xi,xj∈Rn，使得该核函数满足K(xi,xj)＝ΦT(xi)Φ(xj)，这样就可以保证即使在不知道变换函数Φ形式的情况下也能够求得最终的结果，同时能够很好地避免维数灾难的发生，因此在最大间隔情况下，非线性支持向量机的目标函数为：W(a)=Σi=1Nai-12Σi,j=1NaiajyiyjK(xi,xj)---(19)]]>同样可得知，非线性支持向量机的最优化问题如下所示：min12wT·w+CΣi=1Nξist:Σi=1myiai=0,0≤ai≤C,i=1,...,m---(20)]]>它的对偶问题则为：maxW(a)=Σi=1Nai-12Σi,j=1NaiajyiyjK(xi,xj)st:Σi=1Nyiai=0,0≤ai≤C,i=1,...,N---(21)]]>相应得到的分类函数为：f(x)=sgn[(wTΦ(x))+b]=sgn[Σi=1NyiaiK(xi,xj)+b]---(22)]]>(6)最小二乘支持向量机(LS‑SVM)：最小二乘支持向量机算法的描述如下：首先选择一个非线性的变换把n维输入和1维输出的样本向量(x1,y1)，(x2,y2)，…，(xN,yN)，xi∈Rn,yi∈R，i＝1,2,...,N，l从原空间映射到高维特征空间，则分类LS‑SVM的优化问题则变为：即将支持向量机的不等式约束改为等式约束，构造拉格朗日函数进行求解：L(w,b,a)=12||w||2+C2Σi=1Nξi2-Σi=1Nai[yi(wTxi+b)-1+ξi]---(24)]]>其中ai——拉格朗日乘子；最优解满足KTT条件则：进而可以得到下面的方程组：对于i＝1,2,...,N，通过式子(25)消去w和ξ得到如下线性系统：0YTYZZT+C-1Iba=01]]>式中：令Ω＝ZZT+C‑1I，求解可得到：a=Ω-1(Y-b1)b=1TΩ-1Y1TΩ-11---(27)]]>最终得到的决策函数如下：f(x)=Σi=1NaiK(x,xi)+b;---(28)]]>(7)粒子群优化(PSO)算法设置LS‑SVM的模型训练参数C和g：采用最小二乘法支持向量机进行模型建立的过程来需要确定两个参数，即惩罚参数C和核函数参数g，不同参数对模型性能以及泛化能力有不同的影响；因此，在建模的时候，选取适当的参数能够极大地改善识别的准确度，接下来具体研究LS‑SVM的参数影响：惩罚参数C可以理解为对超出样本数据的一个管制制度，他直接影响模型的复杂程度和稳定性，惩罚参数C能够使模型在复杂度和训练误差之间做出一个平衡，使模型的推广能力得到保证，若惩罚参数C过小，则对超出的样本数据惩罚就小，训练的误差就会变大，模型的推广能力就会变差，当使用模型对新的数据进行学习时，错分率就会很高，而当惩罚参数C过大时，训练学习精度相应提高，但会对训练样本产生过度的依赖，对新数据的认识度不够，导致错分率增加，也会影响模型的推广能力，此外，惩罚参数C值也会对算法处理样本中的离群点(数据集合中一些元素和其他元素很不相似或不一致)产生影响，选取合适的值就能在一定程度上抗干扰，使模型稳定；核函数参数g对支持向量的性能有很大的影响，对于高斯核中的核宽度δ而言，它反映了支持向量之间的相关程度，δ过小，支持向量之间的联系比较松弛，机器学习相对复杂，推广能力得不到保证，δ过大，支持向量之间的影响过强，预测模型很难达到足够的精度，所以在选取参数时，应综合考虑他们之间的相互影响；因此，只有恰当的选择LS‑SVM的参数才能得到最小的经验误差，同时提高模型的训练速度以及模型的泛化能力，使得预测结果更为准确；其中传统的模型参数选取主要通过交叉验证法以及经验选择方法等，但是使用经验选择发通常需要获知该重逢的经验知识，在通过判断进行参数的设置，要求比较高，因此一般不采用，而交叉验证法较为普遍，是大多数所采用的方法，便于实现，步骤简单，但是通过这种方法获得的参数精确度一般，并且需要花费较长的时间，本发明采用智能优化算法中的粒子群算法进行对模型参数的寻优，使其获得较好的性能，粒子群算法具体如下：首先在可解空间内初始化一群粒子，而每一个粒子都同时拥有位置、速度和适应度值三个特征，而其中的适应度值的大小则决定该粒子的优劣，同时粒子的移动主要是通过跟踪个体极值Pbest、群体极值Gbest以及过去速度的惯性来更新个体位置，其中个体极值Pbest是指单个粒子所经历的适应度值最优的位置，群体极值Gbest是指种群中所有粒子所经历适应度值最优的位置，粒子每迭代一次就移动一次位置，其适应度值就相应的需要重新计算一次，然后通过对新粒子个体极值的适应度值、群体极值的适应度值与迭代之前的极值比较来确定是否更新粒子的个体极值和群体极值的位置；假设在一个D维的空间中，其中X＝(X1,X2,...,Xn)为由n个粒子组成的种群，其中向量Xi＝(xi1,xi2,...,xiD)T表示第i个粒子，对于第i个粒子，它的适应度值一般就是根据粒子Xi在目标函数中的值，速度为Vi＝(vi1,vi2,...,viD)T，其个体极值为Pi＝(pi1,pi2,...,piD)T，种群的全局极值为Pg＝(pg1,pg2,...,pgD)T在每一次的迭代过程中，每个粒子都是通过个体极值、全局极值以及惯性速度来更新它的速度和位置信息，其更新公式(29)和(30)如下所示：vidk+1=wvidk+c1r1(pidk-xidk)+c2r2(pgdk-xidk)---(29)]]>xidk+1=xidk+vidk+1---(30)]]>式(24)中，w为惯性权重，当w较大时有利于群体在较大的范围内进行搜索，从而避免陷入寻找局部最优，而较小的w能够保证算法能够快速收敛到最优值；d＝1,2,...,D，i＝1,2,...,n，k为当前迭代次数；vid为粒子的速度；c1和c2为非负常数，称为加速度因子，是每一个粒子在Pi和Pg位置的随机加速权重，r1和r2是[0，1]之间的随机数；为了防止粒子群的盲目搜索，一般把粒子群的位置限定在[‑Xmax,Xmax]和速度限定在[‑Vmax,Vmax]的范围内；根据上面的分析，本发明使用粒子群算法优化参数C和g，具体步骤如下：(a)初始化粒子群并随机产生所有粒子的位置和速度并确定粒子的Pi和Pg；(b)将每个粒子当前最优位置Pi及种群中所有粒子所经历的最好位置进Pg行比较，如果这个粒子的位置最优，则将其设置为当前的最好位置Pg，否则，Pg保持不变；(c)更新粒子的速度、位置；(d)若没有达到结束条件，如：预设的运算精度、迭代次数，便返回步骤(b)；(e)开始下一次的迭代计算，否则取当前的Pg为最优解。