CN102512142A - 基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法 - Google Patents

Abstract

基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法，涉及一种脑功能活动信号提取方法，为了解决采用自适应滤波无法有效去除脑功能检测时多种生理干扰且存在需要借助额外设备的问题，它包括如下步骤：步骤一：在待测脑组织头皮放置由双波长光源S、检测器D1和检测器D2构成的近红外探头，获得D1和D2检测的光密度变化量；步骤二：采用修正朗伯比尔定律获取光源S和检测器测得氧合血红蛋白浓度变化量和还原血红蛋白浓度变化量；步骤三：构建自适应滤波脑功能活动信号函数；步骤四：利用最小二乘估计准则求解优化滤波器系数向量，进而求解脑功能活动信号。本发明用于脑功能活动信号提取。它有效去除脑功能检测的生理性干扰且不借助额外设备。

Description

基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法

技术领域

本发明涉及一种脑功能活动信号提取方法，特别涉及一种基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法。

背景技术

近红外光谱技术(NIRS)能提供脑功能活动过程中的大脑皮层血氧代谢信息——氧合血红蛋白浓度变化(Δ[HbO₂])和还原血红蛋白浓度变化(Δ[HHb])，可用于脑功能活动的检测。与其它的脑功能检测方法如：功能性核磁共振、脑磁描记法、正电子放射层扫描术、以及脑电描记法相比，近红外光谱技术具有使用方便、易实施、时间分辨率高、安全、便宜等优点。然而，利用近红外光谱技术进行诱发激励时脑功能活动的检测，会受到人体的生理活动如心脏跳动、呼吸、低频振荡、超低频振荡的影响，称之为生理干扰。这种生理干扰不但出现在头皮、颅骨和脑脊液等外层脑组织中，也出现在脑灰质和脑白质等深层脑组织中，严重影响了脑功能活动信号的准确提取。

在近红外脑功能检测时，除心动周期和呼吸引起的生理干扰外，人体的低频振荡和超低频振荡等干扰信号与脑功能信号的频带严重重叠，很难用一般的滤波技术去除干扰。但采用自适应滤波处理技术，不但可以不必准确知道干扰频率，而且能自动跟踪频率的漂移。但采用自适应滤波技术，参考信号源的选取十分重要。为了选取参考信号，可以借助于脉搏血氧计测量心动周期信号或用呼吸计，肺活量计以及二氧化碳分析仪等仪器测量呼吸信号。但由于该方法需要借助于额外的设备获取参考信号，并且难以获得所有的干扰信号，在实际应用中受到了一定的限制。

发明内容

本发明目的是为了解决采用自适应滤波无法有效的去除脑功能检测时的多种生理干扰并且存在需要借助额外设备的问题，提供了一种基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法。

本发明的基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法，它包括如下步骤：

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

和

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

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

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

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

步骤二：根据步骤一获得的光密度变化量的时间序列采用修正朗伯比尔定律获取双波长光源S和检测器D1测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)，以及双波长光源S和检测器D2测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^F(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^F(k)，Δ[HbO₂]^N(k)表示双波长光源S到检测器D1之间的直线距离为r₁时的氧合血红蛋白浓度变化量的时间序列，Δ[HHb]^N(k)表示双波长光源S到检测器D1之间的直线距离为r₁时的还原血红蛋白浓度变化量的时间序列，Δ[HbO₂]^F(k)表示双波长光源S到检测器D1之间的直线距离为r₂时的氧合血红蛋白浓度变化量的时间序列，Δ[HHb]^F(k)表示双波长光源S到检测器D1之间的直线距离为r₂时的还原血红蛋白浓度变化量的时间序列，

$\mathrm{\Δ}{\left[{\mathrm{HbO}}_2\right]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right)\mathrm{\Δ}{\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{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\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{HbO}}_2\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_1}^N\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{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_2}^N\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为差分路径因子；

步骤三：利用步骤二获得的Δ[HbO₂]^N(k)、Δ[HHb]^N(k)、Δ[HbO₂]^F(k)和Δ[HHb]^F(k)构建自适应滤波的脑功能活动信号，脑功能活动信号表示为

e(k)＝d(k)-x^T(k)w(k)，

其中，d(k)表示由双波长光源S和检测器D2解算的Δ[HbO₂]^F(k)或Δ[HHb]^F(k)，包含生理干扰和诱发的血液动力学变化，k为取样时刻；x(k)是一个列向量，表示为x(k)＝[x(t)x(k-1)…x(k-N)]^T，N为滤波器阶数，x(k-N)为信号x(k)的N个单元延长；x(k)表示由双波长光源S和检测器D2获得的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)，主要由生理干扰构成，作为自适应滤波的参考信号；w(k)＝[w₀(k)w₁(k)...w_N(k)]^T是自适应滤波器的系数向量；

