基于数学形态学和中值滤波的fNIRS信号运动校正算法研究

赵杰; 乔吉日木图; 丁雪桐; 梁晓敏

doi:doi:10.3788/AOS202040.2230002

光学学报, 2020, 40 (22): 2230002, 网络出版: 2020-10-25

基于数学形态学和中值滤波的fNIRS信号运动校正算法研究下载： 1150次

fNIRS Signal Motion Correction Algorithm Based on Mathematical Morphology and Median Filter

论文大纲

赵杰乔吉日木图 ^*丁雪桐梁晓敏

作者单位

河北大学电子信息工程学院, 河北保定 071002

光谱学功能性近红外光谱中值滤波数学形态学运动伪迹 spectroscopy functional near-infrared spectroscopy median filtering mathematical morphology motion artifacts

AI 词云图 AI语音精读 AI语音超短摘要

注：本部分内容由 AI 自动生成，请您知悉。

摘要

功能性近红外光谱技术(fNIRS)作为一种新兴的神经成像技术得到了广泛关注,然而fNIRS信号中运动伪迹的存在会使信号的处理结果产生偏差。提出了一种定向中值滤波和数学形态学相结合的算法——tMedMor算法,并采用该算法对fNIRS信号中的三种运动伪迹(包括尖峰、基线突变和缓慢漂移)进行去除;然后用仿真数据和实验数据进行了验证,并将所提算法与常用的几种算法进行对比,结果表明:tMedMor算法在均方误差、信噪比、皮尔逊相关系数的平方、峰峰误差方面具有良好的表现,说明该算法可以作为一种新方法用于fNIRS信号的预处理阶段。

Abstract

Functional near-infrared spectroscopy (fNIRS) has attracted widespread attention as an emerging neuroimaging technology. However, the existence of motion artifacts in the fNIRS signal leads to bias in its signal processing outcomes. We proposed a tMedMor algorithm that combines the targeted median filtering (tMed) and mathematical morphology (Mor) for the removal of three motion artifacts in the fNIRS signal, namely, spike, baseline shift, and slow drift. Simulated and experimental data were used for verification, and the performance of the proposed algorithm was compared with those of several other common algorithms. Our results revealed that the tMedMor algorithm demonstrates good performance in terms of mean square error, signal-to-noise ratio, square of Pearson correlation coefficient, and peak-to-peak error, which together indicate that tMedMor can be applied as a new approach to the fNIRS signal at the preprocessing stage.

1 引言

功能性近红外光谱技术(fNIRS)是近几年受到广泛关注的一种神经成像技术,其原理是用光源向脑组织发射近红外光,然后用探测器测量来自脑组织的背向散射光的强度,再用修正的比尔-朗伯定律计算血液中有氧血红蛋白(HbO)和脱氧血红蛋白(HbR)的浓度,以反映被试在执行不同任务时的大脑激活状态^[1-2]。与其他成像技术相比,fNIRS具有价格低廉、易携带、无损等优点^[1,3-4],已被广泛应用于步行、舞蹈、音乐表演、言语交流等涉及较多运动的任务范式中^[5-6]。在采用fNIRS采集信号时,被试头部的运动会导致光源和头皮发生分离,致使接收信号的幅度产生陡峭变化,从而产生不同类型的运动伪迹,如尖峰、基线突变和缓慢漂移等^[7-8]。这些运动伪迹会使数据处理结果产生较大偏差^[9],鉴于此,研究人员提出了很多处理运动伪迹的方法。

早期处理运动伪迹的方法是直接丢掉含有运动伪迹的通道或试次,但是当实验试次较少时,这种方法并不适用^[10]。另外还有一种是基于硬件的方法,该方法用短距离通道测量与生理信号无关的纯运动伪迹信号^[11-13],或者用加速度计记录被试在数据采集过程中的运动数据^[14-15],然后用自适应滤波方法或回归方法将运动伪迹从信号中去除。但该类方法对实验设备的要求较高,因此,很多不需要额外测量辅助信号的后处理技术被提出,较常用的有以下几种:运动伪迹去除(MARA)算法、目标主成分分析(tPCA)算法、峰度小波(kWavelet)算法、时域微分分布校正(TDDR)算法。基于小波的方法^[15-16]能够很好地去除尖峰,但是会加剧信号中的基线突变。MARA和tPCA算法^[17-18]依赖于对运动伪迹的检测,需要使用者指定几个参数,而这些参数会随着实验仪器、被试以及通道的变化而改变,给使用者增加了很多困难。TDDR算法^[9]去除尖峰和基线突变的效果较好,但不能很好地去除缓慢漂移。将两种方法组合起来进行研究成为近年来的热点。与单一方法相比,组合方法可以利用两种方法各自的优点,去噪效果更好。例如:将样条插值与SG(Savitzky-Golay)滤波结合的Spline-SG算法使用样条插值去除基线突变,使用SG滤波去除尖峰,去伪迹效果大大提升^[19]。

数学形态学方法具有原理简单、实现效率高^[20]等优点,它是以随机集论为基础建立的。随机集论适合描述信号的随机性质,它通过选择合适的结构元素和形态运算来保留信号的形态信息^[21],更适合处理运动伪迹。而普通的FIR(finite impulse response)或IIR(infinite impulse response)滤波器虽然可以去除生理噪声,但由于运动伪迹的频率范围较宽,甚至与有用信号重叠,所以采用FIR或IIR滤波器去运动伪迹会使信号失真^[22]。中值滤波算法常被用于信号中脉冲噪声的去除,该方法能够很好地保留信号的细节^[23-24]。鉴于数学形态学方法和中值滤波算法各自的优点,本文对中值滤波算法进行改进,提出了定向中值滤波(tMed)算法,然后将其与数学形态学(Mor)算法结合起来用于去除fNIRS信号中的运动伪迹。本文首先通过计算信号的近似梯度和滑动标准差来检测fNIRS信号中的运动伪迹,然后用改进的中值滤波算法去除尖峰,最后用Mor算法去除基线突变和缓慢漂移。本文将所提算法与现有的几种算法进行比较,验证了所提算法——定向中值滤波与数学形态学结合(tMedMor)算法在去除fNIRS信号运动伪迹方面的有效性。

