一种多策略的NIRS干扰检测和去除方法

摘要

一种多策略的NIRS干扰检测和去除方法，包含三个步骤：步骤一，干扰检测，用NIRS信号与加速度信号对干扰进行检测，将干扰分为三种类型，即轻度干扰、中度干扰、重度干扰，并在此基础上确认信号是否发生了基线突变；步骤二，严重干扰处理阶段，将中度干扰、重度干扰和基线突变用插值的方法进行去除；步骤三，小波双阈值干扰去除，得到干净的NIRS信号。本发明采用多种策略来检测NIRS信号中的干扰，适合于短时测量和长时监测等多个应用场景，同时将干扰分为不同的类型并予以不同的处理方式，最后根据检测阈值和处理阈值的双保险策略，可以在有效去除干扰的同时尽可能多地保留有用信息，相比于传统的NIRS干扰去除算法，具有准确性和可靠性都得到了提高。

说明书

技术领域

本发明属于生物医学信号处理技术领域，是针对近红外光谱信号的一种多策略的NIRS干扰检测和去除方法。

背景技术

近红光谱技术(NIRS)是监测脑功能的一个有效方法。NIRS可以反映目标区域的血液动力学信息，通过获取经过人体组织散射的光信号，再利用朗伯-比尔定律进行计算，就可以得到相应组织的含氧血红蛋白的浓度变化(Δ[HbO₂])和脱氧血红蛋白的浓度变化(Δ[Hb])等相关信息。NIRS相对于功能磁共振成像(fMRI)、计算机断层扫描(CT)等传统方法有着无创监测、低成本和易操作性等优点，使之可以应用于许多新的领域，如康复、脑机接口等，同样也被广泛应用于人类睡眠、社会认知等科学研究领域。

然而，NIRS信号(包括含氧血红蛋白浓度变化的信号、脱氧血红蛋白浓度变化的信号等)对于运动干扰非常敏感，这些干扰通常是由光测量探头和被测组织之间的相对位移引起的。监测对象大幅度运动或突然的抖动往往会显著改变探头和被测组织之间的光耦合过程，从而在NIRS光信号中引入噪声，又因为从光信号得到血氧信号是一个非线性的计算过程，最终计算得到的血氧浓度变化信号将会含有更多不可控的干扰信号。而被测对象不同的运动形式会在NIRS信号中造成不同形态的运动伪迹，以前额佩戴的NIRS监测系统为例，快速的头部晃动会在NIRS信号中引入大幅度、高频率的干扰；而缓慢的头部转动会引起持续的低幅低频干扰；同时因为头位的改变，还会伴随有信号的基线突变情况。当受试者的身体状态无法被控制，监测部位的运动不可避免时(如睡眠监测)，运动干扰的问题就显得尤为突出，多种多样的运动干扰会严重污染光信号，使得之后的信号分析工作受到影响。

在分析NIRS信号之前，首要的任务就是找到一种合适的干扰去除算法。目前有一系列的常用干扰去除算法，但是各自都存在不可避免的缺点。一种常用的做法是舍去包含运动干扰的整段信号，但是这种方法在去除干扰的同时，其中包含的有用信号也会被去除。自适应滤波技术可以通过引入一个参考信号来减小运动干扰的影响，但参考信号的选择却增加了额外的难度。维纳滤波和卡尔曼滤波在有些情况下会取得比较好的效果，但是这两种方法都需要对所处理的信号有一个比较好的先验认识。小波分解被广泛尝试应用于NIRS信号的干扰去噪过程，Molavi和Dumont采用离散小波去噪的方法，假设信号的小波系数是正态分布，从而判断出其中的离群值并处理，再进行小波重构得到去噪处理后的NIRS信号，取得了比较好的干扰去除效果；但是因为在离散小波离群值判断时使用单一阈值进行判断，对于较为干净的信号，这种强制去除一定比例小波系数的方法反而会失去一些有用的信息；而且，小波分解对于头动引起的大的基线突变是无能为力的。Virtanen等人提出一种基于加速度信号的干扰去除算法，这种方法可以去除NIRS信号中基线突变的趋势，但是他们使用干扰片段的平均值来代替整段信号的值，毫无疑问会使得干扰段内的有用信息被忽略掉。Scholkmann等人尝试了样条插值的方法对NIRS信号进行干扰去除，但是很难找到一个适合各种各样运动干扰的全局阈值，并且这个方法中还有很多其他的参数需要不断调整，这些都增加了方法使用的复杂性。

