CN102973279A - 独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法 - Google Patents

Abstract

独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法，属于血红蛋白浓度检测技术领域。它解决了近红外脑机接口检测中由于人体生理干扰造成检测获得的氧合血红蛋白浓度变化和还原血红蛋白浓度变化量不准确，而影响脑功能活动信号准确提取的问题。它通过检测器记录大脑安静状态下和处于诱发激励时漫反射光强，获得光密度变化量的时间序列

和

和

再获取Δ[HbO₂]^N(k)、Δ[HHb]^N(k)、Δ[HbO₂]^F(k)Δ[HHb]^F(k)；用x₁(k)表示步骤二中的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)；用x₂(k)表示步骤二中Δ[HbO₂]^F(k)或Δ[HHb]^F(k)；推算出脑功能信号表达式s(k)；求解脑功能信号s(k)。本发明适用于脑机接口的信号检测。

Description

独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法

技术领域

本发明涉及独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法，属于血红蛋白浓度检测技术领域。

背景技术

脑机接口是一种由人脑与计算机或其它电子设备建立起来的，基于大脑功能的电生理测量，而不依赖于外周神经和肌肉组织这些常规的人脑信息输出通道，实现人与外界信息交流和控制的全新通讯系统。通过分析脑信号将用户的运动等意图转换为语言、设备的控制输入量等，使用户直接通过脑信号与外面的环境进行实时的交互，从而绕开了人类神经末梢和肌肉等通常的信息通道，建立能直接“让思想变成行动”的对外信息交流和控制新途径。

目前，脑机接口研究的基本方法是提取并识别特定意识活动的脑功能状态特征信息，现有的主要技术包括脑电图、脑磁图、正电子放射层扫描术、功能磁共振。基于神经元的活动、局部能量代谢与局部血液动力学之间存在着一定的关系，通过测量脑组织对近红外光波的吸收特性，能够提供基于氧合血红蛋白和还原血红蛋白浓度等信息的血液动力学变化。因此，利用近红外光谱技术测量该区域光学参数、血氧及血液动力学参数信息，可以获取大脑皮质在肢体运动、视觉、听觉、触觉以及语言等刺激激励时的功能响应，用于脑机接口的研究。功能近红外光谱技术可以安全、便携、经济以及非侵入式的检测脑活动等特性，基于较高的时间分辨率和合理的空间分辨率，在脑机接口研究中具有一定的发展潜力。

然而，通过近红外光谱技术进行诱发激励时脑功能活动的检测，会受到人体的生理活动如心脏跳动、呼吸、低频振荡、超低频振荡的影响，称之为生理干扰。这种生理干扰不但出现在头皮、颅骨和脑脊液等外层脑组织中，也出现在脑灰质和脑白质等深层脑组织中，这些原因会使得近红外脑机接口检测获得的氧合血红蛋白浓度变化和还原血红蛋白浓度变化量不准确，进而严重影响脑功能活动信号的准确提取。

发明内容

本发明是为了解决近红外脑机接口检测中由于人体生理干扰造成检测获得的氧合血红蛋白浓度变化和还原血红蛋白浓度变化量不准确，而影响脑功能活动信号准确提取的问题，提供了一种独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法。

本发明所述独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法，它包括以下步骤：

步骤一：采用近红外探头靠近待检测头部的头皮表面，使得该近红外探头发射的近红外光入射至待测脑组织，该近红外探头由双波长光源S、检测器D1和检测器D2构成，其中双波长光源S到检测器D1之间的直线距离为r₁，5mm≤r₁≤15mm，双波长光源S到检测器D2之间的直线距离为r₂，30mm≤r₂≤45mm；检测器D1用于感应外层脑组织的血液动力学变化，检测器D2用于感应大脑皮质的血液动力学变化；

步骤二：通过检测器D1记录大脑安静状态下的漫反射光强和大脑处于诱发激励状态下的漫反射光强，获得大脑安静状态下双波长光源S波长为λ₁时、对应的漫反射光的光密度变化量的时间序列

