CN102525422B - 基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法 - Google Patents

Abstract

基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法，涉及脑功能信号提取方法。它解决了当脑组织非均匀性严重时现有技术检测脑功能活动过程中氧合血红蛋白浓度变化Δ[HbO₂]和还原血红蛋白浓度变化Δ[HHb]难以检测的问题。本发明通过检测器记录大脑安静状态下和诱发激励时漫反射光强，以获得两个不同波长λ₁和λ₂时的光密度变化量的时间序列：和

和

采用修正朗伯比尔定律获取r₁测得的Δ[HbO₂]^N(k)和Δ[HHb]^N(k)，r₂测得的Δ[HbO₂]^F(k)Δ[HHb]^F(k)；根据获得的所有参数推算出脑功能信号表达式；求解脑功能信号e(k)。本发明适用于医疗领域。

Description

基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法

技术领域

本发明涉及一种信号提取方法，具体涉及基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法。

背景技术

近红外光谱技术能提供脑功能活动过程中的大脑皮层血氧代谢信息——氧合血红蛋白浓度变化Δ[HbO₂]和还原血红蛋白浓度变化Δ[HHb]，可用于脑功能活动的检测。然而，通过近红外光谱技术进行诱发激励时脑功能活动的检测，会受到人体的生理活动如心脏跳动、呼吸、低频振荡、超低频振荡的影响，称之为生理干扰。这种生理干扰不但出现在头皮、颅骨和脑脊液等外层脑组织中，也出现在脑灰质和脑白质等深层脑组织中，严重影响了脑功能活动信号的准确提取。

由于生理干扰来源于人体不同的生理活动，因而具有多个成分。当脑组织非均匀性严重时，将造成不同生理活动在空间不同位置上对生理干扰的“贡献”不同。对于这种情况，比较可行的办法是对不同类型的干扰进行单独估计。一种方法是通过血压检测仪，呼吸计等仪器获得每个生理干扰的参考信号，然后通过卡尔曼滤波跟踪不同的生理干扰，但这种方法需要借助额外的设备；另一种方法是通过多个先验频率的正弦或余弦信号作为生理干扰的参考信号，通过卡尔曼滤波进行生理干扰的估计，但这需要知道被测者生理干扰频率信息的先验知识，但由于个体差异这在实际应用中往往并不易于实现。

发明内容

本发明的目的是为了解决当脑组织非均匀性严重时采用近红外光谱技术检测脑功能活动过程中氧合血红蛋白浓度变化Δ[HbO₂]和还原血红蛋白浓度变化Δ[HHb]难以检测的问题。

本发明方法包括以下步骤：

步骤一、在待测脑组织的头皮表面放置由双波长光源S和检测器D1和D2构成的近红外探头，其中，双波长光源S到检测器D1之间的直线距离为r₁，5mm≤r₁≤10mm，用于敏感外层脑组织的血液动力学变化；双波长光源S到检测器D2之间的直线距离为r₂，30mm≤r₂≤45mm，能够敏感大脑皮质的血液动力学变化，通过检测器记录大脑安静状态下的漫反射光强和大脑处于诱发激励时的漫反射光强，以获得两个不同波长λ₁和λ₂时的光密度变化量的时间序列：和

，

和

，k为时间，k=1,2,...,N；N为正整数，表示在双波长光源S到检测器D1之间的直线距离为r₁且波长为λ₁时光密度变化量的时间序列，

表示在双波长光源S到检测器D1之间的直线距离为r₁且波长为λ₂时光密度变化量的时间序列，

表示在双波长光源S到检测器D2之间的直线距离为r₂且波长为λ₁时光密度变化量的时间序列，表示在双波长光源S到检测器D2之间的直线距离为r₂且波长为λ₂时光密度变化量的时间序列；

步骤二、根据步骤一获得的光密度变化量的时间序列并采用修正朗伯比尔定律获取r₁测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)，r₂测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^F(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^F(k)；

$\mathrm{\Δ}{\left[\mathrm{Hb}{O}_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))},$

$\mathrm{\Δ}{\left[\mathrm{Hb}{O}_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{Hb}{O}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{H{\mathrm{bO}}_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))},$

其中，ε_HHb(λ₁)为探头光源的波长为λ₁时的消光系数，

为探头光源的波长为λ₂时的消光系数，

DPF为差分路径因子；

步骤三、用x(k)表示步骤二中的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)，用经验模态分解将x(k)分解为N个固态模式函数分量IMF分量，将剩余分量作为最后的IMF分量，则x(k)表示为