2 理论与算法

图1为用tMedMor算法对fNIRS信号进行预处理的流程。首先用tMed算法去除尖峰,再用Mor方法去除基线突变和缓慢漂移。

图 1. 采用tMedMor算法对fNIRS信号进行预处理的流程图

Fig. 1. fNIRS signal preprocessing flowchart using tMedMor algorithm

下载图片查看所有图片

2.1 tMed算法

Mor算法在去除两个相邻较近的尖峰时会引入基线突变,为克服这一缺点,本文提出了tMed算法。在使用Mor方法进行处理之前,先用tMed算法将信号中的尖峰基本去除。tMed算法首先通过计算信号的近似梯度和滑动标准差来检测运动伪迹^[19],然后用改进的中值滤波算法处理运动伪迹。该算法的具体步骤如下:

1) 运动伪迹的检测

先用截止频率为2 Hz的低通滤波器对光密度(OD,optical density)信号I进行滤波得到I_lpf,然后根据采样定理将信号重采样到4 Hz得到I_rs,最后用索贝尔滤波器计算OD信号的近似梯度G。计算公式为

G = [- 1 0 1] * I_{rs}, (1)

式中:*表示卷积。由于对信号进行了重采样,所以(1)式的计算结果与采样率无关,避免了采样率对运动伪迹检测结果的影响。

因为运动伪迹会在梯度信号G中产生异常值,所以可以通过找出G中的异常值来检测运动伪迹。用Q₁、Q₂、Q₃分别表示第一、第二、第三分位点。Q₁、Q₃分别是数据集G从小到大排列后的前一半数据和后一半数据的中值,Q₂为整个数据集G的中值。四分间距(R_IQ)等于Q₃减去Q₁,将梯度值中小于Q₁-1.5×R_IQ或者大于Q₃+1.5×R_IQ的采样点作为异常值点。同理,用5 s的窗口长度来计算OD信号的滑动标准差,然后用同样的过程找出滑动标准差中的异常值点。所有这些异常值点被标记为运动伪迹。

2) 用改进的中值滤波算法去除尖峰

采用中值滤波算法对OD时间序列中有运动伪迹的采样点进行处理。因为运动伪迹中的尖峰幅值较大,普通中值滤波算法的校正效果较差,所以本文提出了一种改进的中值滤波算法。该算法首先以某一个含运动伪迹的采样点为中心,向左右两边逐渐扩大窗口宽度,直至选出W个被标记为无运动伪迹的采样点,然后用这些采样点的中值作为中心采样点的信号值,这样就避免了含运动伪迹的信号重新被计入而造成的不必要的误差。

2.2 Mor算法

数学形态学用腐蚀、膨胀、开运算、闭运算及其组合,以及不同形状和长度的结构元素对信号进行处理。s(x)关于g(x)的腐蚀运算定义为

\begin{array}{l} (s Θ g) (x) = \min_{y = 0,1, \dots, M - 1} [s (x + y) - g (y)], \\ x = 0,1, \dots, N - M 。 (2) \end{array}

s(x)关于g(x)的膨胀运算定义为

\begin{array}{l} (s \oplus g) (x) = \max_{y = 0,1, \dots, M - 1} [s (x - y) + g (y)], \\ x = 0,1, \dots, N + M - 2, (3) \end{array}

式中:s(x)为定义在D₁={0,1,…,N-1}上的一维离散信号;N表示信号的采样点总数;结构元素g(x)是一个定义在D₂={0,1,…,M-1}上的一维离散函数;M为结构元素的长度,N≥M。s(x)关于g(x)的开运算为

(s ⚪ g) (x) = (s Θ g \oplus g) (x), (4)

s(x)关于g(x)的闭运算为

(s • g) (x) = (s \oplus g Θ g) (x) 。 (5)

开闭和闭开运算组合并取平均可构成形态学低通滤波器OC_CO,即

s_{o} = f_{OC_CO} (s_{i}, g) = (s_{i} ⚪ g • g + s_{i} • g ⚪ g) / 2, (6)

式中:s_i为滤波器的输入信号;s_o为滤波器的输出信号;g为滤波器的结构元素。

图2为Mor算法流程图(t为时间)。图2(a)为形态学低通滤波器,用来去除信号中剩余的尖峰;图2(b)为形态学高通滤波器,用来去除信号中的基线突变和缓慢漂移。

图 2. Mor算法流程图

Fig. 2. Flow chart of Mor algorithm

下载图片查看所有图片

该算法的具体步骤如下:

1) 输入信号s₁(t)经过形态学低通滤波处理(去除尖峰)后得到信号s₂(t),即

s_{2} (t) = f_{OC_CO} (s_{1}, g_{1}) 。 (7)

2) 提取基线突变和缓慢漂移。对s₂(t)进行形态学低通滤波处理得到基线突变和缓慢漂移s₃(t),即

s_{3} (t) = f_{OC_CO} (s_{2}, g_{2}) 。 (8)

3) 从s₂(t)中减去s₃(t)得到去运动伪迹后的信号s₄(t),即

s_{4} (t) = s_{2} (t) - s_{3} (t), (9)

式中:g₁、g₂均为直线型结构元素。

下面以形态学低通滤波器去除基线突变和缓慢漂移为例说明以秒(s)为单位的结构元素长度L₁、L₂的选取。

用形态学低通滤波器提取基线突变和缓慢漂移时,信号的衰减幅度由比值K^[25]决定,即

