扩展的卡尔曼滤波在近红外光谱提取脑血流信号中的研究
刘颂阳1, 刘光达1, 刘卓娅2, 邱吉庆3,*, 蔡靖1, 朱展鹏3, 张程3, 齐远3, 张尚1
1. 吉林大学, 吉林 长春 130061
2. 吉林大学珠海学院, 广东 珠海 519041
3. 吉林大学第一医院, 吉林 长春 130061
*通讯联系人 e-mail: songyangl15@mails.jlu.edu.cn

作者简介: 刘颂阳, 1986年生, 吉林大学仪器科学与电气工程学院博士研究生 e-mail: 258727076@qq.com

摘要

脑血流中的血红蛋白有两种: 氧合血红蛋白(HbO2)和还原血红蛋白(HbR)。 这两种血红蛋白在脑血流中浓度的变化可以反应脑部神经活动, 提取其浓度变化信号可以为如癫痫病灶定位、 抑郁等相关疾病的诊断和治疗提供依据和参考。 目前, 使用近红外光谱提取脑血流信号的算法有EEMD-ICA法主成分分析法(PCA)、 独立成分分析法(ICA)、 相干平均法、 自适应滤波等, 这些算法在对近红外脑神经活动信号提取时都有各自的特点和优势, 但都重视如呼吸、 眼动等各种生理干扰, 忽视了测量过程中符合高斯分布的测量干扰, 如仪器精密度、 信号传输中的串扰等。 为了提取脑血流中氧合血红蛋白(HbO2)和还原血红蛋白(HbR)浓度变化信号, 设计了功能性近红外光谱(fNIRS)的脑血流参数采集装置, 选择波长为750和830 nm的二极管近红外光源采集脑部血流变化信号, 采用扩展的卡尔曼滤波(EKF)算法, 把生理干扰和测量干扰建立对应的数学模型, 使用基于误差平方和最小的原理进行递归计算, 通过对下一时刻系统的初步状态估计以及测量得出的反馈相结合, 得到该时刻无限逼近真实值的状态估计, 结合修正的朗伯比尔定律(Lambert-Beer law), 将光密度信号的变化转换为氧合血红蛋白(HbO2)和还原血红蛋白(HbR)浓度变化信号。 结果表明: 所提方法可以有效去除符合高斯分布的测量干扰, 在Valsava实验和视觉诱发实验中, 可以提取出脑血流中氧合血红蛋白(HbO2)和还原血红蛋白(HbR)浓度变化曲线, 和主流的EEMD提取脑信号算法比对其RMSE值提高了0.96%, r值提高了0.6%, 表明提出的方法有一定的优越性。 所提方法为相关脑部疾病诊断等提供了有效的脑神经活动探测方法。

关键词: EKF算法; fNIRS; Valsava实验; 视觉诱发; 血红蛋白
中图分类号:TP274 文献标志码:A
Research on Extended Kalman Filter in Extracting Cerebral Blood Flow Signals by Near Infrared Spectroscopy
LIU Song-yang1, LIU Guang-da1, LIU Zhuo-ya2, QIU Ji-qing3,*, CAI Jing1, ZHU Zhan-peng3, ZHANG Cheng3, QI Yuan3, ZHANG Shang1
1. Jilin University, Changchun 130061, China
2. Zhuhai College of Jilin University, Zhuhai 519041, China
3. The First Hospital of Jilin University, Changchun 130061, China
*Corresponding author
Abstract