综上所述，目前对NIRS信号进行干扰去除的方法不能完全应对所有可能出现的运动干扰情况。在短时间监测中，监测对象的运动状态可以被有效控制，此时离散小波去噪、样条插值等方法可以取得较好的干扰去除结果；然而在长时间监测中，尤其是在睡眠监测这样的条件下，监测对象的运动状态无法预估且不能被控制，NIRS信号中的干扰形式将会多种多样，不仅会有各种幅度大小不一的瞬时尖峰，随着监测对象的转头或翻身还会伴随有信号基线的突变，这就使得原来的一些干扰去除方法无法再满足现有的情况。

发明内容

为了使NIRS技术可以普遍应用于各个领域，需要一种具有广泛适用性的NIRS干扰去除方法，克服上述现有干扰去除技术的缺陷，可以应对NIRS信号中可能存在的各种干扰，因此本发明提出了一种多策略的NIRS干扰检测和去除方法，兼顾了采集到的NIRS信号和加速度信号，考虑了可能包含于NIRS信号中的几乎所有类型的运动干扰，依据信号中干扰的严重程度将干扰分为三种：轻度干扰、中度干扰和重度干扰，同时，在干扰区域内也判断是否存在基线突变的情况。本发明综合了基于加速度的干扰去除、样条插值法去噪、离散小波去噪等方法，并改进了其中的缺陷。本发明提出的干扰去除方法可以很好地应用于长时间监测的信号中，尤其是针对整晚睡眠的情况，对于短时间监测信号，也有显著效果。

为了达到上述目的，本发明采用的技术方案为：

一种多策略的NIRS干扰检测和去除方法，包括以下步骤：

步骤一：干扰检测

本方法在检测干扰类型时，综合考虑了NIRS信号和加速度信号。

首先计算被检测信号即NIRS信号或者加速度信号的移动标准差MSD，然后根据MSD序列可得到一个随着时间变化的动态阈值序列，根据检测信号的MSD序列的值与对应动态阈值的关系，检测出不同类型的干扰；再根据干扰检测的结果及NIRS信号的变化趋势，判断出干扰区域内是否存在基线突变。

①首先，根据NIRS信号对干扰进行初步检测。

由NIRS信号计算得到其MSD序列，根据MSD序列的分布情况得到一个随时间变化的动态阈值序列。然后将MSD序列与动态阈值序列进行比较，如果MSD序列的某个值超过相应动态阈值的4倍，将其标记为疑似重度干扰区；MSD序列在动态阈值的2倍和4倍之间，则对应位置标记为疑似中度干扰区；如果MSD序列在对应动态阈值的1倍和2倍之间，则对应位置标记为轻度干扰区，其余认为无干扰。

紧接着对疑似干扰区进行确认。对于前面标记的疑似重度干扰区和疑似中度干扰区，需要运用原始NIRS信号进行干扰确认。根据斜率的变化对NIRS信号的疑似干扰区进行分割，将其分割成若干个相邻小段，保证分割后的每个小段内有且仅有一个波峰，且不含波谷。如果每个小段内NIRS信号的极差(最大值减去最小值)超过整个干扰区信号标准差的6倍，则确认疑似重度干扰为重度干扰，疑似中度干扰为中度干扰，否则认为该小段信号为轻度干扰；

②其次，判断NIRS信号的干扰区内是否发生基线突变。

在NIRS信号的干扰类型确认之后，需要对NIRS信号的基线突变情况进行检测。使用原始NIRS信号，结合干扰的检测结果，判断在干扰区内是否发生了基线突变：如果干扰区前后非干扰区的信号长度均超过一定时间，则分别计算干扰区前后两段非干扰区NIRS信号的均值，若这两个均值的差超过了干扰区内NIRS信号最大值与最小值之差的一半，则认为这个干扰区的前后存在基线突变，将对应干扰区域标记为存在基线突变。

③再次，对加速度信号进行干扰检测。