及获得大脑处于诱发激励状态下双波长光源S波长为λ₂时、对应的漫反射光的光密度变化量的时间序列

通过检测器D2记录大脑安静状态下的漫反射光强和大脑处于诱发激励状态下的漫反射光强，获得大脑安静状态下双波长光源S波长为λ₁时、对应的漫反射光的光密度变化量的时间序列

及获得大脑处于诱发激励状态下双波长光源S波长为λ₂时、对应的漫反射光的光密度变化量的时间序列

k为获得时间序列的点数，k＝1，2，...，M，M为正整数；

步骤三：根据步骤二中获得的时间序列

和时间序列

采用修正朗伯比尔定律获取与检测器D1的检测信号对应的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)，及获取与检测器D1的检测信号对应的还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)：

$\mathrm{\Δ}{\left[{\mathrm{HbO}}_2\right]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})}{r_1({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))},$

$\mathrm{\Δ}{\left[\mathrm{HHb}\right]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})}{r_1({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))};$

其中，ε_HHb(λ₁)为双波长光源S的波长为λ₁时还原血红蛋白的消光系数，

ε_HHb(λ₂)为双波长光源S的波长为λ₂时还原血红蛋白的消光系数，

为双波长光源S的波长为λ₁时氧合血红蛋白的消光系数，

为双波长光源S的波长为λ₂时氧合血红蛋白的消光系数，

DPF为差分路径因子；

根据步骤二中获得的时间序列

和时间序列

采用修正朗伯比尔定律获取与检测器D2的检测信号对应的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^F(k)，及获取与检测器D2的检测信号对应的还原血红蛋白浓度变化量的时间序列Δ[HHb]^F(k)：

$\mathrm{\Δ}{\left[{\mathrm{HbO}}_2\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)/\mathrm{DPF})}{r_2({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))},$

$\mathrm{\Δ}{\left[\mathrm{HHb}\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)/\mathrm{DPF})}{r_2({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))};$

步骤四：用x₁(k)表示步骤三中的Δ[HbO₂]^F(k)或Δ[HHb]^F(k)，将x₁(k)作为脑机接口检测中的测量通道信号，x₁(k)扩展至加噪模型为：

x₁(k)＝s(k)+n(k)，

其中s(k)为脑机接口信号中的血液动力学信号，n(k)为生理干扰信号；

用x₂(k)表示步骤三中的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)，作为脑机接口检测中的虚拟通道信号，则有观测矩阵x：

$x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}s+n\\ u\end{array}\right]=\left[\begin{array}{cc}1& a\\ 0& 1\end{array}\right]\left[\begin{array}{c}s\\ u\end{array}\right]=\mathrm{Ar},$ $A=\left[\begin{array}{cc}1& a\\ 0& 1\end{array}\right],$ $r=\left[\begin{array}{c}s\\ u\end{array}\right],$

其中，x₁为整个观测的时间序列x₁(k)，是一维相量；

x₂为整个观测的时间序列x₂(k)，是一维相量，即检测器D2检测的生理干扰混合相量u；

s为血液动力学信号的整个时间序列；

n为整个观测的时间序列x₁(k)中混叠的噪声，

a为检测器D2检测获得的干扰信号与检测器D1检测获得信号的比例权值；

A为混合矩阵，

r为分离的脑机接口信号和噪声组成的相量；

步骤五、用独立成分分析ICA算法，在混合矩阵A和血液动力学信号的整个时间序列s未知的情况下，根据观测矩阵x确定分离矩阵W，使得变换后的输出Y＝Wx；

步骤六、根据步骤五中的Y＝Wx推算出二维矩阵，两行分别表示为l和d，l和d分别与s和u满足线性相关条件，用bl(k)+c作为s(k)的最优估计，则残差e(k)为：

e(k)＝x₁(k)-bl(k)-c；

步骤七、利用最小二乘算法，采用残差e(k)的累计平方误差性能函数J作为代价函数：