There are two types of hemoglobin in the cerebral blood stream: oxygenated hemoglobin (HbO2) and reduced hemoglobin (HbR). The changes in the concentration of these two hemoglobins in the cerebral blood flow can reflect the neural activity in the brain. Extracting the signals of concentration changes can provide basis and reference for the diagnosis and treatment of related diseases such as epilepsy focus localization and depression. At present, algorithms for extracting cerebral blood flow signals using near-infrared spectroscopy include the EEMD-ICA method principal component analysis (PCA) , independent component analysis (ICA), the coherent averaging method, Adaptive filtering, etc. The above algorithms have their own characteristics and advantages in the extraction of near-infrared brain neural activity signals. However, the above methods all pay attention to various physiological interferences such as respiration and eye movement and ignore measurement interferences that conform to Gaussian distribution during measurements, such as instrument precision and crosstalk in signal transmission. In order to extract signals of changes in the concentration of oxygenated hemoglobin (HbO2) and reduced hemoglobin (HbR) in cerebral blood flow, a functional near infrared spectroscopy (fNIRS) cerebral blood flow parameter acq uisition device is designed in this article. In the device, a light source Diode near-infrared light sources with wavelengths of 750 and 830 nm were selected to collect brain blood flow changes. The extended Klaman Filter (EKF) algorithm was used to establish a corresponding mathematical model of physiological interference and measurement interference. Perform recursive calculation with the minimum principle, and combine the initial state estimation of the system at the next moment with the measured feedback to obtain a state estimate of infinitely close to the real value at that moment.), The change of the optical density signal is converted into a signal of change in oxygenated hemoglobin (HbO2) and reduced hemoglobin (HbR) concentration. The results show that the method proposed in this paper can effectively remove the measurement interference that conforms to the Gaussian distribution. In the Valsava experiment and the visual evoked experiment, the curve of changes in the concentration of oxygenated hemoglobin (HbO2) and reduced hemoglobin (HbR) in the cerebral blood flow can be extracted. Compared with the mainstream EEMD algorithm for extracting brain signals, its RMSE value is increased by 0.96%, and r value is increased by 0.6%, which indicates that the proposed method has certain advantages. The method proposed in this paper provides an effective method for detecting neural activity in related brain diseases.

Keyword: EKF algorithm; fNIRS; Valsava experiment; Visual induction; Hemoglobin
引言

大脑作为人体的神经中枢, 是人类高级认知功能的载体。 近年来, 用于研究大脑运作机制、 疾病诊断的脑功能成像方法得到了长足的发展, 脑功能成像技术主流的方法有脑电图(electro encephalo gram, EEG)、 脑磁图(magneto encephalo graphy, MEG)、 功能磁共振成像(functional magnetic resonance imaging, FMRI)、 功能性近红外光谱技术(functional near infrared spectroscopy, FNIRS)等。 在上述方法中, 功能性近红外光谱技术(FNIRS)因具有较高的空间分辨率、 适中的时间分辨率、 无创费用低等优点成为脑功能成像技术研究的热点。 各种脑成像技术优劣对比见表1

表1 各种脑成像技术优劣对比 Table 1 Comparison of various brain imaging parameters

功能性近红外光谱技术(functional near infrared spectroscopy, fNIRS)探测脑组织活动的方式是间接方式。 大脑的各种神经活动需要氧气, 脑部血流(Celebral blood flow)做为大脑供氧的重要渠道, 它对大脑的神经活动是十分敏感的, fNIRS正是利用探测脑血流动力学的参数变化, 来判断脑组织的各种神经活动的。 显而易见的是, 大脑作为人的神经中枢, 各种生理活动如心跳、 呼吸、 脉搏等都会引起神经活动, 也自然会引起脑血流动力学的参数变化, 因此, 如何从包含各种生理干扰的信号中提取脑组织神经活动信号是近红外光谱提取脑功能信号难点。

目前, 常见的提取近红外光脑信号的算法有: 脉搏色素谱法[1]、 EEMD-ICA法[3]主成分分析法(PCA)[5]、 独立成分分析法(ICA)[6]、 相干平均法[7]、 自适应滤波[8]等。 脉搏色素谱法采用注射显影剂吲哚青绿的方式进行示踪估计; EEMD-ICA对信号进行希尔伯特变化, 再进行滤波; PCA和ICA这两种方法适用于连续光谱的近红外光谱系统; 自适应滤波法可以消除心跳干扰; 相干平均法可以滤除脉搏等生理干扰。 上述算法在对近红外脑神经活动信号提取时都取得了一定的进展。

