基于总体最小二乘法的近红外脑功能信号处理方法

摘要

基于总体最小二乘法的近红外脑功能信号处理方法，本发明涉及近红外脑功能信号处理方法。本发明目的是为了解决现有技术对脑灰质层中的近红外脑功能活动响应信号检测提取精度低的问题。过程为：一：获得五个不同波长的近红外光在距离检测器D相同距离R下的光密度变化量的时间信号；二：采用修正郎伯比尔定律对步骤一获得的五个不同波长的近红外光在距离检测器D相同距离R下的光密度变化量的时间信号构建方程组；三：将步骤二中的方程组改写为如下的矩阵形式；四：对三中的增广矩阵B进行奇异值分解；五：得到检测器处的氧合血红蛋白浓度变化时间信号和检测器处的还原血红蛋白浓度变化时间信号的总体最小二乘解。本发明用于信号提取领域。

说明书

技术领域

本发明涉及一种信号提取方法，具体涉及近红外脑功能信号处理方法。

背景技术

连续波近红外光谱技术可通过检测大脑皮层中的还原血红蛋白浓度和氧合血红蛋白的变化信息，提供脑功能活动过程中的大脑皮层血氧变化信息，从而用于脑功能活动响应检测。与正电子放射断层扫描、功能性磁共振成像、脑电检测等脑功能检测方法相比，近红外光谱技术具有非侵入、易实施、安全性好等优点。

在利用近红外光谱技术对脑功能活动进行检测时，需要利用修正郎伯比尔定律对光源检测器获取的光密度变化量的时间序列信号进行处理，从而得到氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号。但是，修正郎伯比尔定律中所使用的差分路径因子通常为文献中的参考值，该数值与实际测试对象的真实差分路径因子之间通常是存在误差，同时测量得到的光密度变化量时间序列信号中也存在着测量误差干扰，从而导致通过修正郎伯比尔定律直接解算得到的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号中存在着误差干扰，导致后续对脑灰质层中的近红外脑功能活动响应信号的检测提取精度低。

发明内容

本发明目的是为了解决现有技术对脑灰质层中的近红外脑功能活动响应信号检测提取精度低的问题，从而提出了基于总体最小二乘法的近红外脑功能信号处理方法。

上述的发明目的是通过以下技术方案实现的：

步骤一：在待测脑组织头皮表面放置一个由光源S与检测器D所构成的近红外探头，光源S与检测器D之间的直线距离为R，光源S发出的近红外光的波长分别为λ₁、λ₂、λ₃、λ₄和λ₅，检测器D用于获取大脑安静状态下的漫反射光强和大脑诱发激励状态下的漫反射光强，从而获得五个不同波长的近红外光在距离检测器D相同距离R下的光密度变化量的时间信号：和

其中，t为采样时刻，t＝1,2,…,N，N为正整数(此处表示t的取值范围是从1到N)；

为光源S发出近红外光的波长为λ₁时，检测器D获得的光密度变化量时间信号；

为光源S发出近红外光的波长为λ₂时，检测器D获得的光密度变化量时间信号；

为光源S发出近红外光的波长为λ₃时，检测器D获得的光密度变化量时间信号；

为光源S发出近红外光的波长为λ₄时，检测器D获得的光密度变化量时间信号；

为光源S发出近红外光的波长为λ₅时，检测器D获得的光密度变化量时间信号；

步骤二：采用修正郎伯比尔定律对步骤一获得的五个不同波长的近红外光在距离检测器D相同距离R下的光密度变化量的时间信号构建方程组，具体方程组表示为：