K = \frac{M_{1}}{M_{2}} = \frac{L_{2} f_{s}}{f_{s} \times \frac{1}{f_{L}}} = L_{2} f_{L}, (10)

式中:M₁是以采样点为单位的形态学滤波器结构元素的长度;M₂是以采样点为单位的OD信号的低频周期;f_s为信号的采样频率;f_L为信号的低频主峰频率。当K≥0.5时,信号完全被滤除掉;K越小,基线突变以及缓慢漂移才能越多地被保存下来。因此,K=0.5时,既能完全检测到基线突变和缓慢漂移并将其去除,又能保留信号的低频部分,同时高频部分也更能较好地被保留,此时

L_{2} = \frac{K}{f_{L}} = \frac{0.5}{f_{L}} 。 (11)

同理,可得形态学低通滤波器的结构元素长度L₁为

L_{1} = \frac{0.5}{f_{H}}, (12)

式中:f_H为信号的高频主峰频率。选取参数时,首先对一组无噪声OD信号进行傅里叶变换,然后找到低频主峰频率f_L和高频主峰频率f_H,然后由(11)、(12)式计算得到L₁和L₂。

3 实验部分

3.1 仿真验证

3.1.1 仿真数据

仿真的原始强度信号由几种成分相加而成,即

Φ_{s}^{(λ)} (t) = Φ_{b}^{(λ)} (t) + Φ_{c}^{(λ)} (t) + Φ_{m}^{(λ)} (t) + Φ_{n}^{(λ)} (t) + Φ_{o}^{(λ)} (t), (13)

式中: $Φ_{s}^{(λ)}$ (t)是波长为λ时的强度信号; $Φ_{b}^{(λ)}$ (t)是由气球模型^[26]仿真得到的与血液动力学响应相关的成分; $Φ_{c}^{(λ)}$ (t)是与脉搏相关的成分,可由动态心率信号模拟得到; $Φ_{m}^{(λ)}$ (t)是与血压相关的成分,可由频率随时间变化的低频正弦振荡模拟得到; $Φ_{n}^{(λ)}$ (t)是宽带噪声成分,可由正态分布模拟得到; $Φ_{o}^{(λ)}$ (t)是信号的补偿,为一常数,具体的模型和使用参数见文献[ 27]。仿真数据的通道数为1,波长分别为690 nm和830 nm,采样频率为25 Hz。仿真数据由20个组块的任务组成,每个组块的持续时间为10 s,每两个组块之间的休息时间为10 s,信号总时长约为7 min。

将原始强度信号转化成OD信号,并将三种类型的运动伪迹加到其中。基线突变由阶跃函数生成,幅值在0.3~0.5内随机选取,出现的频率为每分钟1次,出现位置由均匀分布的随机数决定。将生成的基线突变加到仿真信号的前半段。采用指数增长曲线z(t)=b^t/τ模拟尖峰^[9],其中增长因子b从标准差为25的正态分布数据中随机选择,时间常数τ在0~0.5之间随机选择,尖峰出现的频率为每分钟2次。将生成的尖峰加到整个仿真信号中。缓慢漂移采用二次函数模拟,幅值在0.4~0.8之间随机生成,持续时间在5~7 s之间随机生成,每80 s随机生成一段缓慢漂移,4 min共生成3个缓慢漂移。将生成的缓慢漂移加到仿真信号的后半段。采用以上方法生成20组仿真数据。

3.1.2 数据处理流程

分别采用本文所述的tMedMor、tMed、Mor算法和MARA、tPCA、kWavelet、Spline-SG、TDDR 算法对仿真OD信号进行校正,然后使用MATLAB Homer2工具箱中的hmrBandpassFilt函数对信号进行带通滤波处理,函数参数为hpf=0.01和lpf=0.07,再用修正的比尔-朗伯定律将OD信号转化成血红蛋白浓度信号,最后对HbO信号进行块平均,范围为-5~20 s。

3.1.3 tMedMor方法的参数选择

进行中值滤波时的窗口宽度为W=5 s。无噪OD信号的单边功率谱密度如图3所示(d_PS表示功率谱密度)。由图3可得到高频主峰频率f_H=0.25 Hz,从而由(12)式得到形态学低通滤波器的结构元素长度L₁=2 s;由图3还可得到低频主峰频率f_L=0.05 Hz,从而由(11)式得到形态学高通滤波器的结构元素长度L₂=10 s。

图 3. 无噪OD信号的单边功率谱密度

Fig. 3. Unilateral power spectral density of noiseless OD signal

下载图片查看所有图片

3.2 实验验证

3.2.1 实验数据

在真实静息态数据中加入仿真的血红蛋白浓度信号构成实验数据。静息态数据采用NITRC网站中的DATASET2数据集(网址为https:∥www.nitrc.org/projects/fnirsdata/),该数据集是使用TechEn CW6系统采集的5名正常被试在静息态下的fNIRS信号。光纤探头由15个光源、18个长距离探测器和14个短距离探测器组成,长距离探测器和短距离探测器到光源的距离分别为30 mm和8 mm,采样频率为50 Hz。仿真血液动力学响应函数(HRF)由Gamma函数模拟生成,其峰值位置在6 s左右,持续时长为16 s。采用组块设计实验任务,每个组块持续时间为10 s,每两个组块间隔5~10 s,在6.5 min内共生成17~19个组块。将与任务相关的矩形函数与HRF卷积后生成仿真血红蛋白浓度信号;然后将血红蛋白浓度信号加到静息态数据上,使690 nm处的信号相对于基线发生1%的变化,830 nm处的信号相对于基线发生2%的变化,导致HbO浓度增加了0.6 μmol/L,HbR浓度减少了0.2 μmol/L。两种波长的路径长度因子均为6。