但是, 上述算法在提取近红外光谱脑功能信号时, 重视各种生理干扰, 忽视了在提取信号中的测量干扰, 如仪器精密度、 信号传输中的串扰等。 研究[10, 11]表明, 近红外光谱脑信号用非线性算法提取可以去除测量干扰。 如Fritson使用非线性算法将信号Volterra[11]级数化, 大大提高了信号的信噪比。 本工作提出的扩展的卡尔曼(Extend Klaman Filter, EKF)算法是一种非线性算法。 根据近红外光谱提取的脑信号中的干扰成分是多维非平稳随机的, 利用其时变性、 功率谱不固定, 建立了数学模型, 将信号进行泰勒级数分解, 取一阶函数, 再进行滤波和估计, 进行了高斯(Gauss)滤波实验、 Valsava实验、 视觉诱发实验, 并将EKF算法和主流的EEMD算法进行比对。

1 理论基础
1.1 修正的朗伯比尔定律

朗伯比尔定律(Lambert-Beer law)是分光光度法的基本定律, 是描述物质对某一波长光吸收的强弱与吸光物质的浓度及其厚度间的关系。 适用于所有的电磁辐射和所有的吸光物质, 包括气体、 固体、 液体、 分子、 原子和离子。 其公式为

A=lnI0I=ελcd(1)

式(1)中, ε λ 为吸光物质在波长λ (nm)下的消光系数, c为吸光物质的浓度(mol· L-1), d为光在物质中的传输距离, I0I分别代表出射、 入射光强, A为吸光度。

但事实上, 光在生物组织中会有散射、 吸收等现象, 图1给出了光在组织传播的路径。 因此, 光在组织中的实际传输距离远大于物质厚度。 为了消除这一误差, 1998年, Deply D T等对朗伯比尔定律进行了修正[12], 修正的朗伯比尔定律见式(2)

A=lnI0I=1nεiλciβλl+G(2)

图1 光在组织中的路径Fig.1 The path of light through an tissue

式(2)中, β λ 是生物组织在波长为λ 的路径差分因子DPF(differential pathlength factor), 具体数值可以由Monte-Carlo仿真[13]得到, εiλ(i=1, 2)依次代表HbO2和HbR在波长λ 时的消光系数, ci(i=1, 2)分别代表HbO2和HbR在生物组织中的浓度, L是光源和光源接收器的距离, DPF和L的乘积可以近似表示光在组织中走过的实际路径。 A为光密度, G为光衰减强度, G可以用算法消除。

1.2 功能性近红外光谱

近红外光谱是指波长在650~950 nm的不可见光, 科学家Jobsis[2]在1977年发现近红外光对人体组织具有良好的透射性, 首次验证了用近红外光谱技术无创监测组织血氧的可行性。 各种生理活动都离不开氧气, 脑组织所需氧气的输送靠的是脑血流中的两种血红蛋白, 氧合血红蛋白(Hb02)和还原血红蛋白(HbR), 当脑血流流经脑组织时, 血液中HbO2释放氧气, 供脑组织使用, HbO2转化为HbR; 当血液流经肺泡时, HbR又获得氧气转化为HbO2。 两种血红蛋白的总量基本恒定, 由上述可知, 人脑组织的生理活动会引起HbR和HbO2两种血红蛋白含量[6]的变化。

人体大脑组织中, 水的含量占据了大部分的比重, 约为75%, 除此之外, 还有氧合血红蛋白(HbO2)、 还原血红蛋白(HbR)和黑色素等。 在650~950 nm近红外光波段, 水的吸收率很低, 氧合血红蛋白(HbO2)和还原血红蛋白(HbR)的吸收率较高, 不同物质对不同波长的光的吸收度如图2, 这一区间也称为光窗期。

图2 组织的吸光度图谱Fig.2 Absorbance spectra of tissue

近红外光谱脑信号提取技术正是基于待测组织中的水、 氧合血红蛋白(HbO2)、 还原血红蛋白(HbR)在近红外波段具有不同吸收光谱特性, 以及近红外光可以穿过0.5~3 cm的外层生物组织到达脑组织这一特性进行测定的。 若分别选取波长为λ 1λ 2的近红外光源, 利用修正的朗伯比尔定律,

可得入射光和出射光的光密度变化, 见式(3)— 式(6)

ΔAλ1=DPF(λ1)L[εHbO2λ1Δc(HbO2)+εHbRλ1Δc(HbR)](3)

ΔAλ2=DPF(λ2)L[εHbO2λ2Δc(HbO2)+εHbRλ2Δc(HbR)](4)

