论文+代码阅读：基于图像的EEG疲劳分类深度递归神经网络

2023/10/15

该文章主要记录该论文中提到的主要方法及其实现过程。

原文地址：Classification of EEG Signals from Driving Fatigue by Image-Based Deep Recurrent Neural Networks

参考印证论文地址：

Learning Representations from EEG with Deep Recurrent-Convolutional Neural Networks

参考印证论文代码库：

code

p.s 原文未给出代码，相关代码来自于一篇类似论文

---

摘要

        评估疲劳驾驶检测系统一直是研究人员的主要兴趣。随着深度学习算法和多模态技术的发展，越来越多的学者开始从多通道脑电图时间序列数据中寻找新的表示方法。在这里，我们建议将疲劳和正常状态脑电图数据转换为一系列拓扑频谱图像，并伴随眼睑闭合比 (PERCLOS) 标记输入到我们训练的图像分类深度循环卷积神经网络中，以从图像序列中找到新的数据鲁棒表示。我们使用了 2017 年上海交通大学发布的数据集 SEED-VIG 来帮助构建疲劳检测系统。该数据集包含 23 个subject的眼动信号和 16 个 EEG 电极收集的多模态信号。我们使用适当的频带来表征疲劳以形成颜色频带，基于方位等距投影和Clough-Tocher插值算法形成三通道图像，形成二维图像而不是信号作为模型的输入。该方法与传统的脑电分析的特征向量表示不同，旨在保持脑电的空间、频谱和时间结构，并找到最能表征疲劳的维度(频率)范围。在疲劳检测上获得了满意的分类结果(96%)。

---

解决问题

        该论文旨在从从多通道脑电时间序列数据中寻找新的表示方法。将疲劳和正常状态脑电图数据转换为一系列拓扑频谱图像，并伴随眼睑闭合比 (PERCLOS) 标记输入到训练过后的图像分类深度循环卷积神经网络中，以从图像序列中找到新的数据鲁棒表示。

---

工作流程

操作概览

        该论文旨在从从多通道脑电时间序列数据中寻找新的表示方法。将疲劳和正常状态脑电图数据转换为一系列拓扑频谱图像，并伴随眼睑闭合比 (PERCLOS) 标记输入到训练过后的图像分类深度循环卷积神经网络中，以从图像序列中找到新的数据鲁棒表示。

图1 参考论文流程

图2 原文流程

1.图像生成

        受试者的数据被过滤成三个波段，然后对窗口等于0.5秒的时间序列应用快速傅里叶变换(FFT)操作，并使用16帧构建图像(眼动仪每8秒标记一次)PERCLSOS指数)，然后计算每个波段的平方绝对值之和，对于3波段场景(16帧 * 3波段 * 16通道=768)，我们得到了885片(885 * 8=7080)，由16个二维图像电极转换，从而将脑电信号分类的工作转化为图像分类。

2.神经网络

        在卷积操作期间，每个输入图像都与多个 2D 滤波器进行卷积，然后大小变大，然后通过最大池化层将其采样回原始大小。通过反向传播学习神经网络全时间和滤波器参数来减少分类错误。

模型细节

        具有高度可扩展架构的 VGG 网络（Simonyan & Zisserman，2015）用于卷积操作。卷积层的感受野设置为 3*3 的小感受野，步幅为 1 像素，具有 ReLU 激活函数。在多层卷积之后，进入池化层，池化窗口设置为 2*2，步幅设置为 2 像素。通过改变堆叠卷积层的数量、过滤器的总数和大小、损失函数、下降率、学习因子等参数：3 个堆叠卷积层（3-64），然后是最大池化，然后是 2 1 个卷积层（3-128），然后是另一个最大池化层，然后是 2 个卷积层（3-256），然后是 ReLU 层，最后是最大池化层。使用随机梯度下降 (SGD) 优化算法和交叉熵损失函数，批量大小为 16，周期数为 20，可以获得最佳的分类效果。给定EEG的神经元响应特性和动态特性，我们采用循环神经网络(RNN)在卷积操作后捕获激活序列，长短期记忆网络(LSTM, Hochreiter & Schmidhuber, 1997)是一个记忆改进的循环神经网络，它的“遗忘门”结构允许它有选择地输入下一个单元的信息。同样，我们在测试后在 256 个单元的层中获得了最好的结果。

        通过优化交叉熵损失函数进行训练。CONVNET 中的权重共享导致不同层的梯度非常不同，因此在应用 SGD 时通常使用较小的学习率。Adam 已被证明在用于训练神经网络和多层神经网络时具有具有竞争力的快速收敛速度。我们使用 Adam 算法训练循环卷积网络，第一和第二矩的学习率为 10，衰减率分别为 0.9 和 0.999。批量大小设置为 20。

        为了通过优化交叉熵损失函数来训练我们的网络，卷积层之间的权重共享将导致不同层的不同梯度，我们使用一个小的学习率来避免错误。Adam 已被证明在用于训练各种神经网络时具有较快的收敛速度，因此我们使用 Adam 来训练我们的循环卷积神经网络，学习因子设置为 10，第一和第二矩的衰减率分别为 0.9 和 0.999，批量大小设置为 20。此外，由于使用了 VGG 系统，我们设置了 0.5 的 dropout 率以避免过度拟合，同时尝试随机添加噪声生成的增强数据来训练网络。我们比较了疲劳脑电图分类中常用的各种方法，比较了模糊熵和特征表达式中的样本熵，并比较了分类器中的支持向量机、随机森林、稀疏逻辑回归和深度置信网络(DBN)。