其过程和NIRS信号的检测过程相似，同样也会检测出三种类型的干扰：轻度干扰、中度干扰、重度干扰，不同于NIRS信号检测过程的是，利用加速度信号的MSD序列与动态阈值进行比较时，可直接确认轻、中、重度干扰，而不再设立疑似中、重度干扰区。加速度信号的干扰检测结果仅作为辅助信号，用于对NIRS信号的干扰检测结果进行确认。

④综合NIRS信号和加速度信号的干扰检测结果进行干扰确认。

由NIRS信号和加速度信号的干扰检测过程可以得到不同的干扰检测结果。在干扰确认过程中，以NIRS信号的干扰检测结果为基础，用加速度信号的检测结果进行辅助判断。若NIRS信号的干扰检测结果中检测出了重度干扰，则无论对应位置加速度的检测结果如何，均确认NIRS信号中有重度干扰；若NIRS信号中检测出有轻度或中度干扰，此时若对应位置加速度信号的检测结果中有干扰(轻、中、重度干扰)，则认为NIRS信号中的轻度或中度干扰存在(干扰类型不变)，否则认为NIRS信号在此处没有干扰；若NIRS信号中未检测出干扰，则无论此处加速度信号的干扰检测结果如何，均认为此处不含有干扰。

⑤最后，进行基线突变区域的确认。

综合NIRS信号的基线突变检测结果与加速度信号的干扰检测结果，确认NIRS信号发生基线突变的位置。以NIRS信号的基线突变检测结果为基础，若NIRS信号中检测出了基线突变，同时加速度信号在对应位置检测出了干扰(轻度、中度或重度干扰)，则认为该处存在基线突变，若加速度信号在对应位置没有检测出干扰，则认为NIRS信号在该处不存在基线突变(应变更之前的检测结果)。

步骤二：严重干扰处理

通过步骤一的干扰检测，可以检测出轻度干扰、中度干扰和重度干扰，在此基础上还检测出了发生基线突变的位置。在这里定义中度干扰、重度干扰和基线突变属于严重干扰，这些正是这一步骤需要处理的。在步骤二中，使用三次平滑样条插值的方法，首先处理中度干扰和重度干扰，再处理基线突变。对于轻度干扰区域，则直接跳过步骤二而在步骤三中进行去除。

在中度干扰和重度干扰的去除时，对每一个干扰区进行三次平滑样条插值(平滑参数的设置值逼近1，类似三次样条插值)，并用原始NIRS信号减去插值后的信号得到干扰去除后的信号。然后对于基线突变进行去除，对含有基线突变的干扰区同样进行三次平滑样条插值(平滑参数的设置值逼近0，类似线性插值)，用原始NIRS信号减去对干扰区进行插值得到的信号后，还需要在这个信号的基础上加上相对于此基线突变干扰区开始时的直流偏移量，以保证信号的连续性。

至此，已经从原始NIRS信号中去除了严重干扰的部分，在去除干扰引起的大幅度信号波动的同时，也尽可能多地保留了细节的有用信息。

步骤三：小波双阈值处理方法

对于NIRS信号，运动干扰的存在通常会在离散小波系数中呈现出离群值，利用双阈值的离散小波去噪法来处理剩下的运动干扰，采用两个重要参数w₁和w₂来判断离群值，其中w₁称作检测阈值，用来检测哪些小波系数应被视为运动干扰，而w₂称作处理阈值，略小于w₁，用来决定具体哪些小波系数应该被处理。而双阈值的概念即是综合利用检测阈值和处理阈值来检测小波系数中的离群点。

首先对需要进行干扰去除的信号进行离散小波分解得到小波系数，并估计其标准差，然后根据小波系数的标准差及概率参数α来确定双阈值的值。对于离散小波变换得到的小波系数，逐个比较其绝对值是否大于w₁，如果是，则说明该小波系数中包含了运动干扰的情况；然后在此超过检测阈值w₁的小波系数附近，再次检测大于绝对值大于w₂的小波系数，一直到其左右两侧小波系数的绝对值小于w₂为止，标记在这个过程中检测到的所有值大于w₂的小波系数，并将值其置为0，而其余的小波系数的值则保持不变。最后用处理后的小波系数去重构得到希望的NIRS信号。