步骤四：使脑功能活动信号e(k)的累计平方误差性能函数J(k)最小，求解w(k)最优的滤波器系数向量：

$J\left(k\right)=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}{[d\left(i\right)-x^T\left(i\right)w\left(k\right)]}^2$

其中，χ为指数加权因子，χ＝0.99，i＝0，1，…，k，k为正整数；

通过对J(k)相对于w(k)求导，则有

$\frac{\∂J\left(k\right)}{\∂w\left(k\right)}=-2\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}x\left(i\right){[d\left(i\right)-x^T\left(i\right)w\left(k\right)]}^2$

令上式等于零，能够得到最优的滤波器系数向量w^*(k)的表达式为

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

式中，R(k)为参考信号x(k)的确定性相关矩阵，p(k)为参考信号x(k)与期望信号d(k)之间的确定性互相关向量，表示为

$R\left(k\right)=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}x\left(i\right)x^T\left(i\right)$

$p\left(k\right)=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}x\left(i\right)d\left(i\right);$

步骤五：求解脑功能活动信号e(k)：

e(k)＝d(k)-x^T(k)w^*(k)

本发明的优点在于当近红外光无法达到大脑皮质，利用双波长光源S与近端检测器D1测量结果主要包含生理干扰；当近红外光穿透外层脑组织达到大脑皮质，利用双波长光源S与远端检测器D2测量结果包含脑功能活动信号以及生理干扰。为了借助自适应滤波技术消除脑功能检测中的生理性干扰，将双波长光源S与远端检测器D2测量结果视为期望信号，双波长光源S与近端检测器D1测量结果视为参考信号。自适应滤波算法相当于将利用双波长光源S与近端检测器D1测量的外层脑组织血液动力学参数与期望信号中的生理干扰进行线性映射，进而估计生理干扰并将其从双波长光源S与远端检测器D2测量结果中消除。这样有效的删除了脑功能检测时的多种生理干扰，而不需要借助额外测试设备，例如脉搏血氧计、呼吸计、肺活量计或二氧化碳分析仪等仪器，仅仅通过多距测量探头即可实现对脑功能信号动态跟踪。

附图说明

图1是本发明的由双波长光源S、检测器D1和检测器D2构成的近红外探头结构示意图，其中a表示头皮，b表示颅骨，c表示脑脊液，d表示脑灰质，e表示脑白质。

具体实施方式

具体实施方式一：结合图1说明本实施方式，本发明的基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法，它包括如下步骤：

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

和

和k为时间，k＝1，2，...，N，N为正整数；

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

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

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