进而计算, 可得。

Δc(HbO2)=εHbRλ2ΔA(λ1)DPF(λ1)-εHbRλ1ΔA(λ2)DPF(λ2)L(εHbO2λ1εHbRλ2-εHbO2λ2εHbRλ1)(5)

Δc(HbR)=εHbO2λ1ΔA(λ2)DPF(λ2)-εHbO2λ2ΔA(λ1)DPF(λ1)L(εHbO2λ1εHbRλ2-εHbO2λ2εHbRλ1)(6)

式中, Δ A是光密度变化, ε 是某一物质在某一波长的消光系数, L是光源和光源接收器的距离, DPF和L的乘积可以近似表示光在组织中走过的实际路径, Δ c代表的是浓度变化。 进而计算可得氧合血红蛋白(HbO2)和还原血红蛋白(HbR)两种物质的浓度变化。

2 扩展的卡尔曼滤波算法及近红外光谱模型的建立
2.1 卡尔曼滤波算法

卡尔曼滤波(Kalman filtering)是1961年由Rudolf Kalman提出, 是一种利用线性系统状态方程, 通过系统输入输出观测数据, 基于误差平方和最小的原理进行递归计算的方法, 其本质是参数化的beiyesi模型, 通过对下一时刻系统的初步状态估计以及测量得出的反馈相结合, 得到该时刻无限逼近真实值的状态估计。 在现代随机最优控制和随机信号处理技术中, 信号和噪声往往是多维非平稳随机过程, 对于离散域线性系统, 其线性系统状态预测方程式

xk=Axk-1+Bxk-1+wk-1(7)

p(w)~N(0, Q)(8)

cov(w)=E(wwT)=Q(9)

式中, xkxk-1分别为在k时刻、 k-1时刻的状态值, uk-1为在k-1时刻的控制输入, wk-1为过程噪声, 并且服从高斯分布; A为状态转移系数矩阵, B为控制输入的增益矩阵, Q为过程激励噪声的协方差矩阵。

测量方程为

Zk=Hxk+vk(10)

P(v)~(0, R)(11)

cov(v)=E(vvT)=R(12)

式中, Zkk时刻的测量值, Vk为测量噪声, H为测量矩阵, R为测量噪声协方差矩阵。 并且满足: E(w)=E(v)=0。

2.2 扩展的卡尔曼滤波及模型建立

扩展的卡尔曼滤波是在卡尔曼滤波的基础上, 将非线性系统函数作泰勒(Taylor)级数展开, 再进行线性的卡尔曼滤波, 这就是扩展的卡尔曼滤波(extend Kalman filter, EKF)。

考虑离散时间非线性动态系统

xk+1=f(xk, wk)zk+1=h(xk+1, vk+1)(13)

其中, w(k)是过程噪声, 协方差是Qk ; Vk是测量噪声, 其方差是Rk+1f是状态变量x和时间k相关的非线性函数, h是状态方程中与状态变量x和观测量z相关的非线性函数。 进行泰勒(Taylor)展开, 取展开后的一阶线性部分作为非线性模型的近似。

状态方程表示为

xk=f(xk-1|k-1)+fxxk-1|k-1+wk-1=Fk-1xk-1+wk-1+uk-1(14)

式(14)中, f(xk|k)进行泰勒展开,

Fk-1=fxxk-1|k-1(15)

Fk=f1[xk]x1kf1[xk]x2kf1[xk]xnkf2[xk]x1kf2[xk]x2kf2[xk]xnkfn[xk]x1kfn[xk]x2kfn[xk]xnk(16)

测量方程表示为

zk=h(xk|k-1)+hxxk|k-1(xk-xk|k-1)+Vk=Hkxk+Vk(17)

式(17)中, Gh(xk)表示在xk|kG处进行泰勒展开,

Hk=hxxk|k-1Hk=h1[xk]x1kh1[xk]x2kh1[xk]xnkh2[xk]x1kh2[xk]x2kh2[xk]xnkhn[xk]x1khn[xk]x2khn[xk]xnk(18)

参考经典卡尔曼滤波的推导过程, 可以得到扩展卡尔曼滤波的五个方程式(19)— 式(23)