$\left(k\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{c}_i\left(k\right)$

其中，c_i(k)为分解的IMF分量；

步骤四、用d(k)表示步骤二中的Δ[HbO₂]^F(k)或Δ[HHb]^F(k)，d(k)中包含脑功能活动信号r(k)和生理干扰i(k)，即d(k)=r(k)+i(k)，采用线性映射关系，用c_i(k)的线性组合表示i(k)的估计，即

$\hat{i}\left(k\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(k\right)$

其中，

为i(k)的估计，i=1,2,...,N，w_i(k)为第i个IMF分量的权系数；

步骤五、根据步骤二中的d(k)=r(k)+i(k)和

即可推算出脑功能活动信号估计的表达式：

$e\left(k\right)=d\left(k\right)-\hat{i}\left(k\right)=r\left(k\right)+[i\left(k\right)-\hat{i}\left(k\right)]$

其中，e(k)为脑功能活动信号估计；

步骤六、利用加权最小二乘算法作为代价函数，求取优化系数w_i(k)，再将求取优化的系数w_i(k)带入步骤五中的

公式中，即可获得脑功能活动信号估计e(k)，加权最小二乘算法为：

$\left(k\right)\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{e}^2\left(n\right)$

进一步表示为

$J\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{e}^2\left(n\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{[d\left(n\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(n\right)]}^2$

其中，χ为指数加权因子，χ=0.99;n=1,…k,k为正整数，i=1,2,...,N，N为正整数，求解使J(k)最小的w_i(k)，获得脑功能活动信号估计e(k)。

本发明的优点：本发明方法在多距测量方法的基础上，考虑近端检测器D1获得的血液动力学参数和远端检测器D2受到的生理干扰具有相关性以及每一类型的生理干扰对测量结果的影响可能不同的特点，通过经验模态分解对近端测量结果进行分解得到IMF分量，并通过对IMF分量建立线性映射来估计测量信号中的生理干扰。经验模态分解能将复合信号分解为一系列的固态模式函数，并且分解的固态模式函数具有很好的瞬时频率特性，适用于非线性非平稳信号的分析。本发明通过用经验模态分解算法分解近端检测器测得的外层组织血液动力学参数，从而获得表示外层组织血液动力学参数的固态模式函数分量，并通过优化算法调节不同固态模式函数分量来估计期望信号中的生理干扰，实现对脑功能信号的准确提取的目的。

附图说明

图1是基于多距测量方法的近红外脑功能活动检测探头结构，其中a表示头皮，b表示颅骨，c表示脑脊液，d表示脑灰质，e表示脑白质；图2为基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法原理框图，其中f表示经验模态分解，g表示递归最小二乘算法。

具体实施方式

具体实施方式一：下面结合图1说明本实施方式，本实施方式方法包括以下步骤：

步骤一、在待测脑组织的头皮表面放置由双波长光源S和检测器D1和D2构成的近红外探头，其中，双波长光源S到检测器D1之间的直线距离为r₁，5mm≤r₁≤10mm，用于敏感外层脑组织的血液动力学变化；双波长光源S到检测器D2之间的直线距离为r₂，30mm≤r₂≤45mm，能够敏感大脑皮质的血液动力学变化，通过检测器记录大脑安静状态下的漫反射光强和大脑处于诱发激励时的漫反射光强，以获得两个不同波长λ₁和λ₂时的光密度变化量的时间序列：

和

和

k为时间，k=1,2,...,N；N为正整数，

表示在双波长光源S到检测器D1之间的直线距离为r₁且波长为λ₁时光密度变化量的时间序列，

表示在双波长光源S到检测器D1之间的直线距离为r₁且波长为λ₂时光密度变化量的时间序列，

表示在双波长光源S到检测器D1之间的直线距离为r₂且波长为λ₁时光密度变化量的时间序列，

表示在双波长光源S到检测器D1之间的直线距离为r₂且波长为λ₂时光密度变化量的时间序列；

步骤二、根据步骤一获得的光密度变化量的时间序列并采用修正朗伯比尔定律获取r₁测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)，r₂测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^F(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^F(k)；

$\mathrm{\Δ}[\mathrm{Hb}{O}_2{]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPE})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\text{\λ}}_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{\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)\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{Hb}{O}_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{Hb}{O}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{H{\mathrm{bO}}_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))},$

其中，ε_HHb(λ₁)为探头光源的波长为λ₁时的消光系数，