双阈值的确定方法如下：假设某一层离散小波系数的长度为L，将这些小波系数进行升序排列，然后找到距离0.1587L位置点最近的小波系数的值记为表示的则是对这些小波系数的标准差估计值；之后找到概率参数α在正态分布中对应的上分位数的值u，处理阈值w₂的值即等于检测阈值应略大于处理阈值。

所述的利用双阈值的离散小波去噪法来处理剩下的运动干扰，方法如下：对于离散小波变换得到的每一层的小波系数，逐个比较其绝对值是否大于w₁，如果是，则说明该小波系数中包含了运动干扰的情况；然后在此超过检测阈值w₁的小波系数附近，再次检测大于绝对值大于w₂的小波系数，一直到其左右两侧小波系数的绝对值小于w₂为止，标记在这个过程中检测到的所有值大于w₂的小波系数，并将值其置为0，而其余的小波系数的值则保持不变；最后用处理后的每一层的小波系数去重构得到最终干净的NIRS信号。

为了避免边界处出现吉布斯或者突变现象，可借鉴循环平均的思想，这样可以减少由于离散小波变换带来的尖峰干扰。

本发明的优点主要有四点：

一、同时利用了由近红外监测探头得到的加速度信号和NIRS信号进行干扰检测，适用于多种测量场景。本发明在步骤一的干扰检测中，综合利用了加速度信号和NIRS信号的干扰检测结果，而在实际应用中，若无法获取加速度传感器信号，或监测对象的运动状态可以加以控制，也可单独选择NIRS信号进行干扰检测。例如在短时间测量中，加速度信号中的干扰是可控的，所以可以只采用基于NIRS信号的干扰检测算法；在长时间的监护过程中，可以辅助以加速度信号来准确检测运动干扰。

二、在运动干扰检测中采用动态阈值，考虑了NIRS信号中大多数的运动干扰。随时间变化的动态阈值可以更好地适应长时间NIRS信号的干扰检测过程，使得干扰定位过程更合理、准确。根据干扰强度的不同，本发明将NIRS信号中干扰分为三种：轻度干扰、中度干扰、重度干扰，并且还将信号的基线突变情况考虑在内，这些基本囊括了实际测量过程中被试的转头、瞬时抖动、夜间翻身等各种干扰事件。对于不同类型的干扰，本发明规定了不同大小的动态阈值范围来进行检测，以更好地区分不同强度的运动干扰。

三、利用样条平滑插值方法对中度/重度干扰、基线突变进行干扰处理。针对中度/重度干扰的去除，采用平滑参数接近于1的样条插值过程可以拟合出干扰区域内信号的大幅波动情况，用原信号减去拟合结果可以在去除运动干扰引起的大波动的同时保留信号本身的细节变化信息；针对基线突变的去除，则采用平滑参数接近于0的样条插值过程，可以较好拟合出基线突变时的趋势变化情况，用经过中度和重度干扰去除后的NIRS信号减去基线突变的拟合结果，可在去除基线突变的同时，将基线变化过程中的细节波动信息保留下来。

四、创新性地提出了离散小波双阈值去噪的方法。目前常用的离散小波去噪方法中只有一个阈值决定小波系数中的离群值，在这样的情况下如果概率参数α取值不合理，将会导致丢失很多有用的信息。本发明在去除严重干扰之后，提出了双阈值的离散小波，检测阈值略大于处理阈值，这样可以尽可能准确地处理干扰区域，同时保留有效的信号。

附图说明

图1是本发明算法的整体框图。

图2是利用NIRS信号的进行干扰检测过程的详细描述。

图3是利用加速度信号的进行干扰检测过程的详细描述。

图4是综合NIRS信号和加速度信号的干扰检测结果的详细描述。

图5是严重干扰处理过程的详细描述。

图6是小波双阈值去除干扰算法的详细描述。

图7是仿真产生的NIRS信号和加速度信号，其中图7a是仿真产生的NIRS信号，图7b加速度信号。

图8是由NIRS信号计算得到的MSD序列及T₁、T₃的图形描述。

图9是用T₁、T₃处理原始序列并滤波之后的图形描述。

图10是最终得到的动态阈值与初始MSD序列及T₂的图形描述。

图11是对仿真信号进行干扰检测的初步结果。

图12是对仿真信号进行干扰确认后的最终结果。

图13是对仿真NIRS信号进行严重干扰去除后的图示。