步骤二：根据步骤一获得的光密度变化量的时间序列采用修正朗伯比尔定律获取双波长光源S和检测器D1测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)，以及双波长光源S和检测器D2测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^F(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^F(k)，Δ[HbO₂]^N(k)表示双波长光源S到检测器D1之间的直线距离为r₁时的氧合血红蛋白浓度变化量的时间序列，Δ[HHb]^N(k)表示双波长光源S到检测器D1之间的直线距离为r₁时的还原血红蛋白浓度变化量的时间序列，Δ[HbO₂]^F(k)表示双波长光源S到检测器D1之间的直线距离为r₂时的氧合血红蛋白浓度变化量的时间序列，Δ[HHb]^F(k)表示双波长光源S到检测器D1之间的直线距离为r₂时的还原血红蛋白浓度变化量的时间序列，

$\mathrm{\Δ}{\left[{\mathrm{HbO}}_2\right]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right)\mathrm{\Δ}{\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{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\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{HbO}}_2\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_1}^N\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{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_2}^N\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为差分路径因子；

步骤三：利用步骤二获得的Δ[HbO₂]^N(k)、Δ[HHb]^N(k)、Δ[HbO₂]^F(k)和Δ[HHb]^F(k)构建自适应滤波的脑功能活动信号，脑功能活动信号表示为

e(k)＝d(k)-x^T(k)w(k)，

其中，d(k)表示由双波长光源S和检测器D2解算的Δ[HbO₂]^F(k)或Δ[HHb]^F(k)，包含生理干扰和诱发的血液动力学变化，k为取样时刻；x(k)是一个列向量，表示为x(k)＝[x(t)x(k-1)…x(k-N)]^T，N为滤波器阶数，x(k-N)为信号x(k)的N个单元延长；x(k)表示由双波长光源S和检测器D2获得的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)，主要由生理干扰构成，作为自适应滤波的参考信号；w(k)＝[w₀(k)w₁(k)...w_N(k)]^T是自适应滤波器的系数向量；

步骤四：使脑功能活动信号e(k)的累计平方误差性能函数J(k)最小，求解w(k)最优的滤波器系数向量：

$J\left(k\right)=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}{[d\left(i\right)-x^T\left(i\right)w\left(k\right)]}^2$

其中，χ为指数加权因子，χ＝0.99，i＝0，1，…，k，k为正整数；

通过对J(k)相对于w(k)求导，则有

$\frac{\∂J\left(k\right)}{\∂w\left(k\right)}=-2\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}x\left(i\right){[d\left(i\right)-x^T\left(i\right)w\left(k\right)]}^2$

令上式等于零，能够得到最优的滤波器系数向量w^*(k)的表达式为

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

式中，R(k)为参考信号x(k)的确定性相关矩阵，p(k)为参考信号x(k)与期望信号d(k)之间的确定性互相关向量，表示为

$R\left(k\right)=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}x\left(i\right)x^T\left(i\right)$

$p\left(k\right)=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}x\left(i\right)d\left(i\right);$

步骤五：求解脑功能活动信号e(k)：

e(k)＝d(k)-x^T(k)w^*(k)。

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

具体实施方式三：本实施方式与具体实施方式一不同的是步骤一所述的双波长光源S与检测器D1的直线距离为10mm，双波长光源S与检测器D2的直线距离为40mm。

具体实施方式四：本实施方式与具体实施方式一不同的是步骤一中光密度变化量的时间序列和

按如下公式获取：