为探头光源的波长为λ₂时的消光系数，

DPF为差分路径因子；

步骤三、用x(k)表示步骤二中的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)，用经验模态分解将x(k)分解为N个固态模式函数分量IMF分量，将剩余分量作为最后的IMF分量，则x(k)表示为

$x\left(k\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{c}_i\left(k\right)$

其中，c_i(k)为分解的IMF分量；

步骤四、用d(k)表示步骤二中的Δ[HbO₂]^F(k)或Δ[HHb]^F(k)，d(k)中包含脑功能活动信号r(k)和生理干扰i(k)，即d(k)=r(k)+i(k)，采用线性映射关系，用c_i(k)的线性组合表示i(k)的估计，即

$\hat{i}\left(k\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(k\right)$

其中，

为i(k)的估计，i=1,2,...,N，w_i(k)为第i个IMF分量的权系数；

步骤五、根据步骤二中的d(k)=r(k)+i(k)和

即可推算出脑功能信号的表达式：

$e\left(k\right)=d\left(k\right)-\hat{i}\left(k\right)=r\left(k\right)+[i\left(k\right)-\hat{i}\left(k\right)]$

其中，e(k)为脑功能信号，r(k)为e(k)的脑功能信号估计；

步骤六、利用加权最小二乘算法作为代价函数，求取优化系数w_i(k)，再将求取优化的系数w_i(k)带入步骤五中的

公式中，即可获得脑功能信号e(k)，加权最小二乘算法为：

$J\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{e}^2\left(n\right)$

进一步表示为

$J\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{e}^2\left(n\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{[d\left(n\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(n\right)]}^{\text{2}}$

其中，χ为指数加权因子，χ=0.99;n=1,…k,k为正整数，i=1,2,...,N，N为正整数，求解使J(k)最小的w_i(k)，获得脑功能信号e(k)。

具体实施方式二、本实施方式与具体实施方式一的区别在于：步骤一所述的双波长光源发出的两种波长分别为λ₁=760nm，λ₂=850nm。

具体实施方式三、本实施方式与具体实施方式一的区别在于：步骤一所述的光源S与检测器D1的间距为10mm，发光源S与检测器D2的间距为40mm。

具体实施方式四、本实施方式与具体实施方式一的区别在于：步骤一中光密度变化量的时间序列

和

按如下公式获取：

${\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),$

其中：

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

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

为探头光源的波长为λ₁时，大脑处于诱发激励时检测器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),$

其中：

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

为探头光源的波长为λ₂时，大脑处于诱发激励时检测器D1测得的出射光强；

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

为探头光源的波长为λ₂时，大脑处于诱发激励时检测器D2测得的出射光强。

具体实施方式五、本实施方式与具体实施方式一的区别在于：步骤六的脑功能信号e(k)的获得方法为：

步骤六一、通过最小二乘估计准则表示使脑功能信号e(k)的累计平方误差性能函数J(k)最小，J(k)表示为

$J\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{e}^2\left(n\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{[d\left(n\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(n\right)]}^2$

步骤六二、求解最优系数w_i(k)：

通过对J(k)相对于w_i(k)求导，将求导结果置为0，即

$-2\underset{n=1}{\overset{k}{\mathrm{\Σ}}}\left\{{\mathrm{\χ}}^{k-n}\right[d\left(n\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(n\right)\text{]}{c}_j\left(n\right)\text{}=0}$

由上式得到

$\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{X}^{k-n}d\left(n\right)c_j\left(n\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{X}^{k-n}{c}_i\left(n\right)c_j\left(n\right)$

或

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{R}_{\mathrm{ij}}\left(k\right)w_i\left(k\right)=P_j\left(k\right),j=\mathrm{1,2},...,N$

其中，P_j(k)和R_ij(k)的表达式为

$P_j\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{X}^{k-n}d\left(n\right)c_j\left(n\right)$

$R_{\mathrm{ij}}\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{X}^{k-n}{c}_i\left(n\right)c_j\left(n\right)$

其矩阵形式的表示为

可进一步简化为

R(k)w(k)=p(k)

若矩阵R(k)非奇异，最优系数通过下式计算得到

w^*(k)=R^-1(k)p(k)

其中，w*(k)表示为w(k)的最优解，

R^-1(K)为R(K)的逆矩阵，

$w\left(k\right)=\left[\begin{array}{c}{w}_1\left(k\right)\\ w_2\left(k\right)\\ \·\\ \·\\ \·\\ w_N\left(k\right)\end{array}\right],$

$p\left(k\right)=\left[\begin{array}{c}{p}_1\left(k\right)\\ p_2\left(k\right)\\ \·\\ \·\\ \·\\ p_N\left(k\right)\end{array}\right],$

$c\left(k\right)=\left[\begin{array}{c}{c}_1\left(k\right)\\ c_2\left(k\right)\\ \·\\ .\\ .\\ c_N\left(k\right)\end{array}\right];$

步骤六三、求解脑功能信号e(k)：

e(k)=d(k)-c^T(k)w^*(k)，

其中c^T(k)表示的是c(k)的转置矩阵，w^*(k)表示求解的最优系数向量。

通过单光源双检测器的探头结构，光源采用双波长光源λ₁=760nm，λ₂=850nm，光源S到检测器D1的直线距离即光源检测器间距为10mm，光源S到检测器D2的直线距离即光源检测器间距为40mm。光源检测器间距为近红外光探测深度的两倍，这样设置能够使D2检测的近红外光可有效穿入大脑皮层，D1检测的近红外光仅穿头外层脑组织。将获得的光密度变化通过修正朗伯比尔定律转变为氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)、Δ[HbO₂]^F(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)、Δ[HHb]^F(k)。通过经验模态分解算法对近端血液动力学变化Δ[HbO₂]^N(k)或Δ[HHb]^N(k)分解为固态模式函数分量。将IMF分量进行线性组合估计Δ[HbO₂]^F(k)或Δ[HHb]^F(k)中的生理干扰，通过自适应滤波算法将构建脑功能活动信号e(k)。通过最小二乘估计准则求解使脑功能信号e(k)的累计平方误差性能函数J(k)最小，e(k)即是通过自适应滤波剔除生理干扰的脑功能活动信号。

Claims (5)

1.基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法，其特征在于：它包括以下步骤：

步骤一、在待测脑组织的头皮表面放置由双波长光源S和检测器D1和D2构成的近红外探头，其中，双波长光源S到检测器D1之间的直线距离为r₁，5mm≤r₁≤10mm，用于敏感外层脑组织的血液动力学变化；双波长光源S到检测器D2之间的直线距离为r₂，30mm≤r₂≤45mm，能够敏感大脑皮质的血液动力学变化，通过检测器记录大脑安静状态下的漫反射光强和大脑处于诱发激励时的漫反射光强，以获得两个不同波长λ₁和λ₂时的光密度变化量的时间序列：

和，

和

，k为时间，k=1,2,...,N；N为正整数，表示在双波长光源S到检测器D1之间的直线距离为r₁且波长为λ₁时光密度变化量的时间序列，

表示在双波长光源S到检测器D1之间的直线距离为r₁且波长为λ₂时光密度变化量的时间序列，

表示在双波长光源S到检测器D2之间的直线距离为r₂且波长为λ₁时光密度变化量的时间序列，

表示在双波长光源S到检测器D2之间的直线距离为r₂且波长为λ₂时光密度变化量的时间序列；

步骤二、根据步骤一获得的光密度变化量的时间序列并采用修正朗伯比尔定律获取r₁测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)，r₂测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^F(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^F(k)；

$\mathrm{\Δ}[\mathrm{Hb}{O}_2{]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right){\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPE})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\text{\λ}}_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{\mathrm{\ΔOD}}_{{\mathrm{\λ}}_2}^N\left(k\right)\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{Hb}{O}_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{Hb}{O}_2}\left({\mathrm{\λ}}_2\right){\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)-{\mathrm{\ϵ}}_{H{\mathrm{bO}}_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))},$

其中，ε_HHb(λ₁)为探头光源的波长为λ₁时的消光系数，

为探头光源的波长为λ₂时的消光系数，

DPF为差分路径因子；

步骤三、用x(k)表示步骤二中的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)，用经验模态分解将x(k)分解为N个固态模式函数分量IMF分量，将剩余分量作为最后的IMF分量，则x(k)表示为