图14是小波双阈值及原始信号小波系数的图示。

图15是经小波双阈值去噪法处理后的NIRS信号和原始NIRS信号对比图示。

具体实施例

图1-图6的流程图系统地阐述了本发明的操作方法，为了更加清楚说明本发明的操作过程，以下以仿真得到的含氧血红蛋白浓度变化信号(简称：Δ[HbO₂])为例对方法和流程进行详细描述。

1.仿真信号：

仿真信号包括干净的含氧血红蛋白浓度变化信号x(n)、含有干扰的含氧血红蛋白浓度变化信号X(n)，与加速度传感器信号a(n)。时间长度为240s，采样频率为fs＝10Hz。

仿真的不含运动干扰的Δ[HbO₂]信号根据以下公式产生：

其包含四个不同频率的正弦波，其中正弦波表示1Hz的心率成分、0.25Hz的呼吸成分、0.1hz的脉叶波成分和0.04Hz的低频全局振荡成分，而ε(n)表示模拟系统噪声的白噪声信号

含有运动干扰的Δ[HbO₂](记作Δ[HbO₂]-N)可由以下公式表示：

X(n)＝x(n)+δ_peak(n)+δ_shift(n)>

其中δ_peak(n)模拟的是头部抖动引起的瞬时尖峰信号，δ_shift(n)则表示头部转动引起的基线突变信号。

a₁(n)，a₂(n)，a₃(n)分别表示x、y、z三个方向的仿真加速度传感器信号，单位为g，均是是δ_peak(n)和δ_shift(n)的组合结果，在幅值上有差异。

仿真信号的时域波形如图7所示，包括Δ[HbO₂]信号和加速度信号。

一种多策略的NIRS干扰检测和去除方法，包括以下步骤：

步骤一：干扰检测。

本方法在检测干扰类型时，综合考虑了NIRS信号和加速度信号。

首先计算被检测信号(NIRS信号或者加速度信号)的移动标准差(MSD)，然后根据MSD序列可得到一个随着时间变化的动态阈值序列。根据检测信号的MSD序列的值与对应动态阈值的关系，可以检测出不同类型的干扰；再根据干扰检测的结果及NIRS信号的变化趋势，可以判断出干扰区域内是否存在基线突变。

具体为：给定一个信号z(n)，信号长度为N，双边MSD序列s(n)可以按照以下公式计算得到：

这里n＝k+1,k+2,...,N-k，W＝2k+1，表示滑动窗的长度。

首先根据公式(4)将各方向加速度传感器信号整合，得到综合加速度信号序列acc(n)：

其中，M＝3，a_i(n)可代表三个通道的加速度信号。

再根据公式(3)分别计算得到信号X(n)与综合加速度传感器信号acc(n)的MSD序列t(n)和d(n)，其中k＝3fs。

然后，利用序列t(n)和d(n)进行干扰检测。

①Δ[HbO₂]-N信号的干扰检测。

由NIRS信号计算得到其MSD序列，根据MSD序列的分布情况得到一个随时间变化的动态阈值序列。然后将MSD序列与动态阈值序列进行比较，如果MSD序列的某个值超过相应动态阈值的4倍，将其标记为疑似重度干扰区；MSD序列在动态阈值的2倍和4倍之间，则对应位置标记为疑似中度干扰区；如果MSD序列在对应动态阈值的1倍和2倍之间，则对应位置标记为轻度干扰区，其余认为无干扰。

紧接着对疑似干扰区进行确认。对于前面标记的疑似重度干扰区和疑似中度干扰区，需要运用原始NIRS信号进行干扰确认。根据斜率的变化对NIRS信号的疑似干扰区进行分割，将其分割成若干个相邻小段，保证分割后的每个小段内有且仅有一个波峰，且不含波谷。如果每个小段内NIRS信号的极差(最大值减去最小值)超过整个干扰区信号标准差的6倍，则确认疑似重度干扰为重度干扰，疑似中度干扰为中度干扰，否则认为该小段信号为轻度干扰；

具体为：