$J=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{e}^2\left(k\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(x_1\left(k\right)-\mathrm{bl}\left(k\right)-c)}^2,$

求取使J最小的系数b和c，再将求取的系数b和c带入步骤六中的bl(k)+c公式中，即可获得脑机接口信号中的血液动力学信号s(k)。

步骤二中光密度变化量的时间序列和光密度变化量的时间序列

按如下公式获取：

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^N\left(k\right)=\log I_{\mathrm{base}}^N\left({\mathrm{\λ}}_1\right)/I_{\mathrm{stim}}^N\left({\mathrm{\λ}}_1\right),$

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)=\log I_{\mathrm{base}}^F\left({\mathrm{\λ}}_1\right)/I_{\mathrm{stim}}^F\left({\mathrm{\λ}}_1\right),$

其中：

为双波长光源S波长为λ₁时，大脑安静状态下检测器D1测得的出射光强；

为双波长光源S波长为λ₁时，大脑安静状态下检测器D2测得的出射光强；

为双波长光源S波长为λ₁时，大脑处于诱发激励状态下检测器D1测得的出射光强；

为双波长光源S波长为λ₁时，大脑处于诱发激励状态下检测器D2测得的出射光强，

步骤二中光密度变化量的时间序列

和光密度变化量的时间序列

按如下公式获取：

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)=\log I_{\mathrm{base}}^N\left({\mathrm{\λ}}_2\right)/I_{\mathrm{stim}}^N\left({\mathrm{\λ}}_2\right),$

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)=\log I_{\mathrm{base}}^F\left({\mathrm{\λ}}_2\right)/I_{\mathrm{stim}}^F\left({\mathrm{\λ}}_2\right),$

其中：

为双波长光源S波长为λ₂时，大脑安静状态下检测器D1测得的出射光强，

为双波长光源S波长为λ₂时，大脑处于诱发激励状态下检测器D1测得的出射光强；

为双波长光源S波长为λ₂时，大脑安静状态下检测器D2测得的出射光强，

为双波长光源S波长为λ₂时，大脑处于诱发激励状态下检测器D2测得的出射光强。

步骤七中残差e(k)的获得方法为：

首先，通过最小二乘估计准则表示使残差e(k)的累计平方误差性能函数J最小，J表示为：

$J=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{e}^2\left(k\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(x_1\left(k\right)-\mathrm{bl}\left(k\right)-c)}^2;$

其次，求解最优系数b和c：

对J相对于b，c求导，并将求导结果置为0，即：

$\∂J/\∂b=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}2[x_1\left(k\right)-\mathrm{bl}\left(k\right)-c](-l\left(k\right))=2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}(c+\mathrm{bl}\left(k\right)-x_1\left(k\right))l\left(k\right)=0,$

$\∂J/\∂c=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}2[x_1\left(k\right)-\mathrm{bl}\left(k\right)-c](-1)=2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\left[c+\mathrm{bl}\left(k\right)-x_1\left(k\right)\right]=0,$

则：

$b\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{l}^2\left(k\right)-\underset{k=1}{\overset{N}{\mathrm{\Σ}}}l\left(k\right)x_1\left(k\right)+c\underset{k=1}{\overset{N}{\mathrm{\Σ}}}l\left(k\right)=0,$

$b\underset{k=1}{\overset{N}{\mathrm{\Σ}}}l\left(k\right)-\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{x}_1\left(k\right)+\mathrm{Nc}=0;$

最后，求解脑机接口信号中的血液动力学信号s(k)：

s(k)＝bl(k)+c。

本发明的优点：本发明是针对近红外脑机接口研究过程中大脑特定区域氧合血红蛋白浓度变化Δ[HbO₂]和还原血红蛋白浓度变化Δ[HHb]易受到心跳、呼吸及人体自发低频振荡噪声干扰的状况，在多距测量方法的基础上，考虑近端检测器D1获得的血液动力学参数和远端检测器D2受到的生理干扰具有相关性，通过独立成分分析算法对远端和近端检测器测量结果进行分解，并通过对分解的独立成分分量建立线性映射来估计测量信号中的脑功能信号。独立成分分析能将复合信号分解为一系列的独立成分分量，并且分解的独立成分分量具有很好的独立性，适用于盲源信号的分析。本发明通过用独立成分分析算法分解远端检测器和近端检测器测得的外层组织血液动力学参数，从而获得表示大脑皮质血液动力学参数的估计，实现对脑功能信号的准确提取。

