人体疲劳程度检测，生理信号处理

写在前面 ：

参加比赛做了检测人体人体疲劳程度的手环，题主现在已经毕业。检测的生理指标有肌电、脉搏、血氧含量。对时域和频域进行分析，通过支持向量机分级模型得出疲劳等级。

疲劳分级

本系统根据斯坦福嗜睡量表分为四个疲劳等级：清醒(0级)、轻度疲劳(1级)、中度疲劳(2级)、重度疲劳(3级)

​ 表1 斯坦福嗜睡量表(SSS)

评分疲劳程度疲劳等级1感觉非常振奋、充满生机与活力0级2身体机能处于较高水平，不过不在峰值，可以集中注意力0级3很清醒，不过身体及思维都比较放松，能及时反应但不够灵敏1级4有些倦意、松懈1级5充满倦意，不再想保持清醒，非常松懈2级6开始打瞌睡，头晕眼花，不再与睡意作斗争，只想躺下来休息3级7睡眠初期，开始出现梦境3级X沉睡3级

检测指标

根据相关研究，疲劳状态通过选取2个或2个以上检测指标均认为有效，且现有研究中使用频率最高的检测指标主要为肌电、脉搏和血氧饱和度信号，考虑到数据采集的便捷性以及疲劳程度的相关性，本文选取了这3个生理指标作为驾驶疲劳的检测标准。

（1）肌电(EMG)特征参数提取

研究发现EMG信号随着人体疲劳程度的加深，其幅值会增大，频率会降低，因此

通过分析EMG信号来检测驾驶员的疲劳状态。通过对处理后的数据进行时域分析和频域分析，得到肌电信号的特征参数：积分肌电值(IEMG)、均方根值(RMS)、平均功率频率(MPF)、中位频率(MF)。

IEMG是指在一定时间内肌肉中参与活动的运动单位放电总量，通常其幅值随疲劳程度的加深而增加。RMS反映了运动单位数量的募集程度及其放电频率的同步化程度。其计算公式分别为：

​

I

E

M

G

=

∫

t

1

t

2

∣

X

(

t

)

∣

d

t

I E M G=\int_{t_{1}}^{t_{2}}|X(t)| d t

IEMG=∫t1​t2​​∣X(t)∣dt

R

M

S

=

1

T

∑

t

=

1

T

x

2

(

t

)

R M S=\sqrt{\frac{1}{T} \sum_{t=1}^{T} x^{2}(t)}

RMS=T1​t=1∑T​x2(t)

​ 其中，t1为采集肌电信号开始时刻，t2为采集肌电信号结束时刻，*X(t)*为肌电信号大小。

将表面肌电(EMG)信号进行快速傅里叶变换( FFT)，获得EMG信号的功率谱，反映 EMG信号在不同频率范围内的强度。MF和MPF为常用的肌肉疲劳的评价指标。其计算公式为：

M

P

F

=

∫

0

∞

f

⋅

P

S

D

(

f

)

d

f

∫

0

∞

P

S

D

(

f

)

d

f

M P F=\frac{\int_{0}^{\infty} f \cdot P S D(f) d f}{\int_{0}^{\infty} P S D(f) d f}

MPF=∫0∞​PSD(f)df∫0∞​f⋅PSD(f)df​ ​

∫

0

M

F

P

S

D

(

f

)

d

f

=

∫

M

F

f

0

P

S

D

(

f

)

d

f

=

1

2

∫

0

f

0

P

S

D

(

f

)

d

f

\int_{0}^{M F} P S D(f) d f=\int_{M F}^{f_{0}} P S D(f) d f=\frac{1}{2} \int_{0}^{f_{0}} P S D(f) d f

∫0MF​PSD(f)df=∫MFf0​​PSD(f)df=21​∫0f0​​PSD(f)df 其中，PSD(f)为功率谱密度，MF为中位频率，f为肌电信号频率，f0为频率上限，即采样频率的一半，等于500 Hz，MPF与MF下降表明被试者出现了肌肉疲劳。

（2）脉搏特征参数提取

通过对处理后的脉搏信号数据进行时域分析和频域分析，得到脉搏信号的特征参数：

RR间期均值(MEAN)、RR间期标准差(SDNN)、低频功率(LF)、高频功率(HF)。在脉搏信号中，首先运用小波变换检测出各主波的波峰，求出两个相邻主波波峰之间的时间差，记作RRi(i=1,2,3……)，则RRi的均值即为RR间期均值、RRi的标准差即为RR间期标准差，并运用快速傅里叶变换(FFT)将时域指标变为频域指标，得到LF和HF。其计算公式为：

M

E

A

N

=

∑

i

=

1

N

R

R

i

N