将Δ[HbO₂]-N信号的MSD序列t(n)按数值从小到大的方式进行排列可得序列L(n)。记序列L(n)第30％的位置对应的数据为T₁，第50％的位置对应的数据为T₂，第70％的位置对应的数据为T₃(如图8所示)。判断MSD序列t(n)对应位置的值与T₁和T₃值的大小关系，如果t(n)相应的值大于T₃，则将其置为T₃；若t(n)相应的值小于T₁，则将将其置为T₁(如图9所示)。接下来对重置后的t(n)进行一个较低频率的低通滤波，得到t′(n)，为t′(n)序列加上直流偏量3(T₃-T₁)从而得到序列t″(n)。最后将t″(n)序列中所有小于T₂的值置为T₂，即可得到最终的动态阈值A(n)序列(如图10所示)。

然后进行干扰区域的判断。将MSD序列t(n)与动态阈值序列A(n)进行比较，用t_k表示t(n)序列中的任意一点，A_k代表A(n)中与t_k位置对应的一个点。如果t_k超过4A_k，将其标记为疑似重度干扰区；如果t_k在2A_k和4A_k之间，则对应位置标记为疑似中度干扰区；果t_k在A_k和2A_k之间，则对应位置标记为轻度干扰区；其余情况下，则认为此处不含有干扰。

之后对疑似干扰事件进行确认。对于前面标记的疑似重度干扰区和疑似重度干扰区，需要运用原始Δ[HbO₂]-N信号X(n)进行干扰确认。根据斜率的变化对X(n)的疑似干扰区进行分割，将其分割成若干个相邻小段，保证段与段的交接点在波谷位置，分割后的每个小段内只包含一个波峰。如果某个小段内X(n)的极差(最大值减去最小值)超过整个疑似干扰区信号标准差的6倍，则确认此疑似重度干扰区内的该小段信号为重度干扰，相应的疑似中度干扰区内的该小段信号为中度干扰；若不满足超过6倍这一条件，则认为该小段信号为轻度干扰。

②其次，判断NIRS信号的干扰区内是否发生基线突变。

在NIRS信号的干扰类型确认之后，需要对NIRS信号的基线突变情况进行检测。使用原始NIRS信号，结合干扰的检测结果，判断在干扰区内是否发生了基线突变：如果干扰区前后非干扰区的信号长度均超过一定时间，则分别计算干扰区前后两段非干扰区NIRS信号的均值，若这两个均值的差超过了干扰区内NIRS信号最大值与最小值之差的一半，则认为这个干扰区的前后存在基线突变，将对应干扰区域标记为存在基线突变。

具体为：

对Δ[HbO₂]-N信号进行基线突变的检测。

以Δ[HbO₂]-N信号的干扰检测结果为基础，判断在干扰区内是否发生了基线突变。如果干扰区域前后非干扰区的信号长度均超过5s，则分别计算前后两段非干扰区X(n)的均值，若一个干扰事件前后两个均值的差超过了干扰区内X(n)信号最大值与最小值之差的一半，则认为这个干扰区的前后存在基线突变，将对应干扰区域标记为基线突变的发生区域。

③再次，对加速度信号进行干扰检测。

其过程和NIRS信号的检测过程相似，同样也会检测出三种类型的干扰：轻度干扰、中度干扰、重度干扰，不同于NIRS信号检测过程的是，利用加速度信号的MSD序列与动态阈值进行比较时，可直接确认轻、中、重度干扰，而不再设立疑似中、重度干扰区。加速度信号的干扰检测结果仅作为辅助信号，用于对NIRS信号的干扰检测结果进行确认。

具体为：

综合加速度信号的干扰检测。

将综合加速度信号acc(n)的MSD序列d(n)进行升序排列，将大小排列在第30％处的数据点的值记为T，将T的值限定在区间[a,b]中，[a,b]可取为[0.001,0.002]，如果T＜0.001，则令T＝0.001；如果T＞0.002，则令T＝0.002。将d(n)序列的值与T作比较，如果d(n)序列中相应的值大于T，则将其置为T，如果小于则将其置为得到重置之后的序列为d′(n)，最后得到的动态阈值B(n)等于d′(n)的2倍。

用d_k表示d(n)序列中的任意一点，B_k代表B(n)中与d_k位置对应的一个点。如果d_k的值大于4B_k值，则将d_k对应位置的数据标记为重度干扰；若d_k的值在2B_k到4B_k之间，则对应位置标记为中度干扰；若d_k的值在B_k到2B_k之间，则对应位置标记为轻度干扰；否则，认为该位置没有干扰。