xk|k-1=f(xk-1|k-1)+wk-1(19)

pk|k-1=Fpk-1|k-1FT+Q(20)

Kk=pk|k-1HT(Hpk|k-1HT+R)-1(21)

xk|k=xk|k-1+Kk(Zk-h(xk|k-1))(22)

pk|k=(1-kkH)pk|k-1(23)

式(19)是状态预测方程, 式(20)是误差协方差方程, 式(21)是卡尔曼增益方程, 式(22)是滤波校正方程, 式(23)是误差协方差更新方程。

2.3 EKF模型的建立

简而言之, 人脑的神经活动需要氧气, 氧气需要脑血流中的两种血红蛋白进行输送, 必然会引起氧合血红蛋白(HbO2)和还原血红蛋白(HbR)浓度的变化。 由扩展的朗伯比尔定律及近红外光谱的特性可知, 这种变化可以由近红外光谱测量得到, 此时测量的数据是包含了生理噪声和测量噪声, 将测量所得到的数据经扩展的卡尔曼滤波, 取得HbO2和HbR浓度的变化。 图3表述了上述关系。

图3 脑血流动力学数学模型图Fig.3 Overview diagram of the applied brain hemodynamic model

3 实验结果与讨论
3.1 扩展的卡尔曼滤波(EKF)对Gauss噪声滤波有效性的仿真实验

近红外光谱脑信号采集过程中的噪声干扰符合高斯分布, 取干扰信号w(k)和测量噪声信号v(k)的幅值均为0.01的高斯白噪声信号, 输入信号幅值为1.0、 频率为1.5的正弦信号。 使用EKF实现信号的滤波, 取Q=1, R=1。 仿真时间是5 s, 在MATLAB中编程进行仿真, 图4是在信号中叠加了高斯(Gauss)噪声的波形图, 图5是经EKF滤波后的信号。 从MATLAB仿真结果看, 信号平滑、 无失真、 没有畸变, 该算法在去除高斯(Gauss)白噪声方面效果明显。

图4 叠加噪声的原始信号Fig.4 Ideal signal and signal with noise

图5 EKF滤波的信号Fig.5 Ideal signal and EKF filtered signal

3.2 EKF提取HbR和HbO2浓度变化

3.2.1 实验装置

为了测量整个前额叶区域的血红蛋白和还原血红蛋白的变化量, 设计了近红外光谱脑信号采集装置。 其中, 直流稳压电源和光源驱动电路给近红外光光源稳定供电, 光源为波长750和830 nm的二极管近红外光源, 用同相放大和低通滤波做初步处理, 将处理后的信号发送至上位机, 上位机用C#编写, 其框图见图6。

图6 装置设计框图Fig.6 Experimental design block diagram

3.2.2 Valsava实验与分析

为了获得HbO2和HbR在血液中的浓度变化曲线, 使用波长分别为750和830 nm的近红外光源进行Valsava实验, 近红外光源和光源探测器用黑布完全遮起来, 被试者是25岁成年男性, 身体健康, 无烟酒嗜好。 耳朵中放置隔音棉, 尽可能减少外界环境干扰, 进行闭气-呼气实验。 在实验开始的0~60 s, 被试者正常呼吸; 第60~120 s, 被试者闭气; 第120 s之后, 被试者恢复正常呼吸。 采集到的HbO2和HbR经扩展的卡尔曼滤波后变化曲线如图7所示。

图7 Valsava实验血红蛋白变化曲线图Fig.7 The hemoglobin waveform extracted by Valsava experiment

3.2.3 视觉诱发实验与分析

为了验证本方法在提取氧合血红蛋白(HbO2)和还原血红蛋白(HbR)浓度变化的对外界刺激的敏感性, 进行视觉诱发实验。 Arturs发现, 视觉诱发实验会引起脑血流参数的变化[14], Zaletel发现视觉诱发实验和脑血流流速变化是正相关关系[15]。 本工作采用黑白相间的图片作为视觉诱发源, 如图8所示。

图8 视觉诱发源Fig.8 Visual evoked