其中，ε_HHb(λ₁)为光源S发出近红外光的波长为λ₁时的还原血红蛋白消光系数；ε_HHb(λ₂)为光源S发出近红外光的波长为λ₂时的还原血红蛋白消光系数；ε_HHb(λ₃)为光源S发出近红外光的波长为λ₃时的还原血红蛋白消光系数；ε_HHb(λ₄)为光源S发出近红外光的波长为λ₄时的还原血红蛋白消光系数；ε_HHb(λ₅)为光源S发出近红外光的波长为λ₅时的还原血红蛋白消光系数；为光源S发出近红外光的波长为λ₁时的氧合血红蛋白消光系数；为光源S发出近红外光的波长为λ₂时的氧合血红蛋白消光系数；为光源S发出近红外光的波长为λ₃时的氧合血红蛋白消光系数；为光源S发出近红外光的波长为λ₄时的氧合血红蛋白消光系数；为光源S发出近红外光的波长为λ₅时的氧合血红蛋白消光系数；Δ[HbO₂](t)为检测器D处的氧合血红蛋白浓度变化时间信号；Δ[HHb](t)为检测器D处的还原血红蛋白浓度变化时间信号；为光源S发出近红外光的波长为λ₁时的差分路径因子；为光源S发出近红外光的波长为λ₂时的差分路径因子；为光源S发出近红外光的波长为λ₃时的差分路径因子；为光源S发出近红外光的波长为λ₄时的差分路径因子；为光源S发出近红外光的波长为λ₅时差分路径因子；

步骤三：将步骤二中的方程组改写为如下的矩阵形式：

Bz＝0

其中，

B为增广矩阵，z为待求解矩阵；

步骤四：对步骤三中的增广矩阵B进行奇异值分解，即将增广矩阵B分解，表示为：B＝UΣV^H

其中，矩阵为增广矩阵B的左奇异向量矩阵，数据u_ij为矩阵U第i行第j列对应元素，i,j＝1,2,3,4,5；矩阵为增广矩阵B的右奇异向量矩阵，数据v_kl为矩阵V第k行第l列对应元素，k,l＝1,2,3；矩阵为对角矩阵，对角元素σ₁、σ₂和σ₃为增广矩阵B的奇异值；

步骤五：利用步骤四中求取的矩阵V中第三列的对应元素，得到检测器D处的氧合血红蛋白浓度变化时间信号和检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解；分别表示为：

Δ[HbO₂]_TLS(t)为检测器D处的氧合血红蛋白浓度变化时间信号的总体最小二乘解；Δ[HHb]_TLS(t)为检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解。

发明效果

本发明提出基于总体最小二乘法的近红外脑功能信号处理方法，采用修正郎伯比尔定律求解氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号时，不考虑噪声信号干扰和分路径因子偏差影响，通过采用5组不同的近红外波长所得到的光强信号来构建超定方程组，并利用总体最小二乘法来求解该方程，总体最小二乘法同时综合考虑了测量数据误差干扰以及差分路径因子参数偏差对解算结果的影响，可有效抑制误差干扰，实现对误差干扰下的真实近红外脑功能活动信号的高精度检测，更适用于实际应用中的近红外脑功能活动信号处理分析。解决了由于修正郎伯比尔定律中使用的差分路径因子文献参考值与真实值之间存在偏差以及测量光强变化信号中存在噪声信号干扰使得利用修正修正郎伯比尔定律解算得到的脑功能活动响应信号存在误差干扰，进而导致后续对脑灰质层中的近红外脑功能活动响应信号检测提取精度低的问题。

附图说明

图1为本发明中采用的单距离近红外脑功能活动检测探头与五层脑组织模型结构示意图，其中S表示光源，D表示光源检测器，R表示光源S与检测器D之间的直线距离，L₁为头皮，L₂为颅骨，L₃为脑脊髓液，L₄为脑灰质，L₅为脑白质。

图2是本发明流程图。

具体实施方式

具体实施方式一：结合图1、图2说明本实施方式，本实施方式的基于总体最小二乘法的近红外脑功能信号处理方法，具体是按照以下步骤实现的：

步骤一：在待测脑组织头皮表面放置一个由五波长光源S与检测器D所构成的近红外探头，光源S与检测器D之间的直线距离为R，光源S发出的近红外光的波长分别为λ₁、λ₂、λ₃、λ₄和λ₅，检测器D用于获取大脑安静状态下的漫反射光强和大脑诱发激励状态下的漫反射光强，从而获得五个不同波长的近红外光在距离检测器D相同距离R下的光密度变化量的时间信号：和安静状态是指大脑在不进行思考，保持休息情况下的状态；脑诱发激励状态是指大脑在响应外界刺激进行思考动情况下的状态。