附图说明

图1为近红外探头的检测状态示意图，其中a表示头皮，b表示颅骨，c表示脑脊液，d表示脑灰质，e表示脑白质。

具体实施方式

具体实施方式一：下面结合图1说明本实施方式，本实施方式所述独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法，它包括以下步骤：

步骤一：采用近红外探头靠近待检测头部的头皮表面，使得该近红外探头发射的近红外光入射至待测脑组织，该近红外探头由双波长光源S、检测器D1和检测器D2构成，其中双波长光源S到检测器D1之间的直线距离为r₁，5mm≤r₁≤15mm，双波长光源S到检测器D2之间的直线距离为r₂，30mm≤r₂≤45mm；检测器D1用于感应外层脑组织的血液动力学变化，检测器D2用于感应大脑皮质的血液动力学变化；

步骤二：通过检测器D1记录大脑安静状态下的漫反射光强和大脑处于诱发激励状态下的漫反射光强，获得大脑安静状态下双波长光源S波长为λ₁时、对应的漫反射光的光密度变化量的时间序列

及获得大脑处于诱发激励状态下双波长光源S波长为λ₂时、对应的漫反射光的光密度变化量的时间序列

通过检测器D2记录大脑安静状态下的漫反射光强和大脑处于诱发激励状态下的漫反射光强，获得大脑安静状态下双波长光源S波长为λ₁时、对应的漫反射光的光密度变化量的时间序列

及获得大脑处于诱发激励状态下双波长光源S波长为λ₂时、对应的漫反射光的光密度变化量的时间序列

k为获得时间序列的点数，k＝1，2，...，M，M为正整数；

步骤三：根据步骤二中获得的时间序列

和时间序列采用修正朗伯比尔定律获取与检测器D1的检测信号对应的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)，及获取与检测器D1的检测信号对应的还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)：

$\mathrm{\Δ}{\left[{\mathrm{HbO}}_2\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)/\mathrm{DPF})}{r_1({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))},$

$\mathrm{\Δ}{\left[\mathrm{HHb}\right]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})}{r_1({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))};$

其中，ε_HHb(λ₁)为双波长光源S的波长为λ₁时还原血红蛋白的消光系数，

ε_HHb(λ₂)为双波长光源S的波长为λ₂时还原血红蛋白的消光系数，

为双波长光源S的波长为λ₁时氧合血红蛋白的消光系数，

为双波长光源S的波长为λ₂时氧合血红蛋白的消光系数，

DPF为差分路径因子；

根据步骤二中获得的时间序列

和时间序列

采用修正朗伯比尔定律获取与检测器D2的检测信号对应的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^F(k)，及获取与检测器D2的检测信号对应的还原血红蛋白浓度变化量的时间序列Δ[HHb]^F(k)：

$\mathrm{\Δ}{\left[{\mathrm{HbO}}_2\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)/\mathrm{DPF})}{r_2({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))},$

$\mathrm{\Δ}{\left[\mathrm{HHb}\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)/\mathrm{DPF})}{r_2({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))};$

步骤四：用x₁(k)表示步骤三中的Δ[HbO₂]^F(k)或Δ[HHb]^F(k)，将x₁(k)作为脑机接口检测中的测量通道信号，x₁(k)扩展至加噪模型为：

x₁(k)＝s(k)+n(k)，

其中s(k)为脑机接口信号中的血液动力学信号，n(k)为生理干扰信号；

用x₂(k)表示步骤三中的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)，作为脑机接口检测中的虚拟通道信号，则有观测矩阵x：

