首页> 中国专利> 基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法

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

摘要

基于多距测量方法的经验模态分解优化算法的脑功能信号提取方法,涉及脑功能信号提取方法。它解决了当脑组织非均匀性严重时现有技术检测脑功能活动过程中氧合血红蛋白浓度变化δ[hbo

著录项

  • 公开/公告号CN102525422A

    专利类型发明专利

  • 公开/公告日2012-07-04

    原文格式PDF

  • 申请/专利权人 哈尔滨工业大学;

    申请/专利号CN201110442356.7

  • 申请日2011-12-26

  • 分类号A61B5/00;A61B5/1455;

  • 代理机构哈尔滨市松花江专利商标事务所;

  • 代理人牟永林

  • 地址 150001 黑龙江省哈尔滨市南岗区西大直接92号

  • 入库时间 2023-12-18 05:38:43

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2016-02-17

    未缴年费专利权终止 IPC(主分类):A61B5/00 授权公告日:20140402 终止日期:20141226 申请日:20111226

    专利权的终止

  • 2014-04-02

    授权

    授权

  • 2012-09-05

    实质审查的生效 IPC(主分类):A61B5/00 申请日:20111226

    实质审查的生效

  • 2012-07-04

    公开

    公开

说明书

技术领域

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

背景技术

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

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

发明内容

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

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

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

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

Δ[HbO2]N(k)=(ϵHHb(λ1)ΔODλ2N(k)/DPF)-(ϵHHb(λ2)ΔODλ1N(k)/DPE)r1(ϵHbO2(λ2)ϵHHb(λ1)-ϵHbO2(λ1)ϵHHb(λ2)),

Δ[HHb]N(k)=(ϵHbO2(λ2)ΔODλ1N(k)/DPF)-(ϵHbO2(λ1)ΔODλ2N(k)/DPF)r1(ϵHbO2(λ2)ϵHHb(λ1)-ϵHbO2(λ1)ϵHHb(λ2)),

Δ[HbO2]F(k)=(ϵHHb(λ1)ΔODλ2F(k)/DPF)-(ϵHHb(λ2)ΔODλ1F(k)/DPF)r2(ϵHbO2(λ2)ϵHHb(λ1)-ϵHbO2(λ1)ϵHHb(λ2)),

Δ[HHb]F(k)=(ϵHbO2(λ2)ΔODλ1F(k)/DPF)-(ϵHbO2(λ1)ΔODλ2F(k)/DPF)r2(ϵHbO2(λ2)ϵHHb(λ1)-ϵHbO2(λ1)ϵHHb(λ2)),

其中,εHHb1)为探头光源的波长为λ1时的消光系数,

为探头光源的波长为λ2时的消光系数,

DPF为差分路径因子;

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

x(k)=Σi=1Nci(k)

其中,ci(k)为分解的IMF分量;

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

i^(k)=Σi=1Nwi(k)ci(k)

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

步骤五、根据步骤二中的d(k)=r(k)+i(k)和即可推算出 脑功能信号的表达式:

e(k)=d(k)-i^(k)=r(k)+[i(k)-i^(k)]

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

步骤六、利用加权最小二乘算法作为代价函数,求取优化系数wi(k),再 将求取优化的系数wi(k)带入步骤五中的公式 中,即可获得脑功能信号e(k),加权最小二乘算法为:

J(k)=Σn=1kχk-ne2(n)

进一步表示为

J(k)=Σn=1kχk-ne2(n)=Σn=1kχk-n[d(n)-Σi=1Nwi(k)ci(n)]2

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

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

附图说明

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

具体实施方式

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

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

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

Δ[HbO2]N(k)=(ϵHHb(λ1)ΔODλ2N(k)/DPF)-(ϵHHb(λ2)ΔODλ1N(k)/DPF)r1(ϵHbO2(λ2)ϵHHb(λ1)-ϵHbO2(λ1)ϵHHb(λ2)),

Δ[HHb]N(k)=(ϵHbO2(λ2)ΔODλ1N(k)/DPF)-(ϵHbO2(λ1)ΔODλ2N(k)/DPF)r1(ϵHbO2(λ2)ϵHHb(λ1)-ϵHbO2(λ1)ϵHHb(λ2)),