其中，t为采样时刻，t＝1,2,…,N，N为正整数(此处表示t的取值范围是从1到N)；

为光源S发出近红外光的波长为λ₁时，检测器D获得的光密度变化量时间信号；

为光源S发出近红外光的波长为λ₂时，检测器D获得的光密度变化量时间信号；

为光源S发出近红外光的波长为λ₃时，检测器D获得的光密度变化量时间信号；

为光源S发出近红外光的波长为λ₄时，检测器D获得的光密度变化量时间信号；

为光源S发出近红外光的波长为λ₅时，检测器D获得的光密度变化量时间信号；

步骤二：采用修正郎伯比尔定律对步骤一获得的五个不同波长的近红外光在距离检测器D相同距离R下的光密度变化量的时间信号构建方程组，具体方程组表示为：

其中，ε_HHb(λ₁)为光源S发出近红外光的波长为λ₁时的还原血红蛋白消光系数；ε_HHb(λ₂)为光源S发出近红外光的波长为λ₂时的还原血红蛋白消光系数；ε_HHb(λ₃)为光源S发出近红外光的波长为λ₃时的还原血红蛋白消光系数；ε_HHb(λ₄)为光源S发出近红外光的波长为λ₄时的还原血红蛋白消光系数；ε_HHb(λ₅)为光源S发出近红外光的波长为λ₅时的还原血红蛋白消光系数；为光源S发出近红外光的波长为λ₁时的氧合血红蛋白消光系数；为光源S发出近红外光的波长为λ₂时的氧合血红蛋白消光系数；为光源S发出近红外光的波长为λ₃时的氧合血红蛋白消光系数；为光源S发出近红外光的波长为λ₄时的氧合血红蛋白消光系数；为光源S发出近红外光的波长为λ₅时的氧合血红蛋白消光系数；Δ[HbO₂](t)为检测器D处的氧合血红蛋白浓度变化时间信号；Δ[HHb](t)为检测器D处的还原血红蛋白浓度变化时间信号；为光源S发出近红外光的波长为λ₁时的差分路径因子；为光源S发出近红外光的波长为λ₂时的差分路径因子；为光源S发出近红外光的波长为λ₃时的差分路径因子；为光源S发出近红外光的波长为λ₄时的差分路径因子；为光源S发出近红外光的波长为λ₅时差分路径因子；

步骤三：将步骤二中的方程组改写为如下的矩阵形式：

Bz＝0

其中，

B为增广矩阵，z为待求解矩阵；

步骤四：对步骤三中的增广矩阵B进行奇异值分解，即将增广矩阵B分解，表示为：B＝UΣV^H

其中，矩阵为增广矩阵B的左奇异向量矩阵，数据u_ij为矩阵U第i行第j列对应元素，i,j＝1,2,3,4,5；矩阵为增广矩阵B的右奇异向量矩阵，数据v_kl为矩阵V第k行第l列对应元素，k,l＝1,2,3；矩阵为对角矩阵，对角元素σ₁、σ₂和σ₃为增广矩阵B的奇异值；

步骤五：利用步骤四中求取的矩阵V中第三列的对应元素，得到检测器D处的氧合血红蛋白浓度变化时间信号和检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解；分别表示为：

式中，Δ[HbO₂]_TLS(t)为检测器D处的氧合血红蛋白浓度变化时间信号的总体最小二乘解；Δ[HHb]_TLS(t)为检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中5mm≤R≤40mm；

其中，R表示光源S与检测器D之间的直线距离。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一至二之一不同的是：所述R为5mm。

其它步骤及参数与具体实施方式一至二之一相同。

具体实施方式四：本实施方式与具体实施方式一至二之一不同的是：所述R为40mm。

其它步骤及参数与具体实施方式一至二之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤一中光源S发出的近红外光的波长λ₁为670nm、λ₂为770nm、λ₃为810nm、λ₄为850nm、λ₅为950nm。

其它步骤及参数与具体实施方式一至四之一相同。