由于原信号中所含运动伪迹较少,故添加了额外的运动伪迹。首先将原始强度信号转化成OD信号,然后按不同的方式生成三种运动伪迹(即每120 s交替产生一个基线突变或一个缓慢漂移,每120 s生成一个尖峰,运动伪迹模型见3.1.1节),最后将生成的运动伪迹加到OD信号中。

3.2.2 数据处理流程

分别采用tMedMor、tMed、Mor算法以及对比算法(MARA、tPCA、kWavelet、Spline-SG、TDDR算法)对OD信号进行校正,再用截止频率为0.5 Hz的低通滤波器进行滤波,最后用一般线性模型(GLM)估计HRF。GLM采用最小二乘法估计连续时域高斯基函数的权值,高斯基函数的标准差为0.5 s,均值间距为0.5 s。因为数学形态学本身可以去趋势,所以tMedMor、Mor算法没有用多项式拟合去趋势,其他方法均采用三阶多项式拟合去趋势。

3.2.3 tMedMor方法的参数选择

进行中值滤波时的窗口宽度为W=5 s。选择一组无噪OD信号进行傅里叶变换,结果如图4所示。由图4可以得到高频主峰频率f_H=0.25 Hz,从而由(12)式得到形态学低通滤波器的结构元素长度L₁=2 s;由图4还可以得到低频主峰频率f_L=0.04 Hz,从而由(11)式得到形态学高通滤波器的结构元素长度L₂=12.5 s。

图 4. 无噪OD信号的单边功率谱密度

Fig. 4. Unilateral power spectral density of noiseless OD signal

下载图片查看所有图片

3.3 对比方法

对比算法中的MARA、tPCA、kWavelet、TDDR用Homer2工具箱中的函数实现,Spline-SG算法用原论文中提供的函数实现^[19]。MARA算法使用函数hmrMotionArtifactByChannel检测含运动伪迹的信号段,用函数hmrMotionCorrectSpline去除运动伪迹,所用参数如下^[17]:tMotion=0.5,tMask=2,STDEVthresh=20,AMPthresh=0.5,pSpline =0.99。tPCA方法用函数hmrMotion-CorrectPCArecurse实现,检测运动伪迹的参数与MARA算法中对应的参数相同,其他参数如下^[18]:nSV=0.97,maxIter=3。kWavelet算法用hmrMotionCorrectKurtosisWavelet函数实现,参数^[15]kurtosis=3.3。Spline-SG算法用hmrMotion-CorrectSplineSG函数实现,函数参数如下^[19]:p=0.99,iqr=1.5。根据采样频率不同,本文设置如下:仿真数据中的SG_winSize=301,实验数据中的SG_winSize=601。TDDR算法用hmrMotionCorrectTDDR函数实现。

3.4 评价指标

采用块平均HbO浓度信号的均方误差E_MS、信噪比R_SN、皮尔逊相关系数的平方R²、峰峰误差E_p这四个评价指标比较不同算法去运动伪迹的效果。上述指标的计算公式为

\begin{array}{l} E_{MS} = \frac{1}{N} \overset{N}{\sum_{i = 1}} [x (t_{i}) - y (t_{i} {)]}^{2}, (14) \\ R_{SN} = 10 \lg (\frac{\overset{N}{\sum_{i = 1}} x^{2} (t_{i})}{\overset{N}{\sum_{i = 1}} [x (t_{i}) - y (t_{i} {)]}^{2}}), (15) \\ R^{2} = \frac{1}{M} \overset{N}{\sum_{i = 1}} [\frac{x (t_{i}) - < x (t) >}{s_{x}}] [\frac{y (t_{i}) - < y (t) >}{s_{y}}], (16) \\ E_{p} = \frac{(y_{p} - x_{p})^{2}}{x_{p}^{2}} \times 100, (17) \end{array}

式中: $s_{x} = \sqrt[]{\frac{1}{M} \overset{N}{\sum_{i = 1}} [x (t_{i}) - < x {(t) >]}^{2}}; s_{y} = \sqrt[]{\frac{1}{M} \overset{N}{\sum_{i = 1}} [y (t_{i}) - < y {(t) >]}^{2}};$ x(t)和y(t)分别为无噪信号和已去除运动伪迹的信号;x_p、y_p分别表示信号x(t)和y(t)的峰值;N是信号x(t)和y(t)的采样点总数;M=N-1。E_MS和R_SN评估信号x(t)和y(t)的一致性,R²评估两信号之间的相似性。E_MS和E_p越小,去运动伪迹效果越好;R_SN越大,去运动伪迹效果越好;R²越接近1,去运动伪迹效果越好。

4 结果与讨论

4.1 仿真验证

4.1.1 OD信号波形图

采用tMedMor算法处理后的某个通道中波长为830 nm的OD信号的波形图如图5所示,I_OD为OD信号的强度。tMed算法已将大部分尖峰去除,并将部分缓慢漂移转变为基线突变;而Mor算法去除了基线突变、剩余的缓慢漂移和少量幅值较低的尖峰。将去运动伪迹后的OD信号与不含运动伪迹的OD信号进行对比后可以发现,tMedMor算法能较好地去除尖峰、基线突变、缓慢漂移三种运动伪迹,并能将有用信号还原。

采用tMedMor方法处理仿真数据时的波形图,图中数据选自某个通道中波长为830 nm的OD信号。(a)被运动伪迹污染的OD信号;(b)用tMed算法校正后的OD信号;(c)用Mor算法校正后的OD信号;(d)用tMedMor算法校正后的OD信号与无运动伪迹的OD信号的对比

图 5. 采用tMedMor方法处理仿真数据时的波形图,图中数据选自某个通道中波长为830 nm的OD信号。(a)被运动伪迹污染的OD信号;(b)用tMed算法校正后的OD信号;(c)用Mor算法校正后的OD信号;(d)用tMedMor算法校正后的OD信号与无运动伪迹的OD信号的对比