$x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}s+n\\ u\end{array}\right]=\left[\begin{array}{cc}1& a\\ 0& 1\end{array}\right]\left[\begin{array}{c}s\\ u\end{array}\right]=\mathrm{Ar},$ $A=\left[\begin{array}{cc}1& a\\ 0& 1\end{array}\right],$ $r=\left[\begin{array}{c}s\\ u\end{array}\right],$

其中，x₁为整个观测的时间序列x₁(k)，是一维相量；

x₂为整个观测的时间序列x₂(k)，是一维相量，即检测器D2检测的生理干扰混合相量u；

s为血液动力学信号的整个时间序列；

n为整个观测的时间序列x₁(k)中混叠的噪声，

a为检测器D2检测获得的干扰信号与检测器D1检测获得信号的比例权值；

A为混合矩阵，

r为分离的脑机接口信号和噪声组成的相量；

步骤五、用独立成分分析ICA算法，在混合矩阵A和血液动力学信号的整个时间序列s未知的情况下，根据观测矩阵x确定分离矩阵W，使得变换后的输出Y＝Wx；

步骤六、根据步骤五中的Y＝Wx推算出二维矩阵，两行分别表示为l和d，l和d分别与s和u满足线性相关条件，用bl(k)+c作为s(k)的最优估计，则残差e(k)为：

e(k)＝x₁(k)-bl(k)-c；

步骤七、利用最小二乘算法，采用残差e(k)的累计平方误差性能函数J作为代价函数：

$J=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{e}^2\left(k\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(x_1\left(k\right)-\mathrm{bl}\left(k\right)-c)}^2,$

求取使J最小的系数b和c，再将求取的系数b和c带入步骤六中的bl(k)+c公式中，即可获得脑机接口信号中的血液动力学信号s(k)。

本实施方式中，x₁(k)包含s(k)以及n(k)，x₂(k)为近距离光源检测器获得的血液动力学信息，主要为生理干扰信号且与x1(k)中噪声信号n(k)相关，从而获得观测矩阵x。

具体实施方式二：本实施方式对实施方式一的进一步说明，本实施方式所述双波长光源S发出的两种波长分别为λ₁＝760nm，λ₂＝850nm。

具体实施方式三：本实施方式对实施方式一或二的进一步说明，本实施方式所述双波长光源S到检测器D1之间的直线距离r₁为10mm，双波长光源S到检测器D2之间的直线距离r₂为40mm。

本实施方式中设置的两个检测器间距约为近红外光探测深度的两倍，这样设置能够使检测器D2检测的近红外光可有效穿入大脑皮层，检测器D1检测的近红外光仅穿入头外层脑组织。然后将获得的光密度变化通过修正朗伯比尔定律转变为氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)、Δ[HbO₂]^F(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)、Δ[HHb]^F(k)。将时间序列Δ[HbO₂]^N(k)、Δ[HbO₂]^F(k)或时间序列Δ[HHb]^N(k)、Δ[HHb]^F(k)构成二维观测矩阵x。通过独立成分分析算法根据观测矩阵x确定分离矩阵W，使得变换后的输出Y＝Wx。矩阵Y的两行分别表示为l和d，l和d分别与s和u满足线性相关条件。将得到的分量l进行线性组合估计脑功能信号s，构建鳌头e(k)。通过最小二乘估计准则求解使残差e(k)的累计平方误差性能函数J最小，获得线性组合参数b和c，从而解得剔除生理干扰的脑机接口研究中的脑功能信号。

具体实施方式四：本实施方式对实施方式一、二或三的进一步说明，本实施方式所述步骤二中光密度变化量的时间序列

和光密度变化量的时间序列按如下公式获取：

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^N\left(k\right)=\log I_{\mathrm{base}}^N\left({\mathrm{\λ}}_1\right)/I_{\mathrm{stim}}^N\left({\mathrm{\λ}}_1\right),$

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)=\log I_{\mathrm{base}}^F\left({\mathrm{\λ}}_1\right)/I_{\mathrm{stim}}^F\left({\mathrm{\λ}}_1\right),$