Δ[HbO2]F(k)=(ϵHHb(λ1)ΔODλ2F(k)/DPF)-(ϵHHb(λ2)ΔODλ1F(k)/DPF)r2(ϵHbO2(λ2)ϵHHb(λ1)-ϵHbO2(λ1)ϵHHb(λ2)),

Δ[HHb]F(k)=(ϵHbO2(λ2)ΔODλ1F(k)/DPF)-(ϵHbO2(λ1)ΔODλ2F(k)/DPF)r2(ϵHbO2(λ2)ϵHHb(λ1)-ϵHbO2(λ1)ϵHHb(λ2)),

其中,εHHb1)为探头光源的波长为λ1时的消光系数,

为探头光源的波长为λ2时的消光系数,

DPF为差分路径因子;

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

x(k)=Σi=1Nci(k)

其中,ci(k)为分解的IMF分量;

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

i^(k)=Σi=1Nwi(k)ci(k)

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

步骤五、根据步骤二中的d(k)=r(k)+i(k)和即可推算出 脑功能信号的表达式:

e(k)=d(k)-i^(k)=r(k)+[i(k)-i^(k)]

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

步骤六、利用加权最小二乘算法作为代价函数,求取优化系数wi(k),再 将求取优化的系数wi(k)带入步骤五中的公式 中,即可获得脑功能信号e(k),加权最小二乘算法为:

J(k)=Σn=1kχk-ne2(n)

进一步表示为

J(k)=Σn=1kχk-ne2(n)=Σn=1kχk-n[d(n)-Σi=1Nwi(k)ci(n)]2

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

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

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

具体实施方式四、本实施方式与具体实施方式一的区别在于:步骤一中 光密度变化量的时间序列和按如下公式获取:

ΔODλ1N(k)=logIbaseN(λ1)/IstimN(λ1),

ΔODλ1F(k)=logIbaseF(λ1)/IstimF(λ1),

其中:为探头光源的波长为λ1时,大脑处于安静状态下时检测器 D1测得的出射光强;为探头光源的波长为λ1时,大脑处于安静状态 下时检测器D2测得的出射光强;为探头光源的波长为λ1时,大脑处 于诱发激励时检测器D1测得的出射光强;为探头光源的波长为λ1时, 大脑处于诱发激励时检测器D2测得的出射光强。

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

ΔODλ2N(k)=logIbaseN(λ2)/IstimN(λ2),

ΔODλ2F(k)=logIbaseF(λ2)/IstimF(λ2),

其中:为探头光源的波长为λ2时,大脑处于安静状态下时检测 器D1测得的出射光强,为探头光源的波长为λ2时,大脑处于诱发激 励时检测器D1测得的出射光强;为探头光源的波长为λ2时,大脑处 于安静状态下时检测器D2测得的出射光强,为探头光源的波长为λ2时,大脑处于诱发激励时检测器D2测得的出射光强。

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

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

J(k)=Σn=1kχk-ne2(n)=Σn=1kχk-n[d(n)-Σi=1Nwi(k)ci(n)]2

步骤六二、求解最优系数wi(k):

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

-2Σn=1k{χk-n[d(n)-Σi=1Nwi(k)ci(n)]cj(n)}=0

由上式得到

Σn=1kχk-nd(n)cj(n)=Σi=1Nwi(k)Σn=1kχk-nci(n)cj(n)

Σi=1NRij(k)wi(k)=Pj(k),j=1,2,...,N

其中,Pj(k)和Rij(k)的表达式为

Pj(k)=Σn=1kχk-nd(n)cj(n)

Rij(k)=Σn=1kχk-nci(n)cj(n)

其矩阵形式的表示为

可进一步简化为

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(k)=w1(k)w2(k)...wN(k),

p(k)=p1(k)p2(k)...pN(k),

c(k)=c1(k)c2(k)...cN(k);

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

e(k)=d(k)-cT(k)w*(k),

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

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

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号