Fig. 5. Waveform of simulation data after using tMedMor algorithm, the data in the figure was selected from OD signal with a wavelength of 830 nm in a certain channel. (a) OD signals contaminated by motion artifacts; (b) OD signal corrected by tMed algorithm; (c) OD signal corrected by Mor algorithm; (d) comparing OD signal corrected by tMedMor algorithm with OD signal without motion artifacts

下载图片查看所有图片

4.1.2 各算法的对比

为评估tMedMor、tMed、Mor以及MARA、tPCA、kWavelet、Spline-SG、TDDR去运动伪迹的性能,计算经这些算法处理后的仿真块平均HbO信号的E_MS、R_SN、R²和E_p(对20组数据的E_MS、R_SN、R²、E_p进行平均),结果如表1和图6所示。

表 1. 经不同算法处理后,仿真块平均HbO信号的E_MS、R_SN、R²和E_p

Table 1. E_MS, R_SN, R² and E_p of simulated block-average HbO signal after processing by different algorithms

Method	${E_{MS}}_{}$ /10^-15	${R_{SN}}_{}$ /dB	R²	E_p
Uncorrected	3561.7±3182.1	-24.57±3.81	0.53±0.42	29136.8±43210.4
MARA	498.5±666.6	-14.34±6.04	0.74±0.25	3933.7±4644.9
TDDR	38.4±24.4	-5.44±2.45	0.75±0.24	265.9±228.4
Spline-SG	125.3±110.1	-9.67±4.97	0.68±0.27	1468.9±2258.0
kWavelet	1704.5±1414.4	-19.97±6.40	0.49±0.29	17613.0±27820.8
tPCA	3611.4±3844.7	-24.24±3.89	0.57±0.36	27700.4±31193.8
Mor	26.5±34.3	-0.50±7.01	0.63±0.37	227.9±381.3
tMed	3618.8±3917.1	-23.95±4.22	0.40±0.36	28697.4±50120.2
tMedMor	6.9±4.9	3.18±5.25	0.74±0.26	60.5±78.0

查看所有表

采用不同算法进行运动伪迹校正后,计算得到的仿真块平均HbO信号的评价指标。(a) EMS; (b) RSN; (c) R2; (d) Ep

图 6. 采用不同算法进行运动伪迹校正后,计算得到的仿真块平均HbO信号的评价指标。(a) E_MS; (b) R_SN; (c) R²; (d) E_p

Fig. 6. Evaluation indices calculated after correcting motion artifacts of simulated block-average HbO signal with different algorithms. (a) E_MS; (b) R_SN; (c) R²; (d) E_p

下载图片查看所有图片

由表1和图6看出:与校正之前的结果相比,采用tMedMor、Mor、TDDR、Spline-SG算法处理后,块平均HbO信号的E_MS、R_SN、E_p均有很大改善,而且经tMedMor处理后的块平均HbO信号的E_MS、R_SN、E_p最优;经tMedMor、MARA、TDDR处理后的块平均HbO信号的R²与校正之前相比有了一定改善,其中经TDDR处理后的块平均HbO信号的R²最大;经tPCA和tMed处理后的块平均HbO信号的效果较差;与单独用tMed或Mor处理相比,采用两者结合后的tMedMor算法处理后的块平均HbO信号的性能更好。

4.2 实验验证

4.2.1 OD信号波形图

图7所示为采用tMedMor算法处理真实数据的示例。由图7可知:tMed算法可将大部分尖峰去除,并可将较宽的尖峰变窄,使之更适合采用Mor方法进行处理,但该算法没有去除基线突变和缓慢漂移;Mor算法可将基线突变和缓慢漂移去除,并可将剩余的少量尖峰去除。将已去除运动伪迹的信号与无运动伪迹的信号进行对比后发现,tMedMor算法基本去除了三种运动伪迹,成功地还原了任务态信号。

采用tMedMor算法处理真实数据时的波形图,图中数据选自某个通道中波长为830 nm的OD信号。(a)被运动伪迹污染的OD信号;(b)采用tMed算法校正后的OD信号;(c)采用Mor算法校正后的OD信号;(d) 用tMedMor算法校正后的OD信号与无运动伪迹的OD信号的对比

图 7. 采用tMedMor算法处理真实数据时的波形图,图中数据选自某个通道中波长为830 nm的OD信号。(a)被运动伪迹污染的OD信号;(b)采用tMed算法校正后的OD信号;(c)采用Mor算法校正后的OD信号;(d) 用tMedMor算法校正后的OD信号与无运动伪迹的OD信号的对比

Fig. 7. Waveform of real data after using tMedMor algorithm, the data in the figure was selected from OD signal with a wavelength of 830 nm in a certain channel. (a) OD signals contaminated by motion artifacts; (b) OD signal corrected by tMed method; (c) OD signal corrected by Mor algorithm; (d) comparing OD signal corrected by tMedMor algorithm with OD signal without motion artifacts

下载图片查看所有图片

4.2.2 各算法的对比

用tMedMor、tMed、Mor以及MARA、tPCA、kWavelet、Spline-SG、TDDR去运动伪迹后,计算了块平均HbO信号的E_MS、R_SN、R²和E_p,并对全部被试的所有通道进行了平均,结果如表2和图8所示。

由表2和图8可知:与去运动伪迹之前的结果相比,采用tMedMor、Mor、TDDR、Spline-SG算法对真实数据进行处理后,块平均HbO信号的E_MS、R_SN、E_p均有很大改善,其中经tMedMor算法处理后的块平均HbO信号最优;采用tMedMor、Mor、Spline-SG、TDDR算法处理后,块平均HbO信号的R²相比去运动伪迹前也有很大提升,其中经TDDR算法处理后的块平均HbO信号的R²最大;经MARA、tPCA、kWavelet算法处理后的块平均HbO信号的E_MS、R_SN、E_p、R²也有一定改善,经kWavelet算法处理后的块平均HbO信号的R²较差;相比于其他算法,tMed算法的表现最差;与单独用tMed算法或Mor算法处理的结果相比,采用两者结合的tMedMor算法处理后的块平均HbO信号的性能更好。