其中：为双波长光源S波长为λ₁时，大脑安静状态下检测器D1测得的出射光强；

为双波长光源S波长为λ₁时，大脑安静状态下检测器D2测得的出射光强；

为双波长光源S波长为λ₁时，大脑处于诱发激励状态下检测器D1测得的出射光强；

为双波长光源S波长为λ₁时，大脑处于诱发激励状态下检测器D2测得的出射光强，

步骤二中光密度变化量的时间序列和光密度变化量的时间序列

按如下公式获取：

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)=\log I_{\mathrm{base}}^N\left({\mathrm{\λ}}_2\right)/I_{\mathrm{stim}}^N\left({\mathrm{\λ}}_2\right),$

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)=\log I_{\mathrm{base}}^F\left({\mathrm{\λ}}_2\right)/I_{\mathrm{stim}}^F\left({\mathrm{\λ}}_2\right),$

其中：

为双波长光源S波长为λ₂时，大脑安静状态下检测器D1测得的出射光强，

为双波长光源S波长为λ₂时，大脑处于诱发激励状态下检测器D1测得的出射光强；

为双波长光源S波长为λ₂时，大脑安静状态下检测器D2测得的出射光强，

为双波长光源S波长为λ₂时，大脑处于诱发激励状态下检测器D2测得的出射光强。

在近红外波段氧合血红蛋白HbO₂和还原血红蛋白HHb是吸收体，并且其吸收谱存在显著差异。因此，基于连续光谱技术的近红外脑功能检测是HbO₂和HHb的浓度变化。

具体实施方式五：本实施方式对实施方式一、二、三或四的进一步说明，本实施方式所述步骤七中残差e(k)的获得方法为：

首先，通过最小二乘估计准则表示使残差e(k)的累计平方误差性能函数J最小，J表示为：

$J=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{e}^2\left(k\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(x_1\left(k\right)-\mathrm{bl}\left(k\right)-c)}^2;$

其次，求解最优系数b和c：

对J相对于b，c求导，并将求导结果置为0，即：

$\∂J/\∂b=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}2[x_1\left(k\right)-\mathrm{bl}\left(k\right)-c](-l\left(k\right))=2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}(c+\mathrm{bl}\left(k\right)-x_1\left(k\right))l\left(k\right)=0,$

$\∂J/\∂c=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}2[x_1\left(k\right)-\mathrm{bl}\left(k\right)-c](-1)=2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\left[c+\mathrm{bl}\left(k\right)-x_1\left(k\right)\right]=0,$

则：

$b\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{l}^2\left(k\right)-\underset{k=1}{\overset{N}{\mathrm{\Σ}}}l\left(k\right)x_1\left(k\right)+c\underset{k=1}{\overset{N}{\mathrm{\Σ}}}l\left(k\right)=0,$

$b\underset{k=1}{\overset{N}{\mathrm{\Σ}}}l\left(k\right)-\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{x}_1\left(k\right)+\mathrm{Nc}=0;$

最后，求解脑机接口信号中的血液动力学信号s(k)：

s(k)＝bl(k)+c。

Claims (5)

1.一种独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法，其特征在于，它包括以下步骤：

步骤一：采用近红外探头靠近待检测头部的头皮表面，使得该近红外探头发射的近红外光入射至待测脑组织，该近红外探头由双波长光源S、检测器D1和检测器D2构成，其中双波长光源S到检测器D1之间的直线距离为r₁，5mm≤r₁≤15mm，双波长光源S到检测器D2之间的直线距离为r₂，30mm≤r₂≤45mm；检测器D1用于感应外层脑组织的血液动力学变化，检测器D2用于感应大脑皮质的血液动力学变化；

步骤二：通过检测器D1记录大脑安静状态下的漫反射光强和大脑处于诱发激励状态下的漫反射光强，获得大脑安静状态下双波长光源S波长为λ₁时、对应的漫反射光的光密度变化量的时间序列