M E A N=\sum_{i=1}^{N} \frac{R R_{i}}{N}

MEAN=i=1∑N​NRRi​​

S

D

N

N

=

1

N

∑

i

=

1

N

(

R

R

i

−

M

E

A

N

)

2

S D N N=\sqrt{\frac{1}{N} \sum_{i=1}^{N}\left(R R_{i}-M E A N\right)^{2}}

SDNN=N1​i=1∑N​(RRi​−MEAN)2

​ 其中，N表示正常脉搏搏动的总次数，RRi表示第i个相邻主波波峰之间的时间差，MEAN表示RR间期均值，SDNN表示RR间期标准差。

（3）血氧饱和度(SpO2)特征参数提取

当人体处于疲劳状态时，SpO2含量有所下降[9]。通过对处理后的实验数据进行时域分析，得到SpO2的特征参数：血氧标准差(RSD)。根据所得结果的波动程度，便可准确的反映出被试者的疲劳状态。其计算公式为：

R

S

D

=

1

N

−

1

∑

i

=

1

N

(

R

i

−

R

A

V

G

)

2

R S D=\sqrt{\frac{1}{N-1} \sum_{i=1}^{N}\left(R_{i}-R A V G\right)^{2}}

RSD=N−11​i=1∑N​(Ri​−RAVG)2

​ 其中，RAVG表示血氧饱和度平均值，N表示采集样本中血氧数据总个数，Ri表示第i个样本的血氧饱和度值。

程序运行-matlab

程序采用MATLAB运行，运行时候，需 更改文件路径，滤波器用于滤除采集的生理信号的噪声。脉搏、肌电信号在采集过程中均为低频微弱信号，在实验中常会混入各种噪声，包括基线漂移，50Hz工频干扰以及周围环境的噪声干扰等，故采用巴特沃斯滤波器处理信号中的噪声。血氧信号在采集过程中，由于其频率、幅度变化比较稳定，故一般不需对其进行数据预处理

原数据文件：百度网盘 提取码：mmb7

巴特沃斯滤波器

%带通滤波

%x-输入信号

%fs-输入信号采样频率

%M-滤波器半阶数

%fl-下截止频率

%fh-上截止频率

%xl-返回的滤波结果,其长度为原序列长度减少2M,去掉了前M和后M个不正确的点

function xl=beltfilter(x,fs,M,fl,fh)

N=length(x)-2*M;

[i j]=size(x);

if i~=1

x=x';

end

k=1:M;

w=0.5+0.5*cos(pi*k/M);

wl=2*pi*fl/fs;

wh=2*pi*fh/fs;

h(1)=(wh-wl)/pi;

h(2:M+1)=(sin(wh*k)-sin(wl*k))./(pi*k).*w;

for k=1:N

kk=k-1+M;

xl(k)=x(kk+1)*h(1)+sum(h(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));

end

if i~=1

xl=xl';

end

积分肌电值(IEMG)、均方根值(RMS)、平均功率频率(MPF)、中位频率(MF)

clc;

fs=1000;

%%%==================批量读入文件=======================

path = 'C:\Users\hp\Desktop\EMG\数据\';

namelist = dir([path,'*.txt']);

l = length(namelist);

P = cell(1,1);%定义一个细胞数组，用于存放所有txt文件

for i = 1:l

namelist(i).name;%获得该路径下的文件名，如1.txt是相对路径

filename{i} = [path,namelist(i).name];%通过字符串拼接获得的就是绝对路径了

s = load(filename{i});

x=s(:,1);

%%=================批量处理数据======================

x=x/1000;

figure(1);

subplot(311),plot(x);

title("原信号波形");

y=fft(x); %对信号进行快速Fourier变换

mag=abs(y); %求得Fourier变换后的振幅

subplot(3,1,2) ,plot(mag);

axis([0 500 0 1.1*max(mag)]);

title("傅里叶变换");

grid on;

x=beltfilter(x,fs,80,0,250);

L=length(x);f=(0:L-1)*fs/L;

y=fft(x); %对信号进行快速Fourier变换

mag=abs(y); %求得Fourier变换后的振幅

subplot(3,1,3) ,plot(f,mag);

axis([0 500 0 1.1*max(mag)]);

title("带通滤波器");

grid on;

iemg=sum(abs(x))/length(x); %求积分肌电值

rms=sqrt(sum(x.^2)/length(x)); %求均方根值

L=length(x);f=(0:L-1)*fs/L;

cx=xcorr(x,'unbiased');

cxk=fft(cx,L);

px=abs(cxk);%求功率谱密度

pxx=10*log10(px);

f=(0:L-1)*fs/L;

grid on %做功率谱图

df=fs/L;

p=(sum(px(1:L/2-1))+sum(px(1:L/2)))/2.*df;