表 2. 经不同算法处理后,真实块平均HbO信号的E_MS、R_SN、R²和E_p

Table 2. E_MS, R_SN, R² and E_p of real block-average HbO signal after processing by different algorithms

Method	${E_{MS}}_{}$ /10^-14	${R_{SN}}_{}$ /dB	R²	E_p
Uncorrected	3006.9±2588.0	-23.68±4.35	0.44±0.32	20714.2±19455.8
MARA	76.4±72.5	-6.01±7.19	0.77±0.20	603.6±640.7
TDDR	11.9±13.5	1.60±5.74	0.88±0.12	83.4±91.5
Spline-SG	23.6±30.4	-1.37±5.31	0.86±0.10	170.1±206.3
kWavelet	95.2±98.7	-7.44±6.12	0.50±0.33	566.7±622.3
tPCA	78.3±78.6	-6.66±6.05	0.74±0.19	613.8±777.0
Mor	1.7±1.5	8.21±3.43	0.83±0.17	8.2±14.9
tMed	2992.4±2580.1	-23.67±4.32	0.44±0.32	20270.2±19240.6
tMedMor	1.6±1.4	8.66±3.29	0.85±0.16	5.1±6.6

查看所有表

采用不同算法进行运动伪迹校正后,计算得到的真实块平均HbO信号的评价指标。(a) EMS; (b) RSN; (c) R2; (d) Ep

图 8. 采用不同算法进行运动伪迹校正后,计算得到的真实块平均HbO信号的评价指标。(a) E_MS; (b) R_SN; (c) R²; (d) E_p

Fig. 8. Evaluation indices calculated after correcting motion artifacts of true block-average HbO signal with different algorithms. (a) E_MS; (b) R_SN; (c) R²; (d) E_p

下载图片查看所有图片

4.3 讨论

tMed算法适合去除尖峰,但不能有效去除基线突变和缓慢漂移;Mor算法适合去除基线突变和缓慢漂移,但处理距离较近的两个尖峰时会引入基线突变。为充分发挥这两种算法的优势,本文采用两者结合的方法——tMedMor算法去除fNIRS信号中的运动伪迹。采用仿真数据和实验数据进行了验证,结果表明,与MARA、tPCA、kWavelet、Spline-SG、TDDR算法相比,tMedMor算法的去运动伪迹性能良好,可在fNIRS信号预处理时使用。

TDDR算法比tMedMor算法的性能稍差,这是因为TDDR算法适合去除时域微分较大的运动伪迹^[9],而低幅值缓慢漂移的时域微分较小,不适合用TDDR算法进行处理;而tMedMor算法通过选择合适的参数可以将与HRF信号频率范围不重叠的缓慢漂移去除。MARA、Spline-SG算法^[17,19]采用样条曲线拟合运动伪迹,并从信号中减去该曲线(目的是去运动伪迹),这可能会导致部分HRF信号也被去除,使信号校正过度;而本文的数学形态学方法可以很好地跟踪基线突变和缓慢漂移的形状,并保留信号的细节。MARA算法适合去除基线突变,但去尖峰的效果较差^[19];tMedMor算法能够依次去除三种运动伪迹。所以,与MARA相比,tMedMor的性能有了很大提升。tPCA算法^[9]在4个评价指标上的表现稍差,这是因为该算法依赖于运动伪迹的检测,而不同数据集的最佳伪迹检测参数不同,使用原文献中的参数处理本文数据可能欠妥。kWavelet算法得到的R²指标较差,该方法在处理基线突变和缓慢漂移时会使信号拖尾或平滑,而这会使数据分析结果产生偏差^[9]。

当运动伪迹的频率范围与HRF信号重叠时,tMedMor算法的处理效果较差。未来拟将该算法与其他算法相结合,以弥补该算法的不足。Mor算法只用了直线型结构元素,直线型结构元素对基线突变和缓慢漂移的去除效果较好,但在去除尖峰时会使信号变得不平滑,不过这对数据处理结果的影响不大。下一步可以尝试使用不同形状的结构元素去除尖峰,以改善该算法的处理效果。由于实际数据中无法找到完全无噪声的信号,所以Mor算法的参数选取方法仅作为理论参考,实际选取时可以根据实验任务的周期以及运动伪迹的时间宽度进行综合考虑,未来可以在此方面进行进一步研究,以提高算法的实用性。

5 结论

本文对中值滤波算法进行改进,提出了将定向中值滤波(tMed)和数学形态学(Mor)相结合的算法——tMedMor算法,然后用该算法对fNIRS信号中的运动伪迹进行去除,即:用tMed算法去除尖峰,用Mor算法去除基线突变和缓慢漂移。仿真数据和实验数据的验证结果表明,与其他算法相比,tMedMor算法在E_MS、R_SN、E_p三个评价指标上表现得最好,R²与最优值相比也相差不大。总之,tMedMor算法可在fNIRS信号预处理阶段使用。

参考文献

[1] Lloyd-Fox S, Blasi A, Elwell C E. Illuminating the developing brain: the past, present and future of functional near infrared spectroscopy[J]. Neuroscience & Biobehavioral Reviews, 2010, 34(3): 269-284.

[2] Boas D A, Elwell C E, Ferrari M, et al. Twenty years of functional near-infrared spectroscopy: introduction for the special issue[J]. NeuroImage, 2014, 85: 1-5.

[3] Villringer A, Chance B. Non-invasive optical spectroscopy and imaging of human brain function[J]. Trends in Neurosciences, 1997, 20(10): 435-442.