图11是对Δ[HbO₂]-N信号的干扰和基线突变检测结果，以及对综合加速度信号的干扰检测结果，阴影的深度不同代表干扰的强度不同。

④综合NIRS信号和加速度信号的干扰检测结果进行干扰确认。

由NIRS信号和加速度信号的干扰检测过程可以得到不同的干扰检测结果。在干扰确认过程中，以NIRS信号的干扰检测结果为基础，用加速度信号的检测结果进行辅助判断。若NIRS信号的干扰检测结果中检测出了重度干扰，则无论对应位置加速度的检测结果如何，均确认NIRS信号中有重度干扰；若NIRS信号中检测出有轻度或中度干扰，此时若对应位置加速度信号的检测结果中有干扰(轻、中、重度干扰)，则认为NIRS信号中的轻度或中度干扰存在(干扰类型不变)，否则认为NIRS信号在此处没有干扰；若NIRS信号中未检测出干扰，则无论此处加速度信号的干扰检测结果如何，均认为此处不含有干扰。

具体为：

结合Δ[HbO₂]-N信号和综合加速度信号的干扰检测结果进行干扰确认。

由信号(Δ[HbO₂]-N信号和综合加速度信号)的干扰检测过程可以得到不同的干扰检测结果。在干扰确认过程中，以Δ[HbO₂]-N信号X(n)的干扰检测结果为基础，用综合加速度信号acc(n)的检测结果进行辅助判断。若X(n)的干扰检测结果中检测出了有重度干扰，则不管综合加速度信号对应位置的检测结果是什么，确认X(n)中有重度干扰；若X(n)中检测出有轻度或中度干扰，此时若对应位置acc(n)的检测结果中有干扰(轻度、中度或重度干扰)，则认为X(n)中的轻度或中度干扰存在(保持原干扰检测类型不变)，否则认为X(n)在此处没有干扰；若X(n)中未检测出干扰，则无论此处综合加速度信号的干扰检测结果如何，均认为此处不含有干扰。

⑤最后，进行基线突变区域的确认。

综合NIRS信号的基线突变检测结果与加速度信号的干扰检测结果，确认NIRS信号发生基线突变的位置。以NIRS信号的基线突变检测结果为基础，若NIRS信号中检测出了基线突变，同时加速度信号在对应位置检测出了干扰(轻度、中度或重度干扰)，则认为该处存在基线突变，若加速度信号在对应位置没有检测出干扰，则认为NIRS信号在该处不存在基线突变(应变更之前的检测结果)。

具体为：对Δ[HbO₂]-N进行基线突变确认。

若X(n)检测出了基线突变，同时综合加速度信号acc(n)在对应位置检测出了干扰(轻、中、重度干扰)，则认为该处存在基线突变，若acc(n)在对应位置出没有检测出无干扰，则认为X(n)在该处不存基线突变。

图12是整合了Δ[HbO₂]-N信号和综合加速度信号的检测结果，进行干扰和基线突变确认后的结果。

步骤二：严重干扰处理，包括重度干扰、中度干扰和基线突变的处理。

通过步骤一的干扰检测，可以检测出轻度干扰、中度干扰和重度干扰，在此基础上还检测出了发生基线突变的位置。在这里定义中度干扰、重度干扰和基线突变属于严重干扰，这些正是这一步骤需要处理的。在步骤二中，使用三次平滑样条插值的方法，首先处理中度干扰和重度干扰，再处理基线突变。对于轻度干扰区域，则直接跳过步骤二而在步骤三中进行去除。

在中度干扰和重度干扰的去除时，对每一个干扰区进行三次平滑样条插值(平滑参数的设置值逼近1，类似三次样条插值)，并用原始NIRS信号减去插值后的信号得到干扰去除后的信号。然后对于基线突变进行去除，对含有基线突变的干扰区同样进行三次平滑样条插值(平滑参数的设置值逼近0，类似线性插值)，用原始NIRS信号减去对干扰区进行插值得到的信号后，还需要在这个信号的基础上加上相对于此基线突变干扰区开始时的直流偏移量，以保证信号的连续性。