及获得大脑处于诱发激励状态下双波长光源S波长为λ₂时、对应的漫反射光的光密度变化量的时间序列

通过检测器D2记录大脑安静状态下的漫反射光强和大脑处于诱发激励状态下的漫反射光强，获得大脑安静状态下双波长光源S波长为λ₁时、对应的漫反射光的光密度变化量的时间序列

及获得大脑处于诱发激励状态下双波长光源S波长为λ₂时、对应的漫反射光的光密度变化量的时间序列

k为获得时间序列的点数，k＝1，2，...，M，M为正整数；

步骤三：根据步骤二中获得的时间序列

和时间序列

采用修正朗伯比尔定律获取与检测器D1的检测信号对应的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)，及获取与检测器D1的检测信号对应的还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)：

$\mathrm{\Δ}{\left[{\mathrm{HbO}}_2\right]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})}{r_1({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))},$

$\mathrm{\Δ}{\left[\mathrm{HHb}\right]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})}{r_1({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))};$

其中，ε_HHb(λ₁)为双波长光源S的波长为λ₁时还原血红蛋白的消光系数，

ε_HHb(λ₂)为双波长光源S的波长为λ₂时还原血红蛋白的消光系数，

为双波长光源S的波长为λ₁时氧合血红蛋白的消光系数，

为双波长光源S的波长为λ₂时氧合血红蛋白的消光系数，

DPF为差分路径因子；

根据步骤二中获得的时间序列

和时间序列

采用修正朗伯比尔定律获取与检测器D2的检测信号对应的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^F(k)，及获取与检测器D2的检测信号对应的还原血红蛋白浓度变化量的时间序列Δ[HHb]^F(k)：

$\mathrm{\Δ}{\left[{\mathrm{HbO}}_2\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)/\mathrm{DPF})}{r_2({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))},$

$\mathrm{\Δ}{\left[\mathrm{HHb}\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)/\mathrm{DPF})}{r_2({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right))};$

步骤四：用x₁(k)表示步骤三中的Δ[HbO₂]^F(k)或Δ[HHb]^F(k)，将x₁(k)作为脑机接口检测中的测量通道信号，x₁(k)扩展至加噪模型为：

x₁(k)＝s(k)+n(k)，

其中s(k)为脑机接口信号中的血液动力学信号，n(k)为生理干扰信号；

用x₂(k)表示步骤三中的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)，作为脑机接口检测中的虚拟通道信号，则有观测矩阵x：

$x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}s+n\\ u\end{array}\right]=\left[\begin{array}{cc}1& a\\ 0& 1\end{array}\right]\left[\begin{array}{c}s\\ u\end{array}\right]=\mathrm{Ar},$ $A=\left[\begin{array}{cc}1& a\\ 0& 1\end{array}\right],$ $r=\left[\begin{array}{c}s\\ u\end{array}\right],$

其中，x₁为整个观测的时间序列x₁(k)，是一维相量；

x₂为整个观测的时间序列x₂(k)，是一维相量，即检测器D2检测的生理干扰混合相量u；

s为血液动力学信号的整个时间序列；

n为整个观测的时间序列x₁(k)中混叠的噪声，

a为检测器D2检测获得的干扰信号与检测器D1检测获得信号的比例权值；

A为混合矩阵，

r为分离的脑机接口信号和噪声组成的相量；

步骤五、用独立成分分析ICA算法，在混合矩阵A和血液动力学信号的整个时间序列s未知的情况下，根据观测矩阵x确定分离矩阵W，使得变换后的输出Y＝Wx；

步骤六、根据步骤五中的Y＝Wx推算出二维矩阵，两行分别表示为l和d，l和d分别与s和u满足线性相关条件，用bl(k)+c作为s(k)的最优估计，则残差e(k)为：

e(k)＝x₁(k)-bl(k)-c；

步骤七、利用最小二乘算法，采用残差e(k)的累计平方误差性能函数J作为代价函数：

$J=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{e}^2\left(k\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(x_1\left(k\right)-\mathrm{bl}\left(k\right)-c)}^2,$