[4] Boas D A, Dale A M, Franceschini M A. Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy[J]. NeuroImage, 2004, 23: S275-S288.

[5] Karim H, Schmidt B, Dart D, et al. Functional near-infrared spectroscopy (fNIRS) of brain function during active balancing using a video game system[J]. Gait & Posture, 2012, 35(3): 367-372.

[6] Tuscan L A, Herbert J D, Forman E M, et al. Exploring frontal asymmetry using functional near-infrared spectroscopy: a preliminary study of the effects of social anxiety during interaction and performance tasks[J]. Brain Imaging and Behavior, 2013, 7(2): 140-153.

[7] Brigadoi S, Ceccherini L, Cutini S, et al. Motion artifacts in functional near-infrared spectroscopy: a comparison of motion correction techniques applied to real cognitive data[J]. NeuroImage, 2014, 85: 181-191.

[8] Scholkmann F, Metz A J, Wolf M. Measuring tissue hemodynamics and oxygenation by continuous-wave functional near-infrared spectroscopy: how robust are the different calculation methods against movement artifacts?[J]. Physiological Measurement, 2014, 35(4): 717-734.

[9] Fishburn F A, Ludlum R S, Vaidya C J, et al. Temporal derivative distribution repair (TDDR): a motion correction method for fNIRS[J]. NeuroImage, 2019, 184: 171-179.

[10] Cooper R J, Selb J, Gagnon L, et al. A systematic comparison of motion artifact correction techniques for functional near-infrared spectroscopy[J]. Frontiers in Neuroscience, 2012, 6: 147.

[11] Izzetoglu M, Chitrapu P, Bunce S, et al. Motion artifact cancellation in NIR spectroscopy using discrete Kalman filtering[J]. BioMedical Engineering OnLine, 2010, 9(1): 1-10.

[12] Robertson F C, Douglas T S, Meintjes E M. Motion artifact removal for functional near infrared spectroscopy: a comparison of methods[J]. IEEE Transactions on Biomedical Engineering, 2010, 57(6): 1377-1387.

[13] Gagnon L, Yücel M A, Boas D A, et al. Further improvement in reducing superficial contamination in NIRS using double short separation measurements[J]. NeuroImage, 2014, 85: 127-135.

[14] BlasiA, PhillipsD, Lloyd-FoxS, et al. and Biology. [S.l.:s.n.], 2010: 279-284. https:∥link.springer.com/chapter/10.1007%2F978-1-4419-1241-1_40.

[15] Chiarelli A M, Maclin E L, Fabiani M, et al. A kurtosis-based wavelet algorithm for motion artifact correction of fNIRS data[J]. NeuroImage, 2015, 112: 128-137.

[16] Molavi B, Dumont G A. Wavelet-based motion artifact removal for functional near-infrared spectroscopy[J]. Physiological Measurement, 2012, 33(2): 259-270.

[17] Scholkmann F, Spichtig S, Muehlemann T, et al. How to detect and reduce movement artifacts in near-infrared imaging using moving standard deviation and spline interpolation[J]. Physiological Measurement, 2010, 31(5): 649-662.

[18] Yucel M A, Selb J, Cooper R J, et al. Targeted principle component analysis: a new motion artifact correction approach for near-infrared spectroscopy[J]. Journal of Innovative Optical Health Sciences, 2014, 7(2): 1350066.

[19] Jahani S, Setarehdan S K, Boas D A, et al. Motion artifact detection and correction in functional near-infrared spectroscopy: a new hybrid method based on spline interpolation method and Savitzky-Golay filtering[J]. Neurophotonics, 2018, 5(1): 015003.

[20] 邢承滨, 邓兴升, 徐康. 形态学滤波关键参数阈值的等值线确定方法[J]. 激光与光电子学进展, 2019, 56(16): 162802.

Xing C B, Deng X S, Xu K. Contour determination method for threshold of morphological filtering key parameters[J]. Laser & Optoelectronics Progress, 2019, 56(16): 162802.

[21] 刘姝. 数学形态学在信号处理方面的应用研究[D]. 大连: 大连理工大学, 2006: 4- 5.

LiuS. Mathematical morphology and its application on signal processing[D]. Dalian: Dalian University of Technology, 2006: 4- 5.

[22] 季虎, 孙即祥, 毛玲. 基于小波变换与形态学运算的ECG自适应滤波算法[J]. 信号处理, 2006, 22(3): 333-337.

Ji H, Sun J X, Mao L. An ECG adaptive filter algorithm based on wavelet transform and mathematical morphology[J]. Signal Processing, 2006, 22(3): 333-337.

[23] Wang Z, Zhang D. Progressive switching median filter for the removal of impulse noise from highly corrupted images[J]. IEEE Transactions on Circuits and Systems, 1999, 46(1): 78-80.

[24] Wang J H, Lin L D. Improved median filter using minmax algorithm for image processing[J]. Electronics Letters, 1997, 33(16): 1362-1363.

[25] 朱湘临, 杨建宁, 唐平, 等. 心电图波形基线漂移的数学形态滤波在线实时处理[J]. 生物医学工程学杂志, 2010, 27(1): 48-52.

Zhu X L, Yang J N, Tang P, et al. Baseline drift of electrocardiograph dealt with mathematical morphology filter in real-time[J]. Journal of Biomedical Engineering, 2010, 27(1): 48-52.

[26] Buxton R B, Wong E C, Frank L R. Dynamics of blood flow and oxygenation changes duringbrain activation: the balloon model[J]. Magnetic Resonance in Medicine, 1998, 39(6): 855-864.

[27] Leamy DJ, Ward TE, Sweeney KT. Functional near infrared spectroscopy (fNIRS) synthetic data generation[C]∥2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society. August 30-September 3, 2011, Boston, MA, USA.New York: IEEE Press, 2011: 6589- 6592.