至此，已经从原始NIRS信号中去除了严重干扰的部分，在去除干扰引起的大幅度信号波动的同时，也尽可能多地保留了细节的有用信息。

具体为：

利用三次平滑样条插值对NIRS信号的干扰区域进行拟合，再用X(n)减去样条插值的拟合结果可以去除严重干扰。三次平滑样条插值函数用f来表示，通过求解下式的最小值来实现：

其中X(n)为原始Δ[HbO₂]-N信号，N是信号长度，D²f表示函数f的二阶导数，w(n)是权重向量，默认值为全1。p是平滑参数，当p值为1则是三次样条插值，p值为0则是线性插值。

在进行中度干扰和重度干扰的去除时，将p值取为0.9，对每一个中度或重度干扰区间进行样条插值，并用X(n)减去插值后的信号即可去除中度与重度干扰。

每一个基线突变干扰进行去除时，将p值取为1/(fs)³，对于一个存在基线突变的干扰区间进行样条插值，用去除中、重度干扰后的X(n)减去插值后的结果，并加上此干扰区开始时的直流偏移量，即完成一个基线突变干扰的去除。

将完成步骤二严重干扰处理后的Δ[HbO₂]-N信号记为y(n)，其结果如图13所示。

步骤三：小波双阈值处理方法。

对于NIRS信号，运动干扰的存在通常会在离散小波系数中呈现出离群值，根据这一点，本发明提出了双阈值的离散小波去噪法来处理剩下的一些运动干扰。采用两个重要参数w₁和w₂来判断离群值，其中w₁称作检测阈值，用来检测哪些小波系数应被视为运动干扰，而w₂称作处理阈值，略小于w₁，用来决定具体哪些小波系数应该被处理。而双阈值的概念即是综合利用检测阈值和处理阈值来检测小波系数中的离群点。

首先对需要进行干扰去除的信号进行离散小波分解得到每一层的小波系数，并对每一层的小波系数分别进行干扰去除。首先估计某一层小波系数的标准差，然后根据小波系数的标准差及概率参数α来确定双阈值的值。对于此层离散小波变换得到的小波系数，逐个比较其绝对值是否大于w₁，如果是，则说明该小波系数中包含了运动干扰的情况；然后在此超过检测阈值w₁的小波系数附近，再次检测大于绝对值大于w₂的小波系数，一直到其左右两侧小波系数的绝对值小于w₂为止，标记在这个过程中检测到的所有值大于w₂的小波系数，并将值其置为0，而其余的小波系数的值则保持不变。最后用处理后的所有层的小波系数去重构得到希望的NIRS信号。

为了避免边界处出现吉布斯或者突变现象，可借鉴循环平均的思想，这样可以减少由于离散小波变换带来的尖峰干扰。

具体为：

首先对信号y(n)进行离散小波分解，得到相应的小波系数。

w_jk＝w^*_jk+n_jk>

其中是有用信号的离散小波系数，n_jk是噪声的离散小波系数。服从一个正态分布，即而n_jk只占小波系数中极小的一部分。假设第j层所有小波系数w_jk的长度为L，将其进行升序排列，然后找到距离0.1587L位置点最近的小波系数的值，记为表示的是对所有小波系数的方差估计值，可以用来描述小波系数w_jk的分布。

根据公式(7)，可从α中得到分布参数u，之后即可将w₁的参考值设置为将w₂的参考值设为

α＝2[1-Φ(-u)] (7)

此处Φ指标准正态分布的分布函数。令概率参数α＝0.01，即确定相应的双阈值w₁和w₂。

按照小波双阈值处理方法中的描述，将小波系数与双阈值进行比较，标记需要处理的小波系数并将其值置为0，如图14所示。

最后用处理后的所有层的小波系数去重构得到经过干扰去除的NIRS信号y′(n)，在小波分解和重构过程利用公式(8)进行循环平均。

其中y(n)是由步骤二得到的去除严重干扰后的Δ[HbO₂]-N信号，S是移位操作，T代表离散小波的分解和重构过程，M为循环的次数且M＝16。

对于血氧信号，感兴趣的频段一般为0.003-0.04Hz，因此可对信号y′(n)进行截止频率为0.003Hz的高通滤波并最终得到经过MSBAR算法去除干扰后的较为纯净的NIRS信号Y(n)，如图15所示。