pf=(sum(px(1:L/2-1).*[1:L/2-1]'.*df)+sum(px(1:L/2).*[1:L/2]'.*df))/2*df;

MPF=pf/p; %求平均功率频率

N=1;pp=0;

while abs(pp-p/2)>(px(N)+px(N+1))/2*df

pp=pp+(px(N)+px(N+1))/2*df;

N=N+1;

end

n_1=(N+N+1)/2;

MF=df*n_1; %求中值频率

fprintf("\n积分肌电值 %.4f\n",iemg);

fprintf("均分根值 %.4f\n",rms);

fprintf("平均功率频率 %.4f\n",MPF);

fprintf("中值频率 %.4f\n\n",MF);

%%=================记录数据=====================

result1(i)=iemg;

result2(i)=rms;

result3(i)=MPF;

result4(i)=MF;

end

%%%==================数据文件的写入=======================

fid = fopen('C:\Users\hp\Desktop\EMG\数据\结果\iemg.txt','wt');

fprintf(fid,'%g\n',result1);

fclose(fid);

fid = fopen('C:\Users\hp\Desktop\EMG\数据\结果\rms.txt','wt');

fprintf(fid,'%g\n',result2);

fclose(fid);

fid = fopen('C:\Users\hp\Desktop\EMG\数据\结果\MPF.txt','wt');

fprintf(fid,'%g\n',result3);

fclose(fid);

fid = fopen('C:\Users\hp\Desktop\EMG\数据\结果\MF.txt','wt');

fprintf(fid,'%g\n',result4);

fclose(fid);

%%======================绘图===========================

figure(2)

n=1:l;

subplot(411)

plot(n,result1)

axis([0 1.1*l 0 1.1*max(result1)]);

%suptitle('肌电的4种指标')

title("肌电积分值(IEMG)");

xlabel("时间(min)");

subplot(412)

plot(n,result2)

axis([0 1.1*l 0 1.1*max(result2)]);

title("均方根值(RMS)");

xlabel("时间(min)");

subplot(413)

plot(n,result3)

axis([0 1.1*l 0 1.1*max(result3)]);

title("平均功率频度(MPF)");

xlabel("时间(min)");

subplot(414)

plot(n,result4)

axis([0 1.1*l 0 1.1*max(result4)]);

title("中位频率(MF)");

xlabel("时间(min)");

血氧均值、血氧均方差（RSD）

clc;

%%%==================批量读入文件=======================

path = 'C:\Users\hp\Desktop\SPO2\数据\';

namelist = dir([path,'*.txt']);

l = length(namelist);

P = cell(1,1);%定义一个细胞数组，用于存放所有txt文件

for i = 1:l

namelist(i).name;%获得该路径下的文件名，如1.txt是相对路径

filename{i} = [path,namelist(i).name];%通过字符串拼接获得的就是绝对路径了

s = load(filename{i});

x=s(:,1);

figure(1)

plot(x);

title("SPO2数据");

axis([0 100 0 1.1*max(x)]);

%%=================批量处理数据======================

L=length(x);

RAVG=sum(x)/L; %%血氧均值

RSD=sqrt( sum((x-RAVG).^2)/L ); %%血氧均方差

fprintf("\n 血氧均方差: %.4f\n\n",RSD);

result(i)=RSD;

RAVG1(i)=RAVG;

end

%%%==================数据文件的写入=======================

fid = fopen('C:\Users\hp\Desktop\SPO2\数据\结果\result.txt','wt');

fprintf(fid,'%g\n',result);

fclose(fid);

%%======================绘图===========================

figure(2)

n=1:l;

subplot(211)

plot(n,RAVG1)

title("血氧平均值");

grid on;

% set(gca,'xminortick','on'); %坐标刻度显示

% set(gca,'ticklength',[0.02 0.025]);

% set(gca,'tickdir','out');

xlabel("时间(min)");

ylabel("含量(%)");

subplot(212)

plot(n,result)

title("血氧均方差")

axis([0 1.1*l 0 1.1*max(result)]);

xlabel("时间(min)");

grid on;

实验设备

表2 实验设备简介

设备名称数量功能简介SEN0240肌电传感器1个采集驾驶员的肌电信号MAX30102脉搏血氧传感器1个采集驾驶员的脉搏信号和血氧信号STM32F103单片机2个设备的控制和数据处理单元HC-05蓝牙模块2个实现手环系统和车外警示系统的连接与通信RGB七彩LED灯1个依据不同颜色，显示出驾驶员当前所处的疲劳状态车外警示灯1个给予道路上其他驾驶员警示作用振动模块1个通过振动，提醒驾驶员应注意休息语音模块1个通过语音提示，提醒驾驶员应注意休息