赵杰, 乔吉日木图, 丁雪桐, 梁晓敏. 基于数学形态学和中值滤波的fNIRS信号运动校正算法研究[J]. 光学学报, 2020, 40(22): 2230002. Jie Zhao, Jirimutu Qiao, Xuetong Ding, Xiaomin Liang. fNIRS Signal Motion Correction Algorithm Based on Mathematical Morphology and Median Filter[J]. Acta Optica Sinica, 2020, 40(22): 2230002.

基于数学形态学和中值滤波的fNIRS信号运动校正算法研究下载： 1150次

1 引言

2 理论与算法

图 1. 采用tMedMor算法对fNIRS信号进行预处理的流程图

Fig. 1. fNIRS signal preprocessing flowchart using tMedMor algorithm

2.1 tMed算法

2.2 Mor算法

图 2. Mor算法流程图

Fig. 2. Flow chart of Mor algorithm

3 实验部分

3.1 仿真验证

图 3. 无噪OD信号的单边功率谱密度

Fig. 3. Unilateral power spectral density of noiseless OD signal

3.2 实验验证

图 4. 无噪OD信号的单边功率谱密度

Fig. 4. Unilateral power spectral density of noiseless OD signal

3.3 对比方法

3.4 评价指标

4 结果与讨论

4.1 仿真验证

表 1. 经不同算法处理后,仿真块平均HbO信号的E_MS、R_SN、R²和E_p

Table 1. E_MS, R_SN, R² and E_p of simulated block-average HbO signal after processing by different algorithms

图 6. 采用不同算法进行运动伪迹校正后,计算得到的仿真块平均HbO信号的评价指标。(a) E_MS; (b) R_SN; (c) R²; (d) E_p

Fig. 6. Evaluation indices calculated after correcting motion artifacts of simulated block-average HbO signal with different algorithms. (a) E_MS; (b) R_SN; (c) R²; (d) E_p

4.2 实验验证

表 2. 经不同算法处理后,真实块平均HbO信号的E_MS、R_SN、R²和E_p

Table 2. E_MS, R_SN, R² and E_p of real block-average HbO signal after processing by different algorithms

图 8. 采用不同算法进行运动伪迹校正后,计算得到的真实块平均HbO信号的评价指标。(a) E_MS; (b) R_SN; (c) R²; (d) E_p

Fig. 8. Evaluation indices calculated after correcting motion artifacts of true block-average HbO signal with different algorithms. (a) E_MS; (b) R_SN; (c) R²; (d) E_p

4.3 讨论

5 结论

Article Outline

关于本站 Cookie 的使用提示

全站搜索

基于数学形态学和中值滤波的fNIRS信号运动校正算法研究 下载： 1150次

1 引言

2 理论与算法

图 1. 采用tMedMor算法对fNIRS信号进行预处理的流程图

Fig. 1. fNIRS signal preprocessing flowchart using tMedMor algorithm

2.1 tMed算法

2.2 Mor算法

图 2. Mor算法流程图

Fig. 2. Flow chart of Mor algorithm

3 实验部分

3.1 仿真验证

图 3. 无噪OD信号的单边功率谱密度

Fig. 3. Unilateral power spectral density of noiseless OD signal

3.2 实验验证

图 4. 无噪OD信号的单边功率谱密度

Fig. 4. Unilateral power spectral density of noiseless OD signal

3.3 对比方法

3.4 评价指标

4 结果与讨论

4.1 仿真验证

表 1. 经不同算法处理后,仿真块平均HbO信号的EMS、RSN、R2和Ep

Table 1. EMS, RSN, R2 and Ep of simulated block-average HbO signal after processing by different algorithms

图 6. 采用不同算法进行运动伪迹校正后,计算得到的仿真块平均HbO信号的评价指标。(a) EMS; (b) RSN; (c) R2; (d) Ep

Fig. 6. Evaluation indices calculated after correcting motion artifacts of simulated block-average HbO signal with different algorithms. (a) EMS; (b) RSN; (c) R2; (d) Ep

4.2 实验验证

表 2. 经不同算法处理后,真实块平均HbO信号的EMS、RSN、R2和Ep

Table 2. EMS, RSN, R2 and Ep of real block-average HbO signal after processing by different algorithms

图 8. 采用不同算法进行运动伪迹校正后,计算得到的真实块平均HbO信号的评价指标。(a) EMS; (b) RSN; (c) R2; (d) Ep

Fig. 8. Evaluation indices calculated after correcting motion artifacts of true block-average HbO signal with different algorithms. (a) EMS; (b) RSN; (c) R2; (d) Ep

4.3 讨论

5 结论

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

基于数学形态学和中值滤波的fNIRS信号运动校正算法研究下载： 1150次

表 1. 经不同算法处理后,仿真块平均HbO信号的E_MS、R_SN、R²和E_p

Table 1. E_MS, R_SN, R² and E_p of simulated block-average HbO signal after processing by different algorithms

图 6. 采用不同算法进行运动伪迹校正后,计算得到的仿真块平均HbO信号的评价指标。(a) E_MS; (b) R_SN; (c) R²; (d) E_p

Fig. 6. Evaluation indices calculated after correcting motion artifacts of simulated block-average HbO signal with different algorithms. (a) E_MS; (b) R_SN; (c) R²; (d) E_p

表 2. 经不同算法处理后,真实块平均HbO信号的E_MS、R_SN、R²和E_p

Table 2. E_MS, R_SN, R² and E_p of real block-average HbO signal after processing by different algorithms

图 8. 采用不同算法进行运动伪迹校正后,计算得到的真实块平均HbO信号的评价指标。(a) E_MS; (b) R_SN; (c) R²; (d) E_p

Fig. 8. Evaluation indices calculated after correcting motion artifacts of true block-average HbO signal with different algorithms. (a) E_MS; (b) R_SN; (c) R²; (d) E_p