$\left(k\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{c}_i\left(k\right)$

其中，c_i(k)为分解的IMF分量；

步骤四、用d(k)表示步骤二中的Δ[HbO₂]^F(k)或2[HHb]^F(k)，d(k)中包含脑功能活动信号r(k)和生理干扰i(k)，即d(k)=r(k)+i(k)，采用线性映射关系，用c_i(k)的线性组合表示i(k)的估计，即

$\hat{i}\left(k\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(k\right)$

其中，为i(k)的估计，i=1,2,...,N，wi(k)为第i个IMF分量的权系数；

步骤五、根据步骤二中的d(k)=r(k)+i(k)和

即可推算出脑功能活动信号估计的表达式：

$e\left(k\right)=d\left(k\right)-\hat{i}\left(k\right)+[i\left(k\right)-\hat{i}\left(k\right)]$

其中，e(k)为脑功能活动信号估计；

步骤六、利用加权最小二乘算法作为代价函数，求取优化系数w_i(k)，再将求取优化的系数w_i(k)带入步骤五中的

公式中，即可获得脑功能活动信号估计e(k)，加权最小二乘算法为：

$J\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{e}^2\left(n\right)$

进一步表示为

$J\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{e}^2\left(n\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{[d\left(\text{n}\right)\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(n\right)\text{]}}^2$

其中，χ为指数加权因子，χ=0.99;n=1,…k,k为正整数，i=1,2,...,N，N为正整数，求解使J(k)最小的w_i(k)，获得脑功能活动信号估计e(k)。

2.根据权利要求1所述的基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法，其特征在于：步骤一所述的双波长光源发出的两种波长分别为λ₁=760nm，λ₂=850nm。

3.根据权利要求1所述的基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法，其特征在于：步骤一所述的光源S与检测器D1的间距为10mm，光源S与检测器D2的间距为40mm。

4.根据权利要求1所述的基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法，其特征在于：步骤一中光密度变化量的时间序列和

按如下公式获取：

$\mathrm{\Δ}{\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),$

其中：

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

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

为探头光源的波长为λ₁时，大脑处于诱发激励时检测器D1测得的出射光强；

为探头光源的波长为λ₁时，大脑处于诱发激励时检测器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),$

其中：

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

为探头光源的波长为λ₂时，大脑处于诱发激励时检测器D1测得的出射光强；

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

为探头光源的波长为λ₂时，大脑处于诱发激励时检测器D2测得的出射光强。

5.根据权利要求1所述的基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法，其特征在于：步骤六的脑功能活动信号估计e(k)的获得方法为：

步骤六一、通过最小二乘估计准则表示使脑功能活动信号估计e(k)的累计平方误差性能函数J(k)最小，J(k)表示为

$J\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{e}^2\left(n\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{[d\left(n\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(n\right)]}^2$

步骤六二、求解最优系数w_i(k)：

通过对J(k)相对于w_i(k)求导，将求导结果置为0，即

$-2\underset{n=1}{\overset{k}{\mathrm{\Σ}}}\left\{{\mathrm{\χ}}^{k-n}\right[d\left(n\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)c_i\left(n\right)\left]c_j\left(n\right)\right\}=0$

由上式得到

$\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}d\left(n\right)c_j\left(n\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_i\left(k\right)\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{c}_i\left(n\right)c_j\left(n\right)$

或

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{R}_{\mathrm{ij}}\left(k\right)w_i\left(k\right)=P_j\left(k\right),j=\mathrm{1,2},...,N$

其中，P_j(k)和R_ij(k)的表达式为

$P_j\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}d\left(n\right)c_j\left(n\right)$

$R_{\mathrm{ij}}\left(k\right)=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-n}{c}_i\left(n\right)c_j\left(n\right)$

其矩阵形式的表示为

可进一步简化为

R(k)w(k)=p(k)

若矩阵R(k)非奇异，最优系数通过下式计算得到

w^*(k)=R^-1(k)p(k)

其中，w*(k)表示为w(k)的最优解，

R^-1(K)为R(K)的逆矩阵，

$w\left(k\right)=\left[\begin{array}{c}{w}_1\left(k\right)\\ w_2\left(k\right)\\ .\\ .\\ .\\ w_N\left(k\right)\end{array}\right],$

$p\left(k\right)=\left[\begin{array}{c}{p}_1\left(k\right)\\ p_2\left(k\right)\\ .\\ .\\ .\\ P_N\left(k\right)\end{array}\right],$

$c\left(k\right)=\left[\begin{array}{c}{c}_1\left(k\right)\\ c_2\left(k\right)\\ .\\ .\\ .\\ c_N\left(k\right)\end{array}\right];$

步骤六三、求解脑功能活动信号估计e(k)：

e(k)=d(k)-c^T(k)w^*(k)，

其中c^T(k)表示的是c(k)的转置矩阵，w^*(k)表示求解的最优系数向量。

Images