---

代码分析

        在解决问题之前，需要先搞明白文中用于分类的图像是如何生成的。

import numpy as np
import scipy.io as sio
from Utils_Bashivan import *
import matplotlib.pyplot as plt


feats = sio.loadmat('Sample Data/FeatureMat_timeWin.mat')['features']
print(feats.shape)
# (2670, 1345)
locs = sio.loadmat('Sample Data/Neuroscan_locs_orig.mat')
print(locs)
locs_3d = locs['A']
locs_2d = []
# 三维坐标转二维坐标
for e in locs_3d:
    locs_2d.append(azim_proj(e))
"""
    gen_images:
    给定2D空间中的电极位置和每个电极的多个特征值，生成EEG图像
    参数 
    locs: [n_electrodes, 2] 包含每个电极的X，Y坐标
    features: [n_samples, n_features] 特征矩阵
    n_gridpoints: 输出图像中的像素数
    normalize: 是否标准化每个频带的所有样本
    augment: 是否生成增强图像
    pca: 是否基于PCA进行数据增强
    std_mult: 添加噪声的标准差乘数
    n_components: 保留的PCA中的组件数
    edgeless: 如果为True，则通过在图像的四个角添加人工通道来生成无边缘图像
"""
images_timewin = np.array([gen_images(np.array(locs_2d),
                                      feats[:, i * 192:(i + 1) * 192], 32, normalize=True) for i in
                           range(int(feats.shape[1] / 192))
                           ])
# 返回数据为[samples, colors, W, H]color=0为red, color=1为green, color=2为blue，每个sample对应一个图像
print(images_timewin.shape)
# (7, 2670, 3, 32, 32)[]
# 提取图片并plt第一张，包括三种颜色
rgb_image = np.zeros((32, 32, 3), dtype=np.uint8)
rgb_image[:, :, 0] = (images_timewin[0, 0, 0, :, :] * 255).astype(np.uint8)
rgb_image[:, :, 1] = (images_timewin[0, 0, 1, :, :] * 255).astype(np.uint8)
rgb_image[:, :, 2] = (images_timewin[0, 0, 2, :, :] * 255).astype(np.uint8)
plt.imshow(rgb_image)
# 存储图片
plt.savefig("Sample Data/image_exem_1.png")
plt.show()


sio.savemat("Sample Data/images_time.mat",{"img":images_timewin})
print("Images Created and Save in Sample Dat/images_time")

         根据前文中给出的代码库的代码，本人增加了一个test.py文件来测试图像生成的功能，其中用到的几个主要函数如下：

def azim_proj(pos):
    """
    计算输入点在3D笛卡尔坐标中的方位等距投影。
    想象一个平面被放置在（切线）球体上。
    如果球体内的光源将经线投影到平面
    结果将是一个平面或方位角地图投影。

    :param pos: position in 3D Cartesian coordinates
    :return: projected coordinates using Azimuthal Equidistant Projection
    """
    # 球坐标
    [r, elev, az] = cart2sph(pos[0], pos[1], pos[2])
    return pol2cart(az, m.pi / 2 - elev)

def cart2sph(x, y, z):
    """
    将笛卡尔坐标转换为球坐标
    Transform Cartesian coordinates to spherical
    :param x: X coordinate
    :param y: Y coordinate
    :param z: Z coordinate
    :return: radius, elevation, azimuth
    """
    x2_y2 = x**2 + y**2
    # 计算球面半径
    r = m.sqrt(x2_y2 + z**2)                    # r
    # 计算极角和方位角
    # 极角是从x轴到矢量的角度
    # 方位角是从x轴到矢量投影的角度
    elev = m.atan2(z, m.sqrt(x2_y2))            # Elevation
    az = m.atan2(y, x)                          # Azimuth
    return r, elev, az