$\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{\Δ}{\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{\Δ}{\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{\Δ}{\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测得的出射光强。

大脑活动过程中，人的思维活动以及肢体运动会引起脑灰质d内血氧含量变化，并且可以认为外层脑组织的血液动力学参数不发生变化。然而，由于在实际测量中存在生理干扰，D1检测器检测到的光密度变化是波动的信号，来源于人体的生理活动引起的外层组织的血液动力学变化；检测器D2测得的光密度变化是脑灰质的血液动力学变化以及各层生理干扰共同作用的结果。因此，通常情况下双波长光源S和检测器D1测得的血液动力学变化与双波长光源S和检测器D2中的生理干扰具有很强的相关性。为了借助自适应滤波技术消除脑功能检测中的生理性干扰，将双波长光源S和检测器D2测量结果视为期望信号，双波长光源S和检测器D1测量结果视为参考信号。自适应滤波相当于将外层脑组织的血液动力学参数与期望信号中的生理干扰进行线性映射，进而估计脑功能活动信号。滤波的过程就是对整个时间序列逐点的进行线性映射，即使外层组织的血液动力学参数变化是非稳态的，只要这种变化的速度相对于自适应滤波收敛速度相比是缓慢的，依然可以提取脑功能信号。

Claims (4)

1.基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法，其特征在于它包括如下步骤：

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

和

和

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

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

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

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

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

步骤二：根据步骤一获得的光密度变化量的时间序列采用修正朗伯比尔定律获取双波长光源S和检测器D1测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^N(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^N(k)，以及双波长光源S和检测器D2测得的氧合血红蛋白浓度变化量的时间序列Δ[HbO₂]^F(k)和还原血红蛋白浓度变化量的时间序列Δ[HHb]^F(k)，Δ[HbO₂]^N(k)表示双波长光源S到检测器D1之间的直线距离为r₁时的氧合血红蛋白浓度变化量的时间序列，Δ[HHb]^N(k)表示双波长光源S到检测器D1之间的直线距离为r₁时的还原血红蛋白浓度变化量的时间序列，Δ[HbO₂]^F(k)表示双波长光源S到检测器D1之间的直线距离为r₂时的氧合血红蛋白浓度变化量的时间序列，Δ[HHb]^F(k)表示双波长光源S到检测器D1之间的直线距离为r₂时的还原血红蛋白浓度变化量的时间序列，

$\mathrm{\Δ}{\left[{\mathrm{HbO}}_2\right]}^N\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right)\mathrm{\Δ}{\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{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\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{HbO}}_2\right]}^F\left(k\right)=\frac{({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_2}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{\mathrm{HHb}}\left({\mathrm{\λ}}_2\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_1}^N\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{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_1}^N\left(k\right)/\mathrm{DPF})-({\mathrm{\ϵ}}_{{\mathrm{HbO}}_2}\left({\mathrm{\λ}}_1\right)\mathrm{\Δ}{\mathrm{OD}}_{{\mathrm{\λ}}_2}^N\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为差分路径因子；

步骤三：利用步骤二获得的Δ[HbO₂]^N(k)、Δ[HHb]^N(k)、Δ[HbO₂]^F(k)和Δ[HHb]^F(k)构建自适应滤波的脑功能活动信号，脑功能活动信号表示为

e(k)＝d(k)-x^T(k)w(k)，

其中，d(k)表示由双波长光源S和检测器D2解算的Δ[HbO₂]^F(k)或Δ[HHb]^F(k)，包含生理干扰和诱发的血液动力学变化，k为取样时刻；x(k)是一个列向量，表示为x(k)＝[x(t)x(k-1)…x(k-N)]^T，N为滤波器阶数，x(k-N)为信号x(k)的N个单元延长；x(k)表示由双波长光源S和检测器D2获得的Δ[HbO₂]^N(k)或Δ[HHb]^N(k)，主要由生理干扰构成，作为自适应滤波的参考信号；w(k)＝[w₀(k)w₁(k)...w_N(k)]^T是自适应滤波器的系数向量；

步骤四：使脑功能活动信号e(k)的累计平方误差性能函数J(k)最小，求解w(k)最优的滤波器系数向量：

$J\left(k\right)=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}{[d\left(i\right)-x^T\left(i\right)w\left(k\right)]}^2$

其中，χ为指数加权因子，χ＝0.99，i＝0，1，…，k，k为正整数；

通过对J(k)相对于w(k)求导，则有

$\frac{\∂J\left(k\right)}{\∂w\left(k\right)}=-2\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}x\left(i\right){[d\left(i\right)-x^T\left(i\right)w\left(k\right)]}^2$

令上式等于零，能够得到最优的滤波器系数向量w^*(k)的表达式为

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

式中，R(k)为参考信号x(k)的确定性相关矩阵，p(k)为参考信号x(k)与期望信号d(k)之间的确定性互相关向量，表示为

$R\left(k\right)=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}x\left(i\right)x^T\left(i\right)$

$p\left(k\right)=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\χ}}^{k-i}x\left(i\right)d\left(i\right);$

步骤五：求解脑功能活动信号e(k)：

e(k)＝d(k)-x^T(k)w^*(k)。

2.根据权利要求1所述的基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法，其特征在于步骤一所述的双波长光源S发出的两种波长分别为λ₁＝760nm，λ₂＝850nm。

3.根据权利要求1所述的基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法，其特征在于步骤一所述的双波长光源S与检测器D1的直线距离为10mm，双波长光源S与检测器D2的直线距离为40mm。

4.根据权利要求1所述的基于多距测量方法的递归最小二乘自适应滤波近红外脑功能活动信号提取方法，其特征在于步骤一中光密度变化量的时间序列

(k)和(k)、按如下公式获取：

$\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{\Δ}{\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{\Δ}{\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{\Δ}{\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测得的出射光强。

Images