将视觉诱发源信号在电脑屏幕上交替播放, 交替频率是5 Hz, 近红外光源和光源探测器用黑布完全遮起来, 被试者是25岁成年男性, 身体健康, 无烟酒嗜好。 耳朵中放置隔音棉, 尽可能减少外界环境干扰。 让被试者观看图片, 测得氧合血红蛋白(HbO2)和还原血红蛋白(HbR)浓度变化如图9所示。

图9 视觉诱发实验血红蛋白浓度变化曲线图Fig.9 Hemoglobin concertration extracted by visual evoked

3.3 EKF和EEMD算法比对

EKF算法在提取近红外光脑信号的有效性得到了验证, 进一步对其滤波效果进行评价, 比对方法是利用均方根误差RMSE (root mean square error)、 相关系数r(correlation coefficient)。 RMSE可以表示测量值和真值的离散情况, 参数越小, 表明滤波效果越好。 r表示的是相关系数, r越接近1, 表示效果越好。 对同一被试者采集10组数据, 分别计算EEMD算法和EKF算法的RMSE值、 r值, 结果如表2。 可以看出, EKF算法对比EEMD算法, 其RMSE值提高了0.96%, r值提高了0.6%。

表2 EKF和EEMD算法比对 Table 2 Comparison of EKF and EEMD algorithm
4 结论

设计了功能性近红外光谱(fNIRS)脑血流信号采集装置, 进行了Valsava和视觉诱发实验, 采用扩展的卡尔曼滤波(EKF)建立数学模型, 对所采集信号进行了EKF处理。 通过实验表明: 本文所提出的装置和算法可以有效去除符合高斯分布的干扰, 提取出脑血流中氧合血红蛋白(HbO2)和还原血红蛋白(HbR)变化曲线图; 和主流的EEMD提取脑信号算法比对其RMSE值提高了0.96%, r值提高了0.6%。 表明本文所提出的方法有一定的优越性。 并可以为相关脑部疾病诊断如癫痫病灶定位、 抑郁研究等提供了一种有效的脑信号提取途径。

参考文献
[1] Yi M H, Qiu T, Okubo M, et al. Vibrational Spectroscopy, 2018, 95: 23. [本文引用:1]
[2] Jobsis F F. Science, 1977, 198(4323): 1264. [本文引用:1]
[3] ZHA Yu-tong, LIU Guang-da(查雨彤, 刘光达). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2015, 35(10): 2746. [本文引用:1]
[4] Gosselin R, Rodrigue D, Duchesne C. Computers & Chemical Engineering, 2011, 35(2): 296. [本文引用:1]
[5] Liu G D, Wei X, Zhang S, et al. Optik, 2018, 160: 168. [本文引用:1]
[6] Zhan Y H, Brooks D H, Franceschini M A, et al. Journal of Biomedical Optics, 2005, 10(1): 011014. [本文引用:2]
[7] Virtanen J, Noponen T. Journal of Biomedical Optics, 2009, 14: 054032. [本文引用:1]
[8] Gratton E, Toronov V, Wolf U, et al. Journal of Biomedical Optics, 2005, 10(1): 011008. [本文引用:1]
[9] Brunner S, Fomin P, Kargel Ch. Waste Management, 2015, 38: 49. [本文引用:1]
[10] Obata T, Liu T T, Miller K L, et al. Neuroimage, 2003, 21: 144. [本文引用:1]
[11] Fritson K J, Josephs O, Rees G. Magnetic Resonance in Medicine 1998, 39: 41. [本文引用:2]
[12] Delpy D T, Cope M, Zee P V D, et al. Physics in Medicine & Biology, 1988, 33(12): 1433. [本文引用:1]
[13] Patterson M S, Wilson B C, Wyman D R. Lasers in Medical Science, 1991, 6(2): 155. [本文引用:1]
[14] Zhang Tao, Shao Changyan, Wang Xinnian. Atmospheric Scattering-Based Multiple Images Fog Removal, 4th International Congress on Image & Signal Processing, 2011, 10.1109/CISP.2011.6100007. [本文引用:1]
[15] LIU Guo, Qun-bo, LIU Yang-yang(刘国, 吕群波, 刘扬阳). Acta Photonica Sinica(光子学报), 2018, 47(2): 210001. [本文引用:1]