求取使J最小的系数b和c，再将求取的系数b和c带入步骤六中的bl(k)+c公式中，即可获得脑机接口信号中的血液动力学信号s(k)。

2.根据权利要求1所述的独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法，其特征在于，双波长光源S发出的两种波长分别为λ₁＝760nm，λ₂＝850nm。

3.根据权利要求1所述的独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法，其特征在于，双波长光源S到检测器D1之间的直线距离r₁为10mm，双波长光源S到检测器D2之间的直线距离r₂为40mm。

4.根据权利要求1所述的独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法，其特征在于，

步骤二中光密度变化量的时间序列

和光密度变化量的时间序列

按如下公式获取：

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^N\left(k\right)=\log I_{\mathrm{base}}^N\left({\mathrm{\λ}}_1\right)/I_{\mathrm{stim}}^N\left({\mathrm{\λ}}_1\right),$

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_1}^F\left(k\right)=\log I_{\mathrm{base}}^F\left({\mathrm{\λ}}_1\right)/I_{\mathrm{stim}}^F\left({\mathrm{\λ}}_1\right),$

其中：

为双波长光源S波长为λ₁时，大脑安静状态下检测器D1测得的出射光强；

为双波长光源S波长为λ₁时，大脑安静状态下检测器D2测得的出射光强；

为双波长光源S波长为λ₁时，大脑处于诱发激励状态下检测器D1测得的出射光强；

为双波长光源S波长为λ₁时，大脑处于诱发激励状态下检测器D2测得的出射光强，

步骤二中光密度变化量的时间序列

和光密度变化量的时间序列

按如下公式获取：

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)=\log I_{\mathrm{base}}^N\left({\mathrm{\λ}}_2\right)/I_{\mathrm{stim}}^N\left({\mathrm{\λ}}_2\right),$

${\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^F\left(k\right)=\log I_{\mathrm{base}}^F\left({\mathrm{\λ}}_2\right)/I_{\mathrm{stim}}^F\left({\mathrm{\λ}}_2\right),$

其中：为双波长光源S波长为λ₂时，大脑安静状态下检测器D1测得的出射光强，

为双波长光源S波长为λ₂时，大脑处于诱发激励状态下检测器D1测得的出射光强；

为双波长光源S波长为λ₂时，大脑安静状态下检测器D2测得的出射光强，

为双波长光源S波长为λ₂时，大脑处于诱发激励状态下检测器D2测得的出射光强。

5.根据权利要求4所述的独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法，其特征在于，步骤七中残差e(k)的获得方法为：

首先，通过最小二乘估计准则表示使残差e(k)的累计平方误差性能函数J最小，J表示为：

$J=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{e}^2\left(k\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{(x_1\left(k\right)-\mathrm{bl}\left(k\right)-c)}^2;$

其次，求解最优系数b和c：

对J相对于b，c求导，并将求导结果置为0，即：

$\∂J/\∂b=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}2[x_1\left(k\right)-\mathrm{bl}\left(k\right)-c](-l\left(k\right))=2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}(c+\mathrm{bl}\left(k\right)-x_1\left(k\right))l\left(k\right)=0,$

$\∂J/\∂c=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}2[x_1\left(k\right)-\mathrm{bl}\left(k\right)-c](-1)=2\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\left[c+\mathrm{bl}\left(k\right)-x_1\left(k\right)\right]=0,$

则：

$b\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{l}^2\left(k\right)-\underset{k=1}{\overset{N}{\mathrm{\Σ}}}l\left(k\right)x_1\left(k\right)+c\underset{k=1}{\overset{N}{\mathrm{\Σ}}}l\left(k\right)=0,$

$b\underset{k=1}{\overset{N}{\mathrm{\Σ}}}l\left(k\right)-\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{x}_1\left(k\right)+\mathrm{Nc}=0;$

最后，求解脑机接口信号中的血液动力学信号s(k)：

s(k)＝bl(k)+c。

Images