def pol2cart(theta, rho):
    """
    将极坐标转换为二维平面坐标
    Transform polar coordinates to Cartesian
    :param theta: angle value
    :param rho: radius value
    :return: X, Y
    """
    return rho * m.cos(theta), rho * m.sin(theta)

def gen_images(locs, features, n_gridpoints, normalize=True,
               augment=False, pca=False, std_mult=0.1, n_components=2, edgeless=False):
    """
    给定2D空间中的电极位置和每个电极的多个特征值，生成EEG图像
    参数
    locs: [n_electrodes, 2] 包含每个电极的X，Y坐标
    features: [n_samples, n_features] 特征矩阵
    n_gridpoints: 输出图像中的像素数
    normalize: 是否标准化每个频带的所有样本
    augment: 是否生成增强图像
    pca: 是否基于PCA进行数据增强
    std_mult: 添加噪声的标准差乘数
    n_components: 保留的PCA中的组件数
    edgeless: 如果为True，则通过在图像的四个角添加人工通道来生成无边缘图像
    """
    """
    images_timewin = np.array([gen_images(np.array(locs_2d),
                                      feats[:, i * 192:(i + 1) * 192], 32, normalize=True) for i in
                           range(int(feats.shape[1] / 192))
                           ])
    """
    feat_array_temp = []
    # 电极数量
    nElectrodes = locs.shape[0]     # Number of electrodes
    # 测试特征向量长度是否可被电极数量整除
    # Test whether the feature vector length is divisible by number of electrodes
    assert features.shape[1] % nElectrodes == 0
    # 颜色数量为特征向量长度除以电极数量（每次提取出的一种特征包括全通道的数据，192则说明提取出了三种特征）
    n_colors = int(features.shape[1] / nElectrodes)
    for c in range(n_colors):
        # 每次提取出一种特征
        feat_array_temp.append(features[:, c * nElectrodes : nElectrodes * (c+1)])
    if augment:
        if pca:
            for c in range(n_colors):
                feat_array_temp[c] = augment_EEG(feat_array_temp[c], std_mult, pca=True, n_components=n_components)
        else:
            for c in range(n_colors):
                feat_array_temp[c] = augment_EEG(feat_array_temp[c], std_mult, pca=False, n_components=n_components)
    n_samples = features.shape[0]
    # 插值,mgrid函数用于生成多维结构化的网格点
    """
    在 locs 数组的 x 和 y 坐标范围内创建一个 32x32 的二维网格。
    这个网格包含了 (x, y) 坐标对，其中 x 和 y 分别在范围内均匀分布了 32 个点。
    """
    grid_x, grid_y = np.mgrid[
                     min(locs[:, 0]):max(locs[:, 0]):n_gridpoints*1j,
                     min(locs[:, 1]):max(locs[:, 1]):n_gridpoints*1j
                     ]
    # 生成插值后的图像
    temp_interp = []
    # 生成空图像（颜色数，样本数，网格数，网格数）
    for c in range(n_colors):
        temp_interp.append(np.zeros([n_samples, n_gridpoints, n_gridpoints]))

    # Generate edgeless images
    if edgeless:
        min_x, min_y = np.min(locs, axis=0)
        max_x, max_y = np.max(locs, axis=0)
        locs = np.append(locs, np.array([[min_x, min_y], [min_x, max_y], [max_x, min_y], [max_x, max_y]]), axis=0)
        for c in range(n_colors):
            feat_array_temp[c] = np.append(feat_array_temp[c], np.zeros((n_samples, 4)), axis=1)

    # Interpolating
    # 插值
    # 遍历每个样本，每个样本的每个特征，每个特征的每个电极
    # griddata函数用于插值，实现方法是基于三角剖分的线性插值
    # 这里griddata的第一个参数是locs，第二个参数是feat_array_temp[c][i, :]，第三个参数是grid_x, grid_y
    for i in range(n_samples):
        for c in range(n_colors):
            temp_interp[c][i, :, :] = griddata(locs, feat_array_temp[c][i, :], (grid_x, grid_y),
                                               method='cubic', fill_value=np.nan)
        print('Interpolating {0}/{1}\r'.format(i + 1, n_samples), end='\r')

    # Normalizing
    # 标准化
    for c in range(n_colors):
        if normalize:
            temp_interp[c][~np.isnan(temp_interp[c])] = \
                scale(temp_interp[c][~np.isnan(temp_interp[c])])
        temp_interp[c] = np.nan_to_num(temp_interp[c])
    return np.swapaxes(np.asarray(temp_interp), 0, 1)     # swap axes to have [samples, colors, W, H]

        最后运行test.py文件得到图像如下：

 图3 测试生成图像

        测试成功。