[发明专利]基于差分路径因子估计的近红外脑功能信号处理方法

摘要

基于差分路径因子估计的近红外脑功能信号处理方法，本发明涉及近红外脑功能信号处理方法。本发明目的是为了解决现有技术中修正郎伯比尔定律所使用的差分路径因子参考值与实际测量对象的真实差分路径因子存在较大的差异，同时光源检测器所采集到的光密度变化量的时间序列信号中也存在着测量误差干扰，导致对连续波近红外脑功能活动响应信号测量提取精度低的问题。获得不同波长的近红外光在距离检测器相同距离下的光密度变化量的时间信号；采用修正郎伯比尔定律对信号构建方程；将方程组改写为矩阵形式；对增广矩阵进行奇异值分解；得到检测器处的氧合血红蛋白和还原血红蛋白浓度变化时间信号的总体最小二乘解。本发明用于脑功能信号领域。

主权项

基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：所述方法具体过程为：步骤一：在待测脑组织头皮表面放置一个由五波长光源S与检测器D所构成的近红外探头，五波长光源S与检测器D之间的直线距离为R，五波长光源S发出的近红外光的波长分别为λ1、λ2、λ3、λ4和λ5，检测器D用于获取大脑安静状态下的漫反射光强和大脑诱发激励状态下的漫反射光强，从而获得五个不同波长的近红外光在距检测器D相同距离R下的光密度变化量的时间信号：和其中，t为采样时刻，t＝1,2,…,N，N为正整数；为五波长光源S发出近红外光的波长为λ1时，检测器D获得的光密度变化量时间信号；为五波长光源S发出近红外光的波长为λ2时，检测器D获得的光密度变化量时间信号；为五波长光源S发出近红外光的波长为λ3时，检测器D获得的光密度变化量时间信号；为五波长光源S发出近红外光的波长为λ4时，检测器D获得的光密度变化量时间信号；为五波长光源S发出近红外光的波长为λ5时，检测器D获得的光密度变化量时间信号；步骤二：对步骤一获得的五个不同波长下的近红外光在距检测器D相同距离R下的光密度变化量时间信号中每两组采用修正郎伯比尔定律构建一个方程组，总计构建十个方程组，具体的方程组表示为：ΔODλi(t)ΔODλj(t)=ϵHbO2(λi)DPF(λi)RϵHHb(λi)DPF(λi)RϵHbO2(λj)DPF(λj)RϵHHb(λj)DPF(λj)RΔ[HbO2]λi,λj(t)Δ[HHb]λi,λj(t)]]>其中，i为数字标号，取值范围是i＝1,2,3,4,5；j为数字标号，取值范围是j＝1,2,3,4,5；λi或λj为五波长光源S发出的近红外光的波长λ1、λ2、λ3、λ4和λ5中的某一个，同时方程组中的数字标号i和j不能取相同值；和表示五波长光源S发出近红外光的波长为λi或λj时，检测器D获得的光密度变化量时间信号；εHHb(λi)或εHHb(λj)为五波长光源S发出近红外光的波长为λi或λj时的还原血红蛋白消光系数；或为五波长光源S发出近红外光的波长为λi或λj时的氧合血红蛋白消光系数；DPF(λi)或DPF(λj)为五波长光源S发出近红外光的波长为λi或λj时的差分路径因子；为采用近红外光波长组合(λi,λj)时，检测器D处的氧合血红蛋白浓度变化时间信号；为采用近红外光波长组合(λi,λj)时，检测器D处的还原血红蛋白浓度变化时间信号；步骤三：对步骤二中得到的十个方程组分别进行求解，得到十组采用不同近红外波长的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号分别表示如下：Δ[HbO2]λi,λj(t)=ϵHHb(λj)ΔODλi(t)DPF(λi)·R-ϵHHb(λi)ΔODλj(t)DPF(λj)·RϵHbO2(λi)ϵHHb(λj)-ϵHHb(λi)ϵHbO2(λj)]]>Δ[HHb]λi,λj(t)=ϵHbO2(λj)ΔODλi(t)DPF(λi)·R-ϵHbO2(λi)ΔODλj(t)DPF(λj)·RϵHHb(λi)ϵHbO2(λj)-ϵHbO2(λi)ϵHHb(λj)]]>步骤四：生成初始化的粒子群，其中第m个粒子的初始化位置向量表示如下：xm(0)＝(xm1(0),xm2(0),xm3(0),xm4(0),xm5(0))其中，m为每个粒子的标号，m＝1,2,…,M，M为正整数，代表步骤四中生成的初始化粒子群的粒子数量；xm1(0)为初始化位置向量中的第一个元素，代表初始化生成的DPF(λ1)数据；xm2(0)为初始化位置向量中的第二个元素，代表初始化生成的DPF(λ2)数据，xm3(0)为初始化位置向量中的第三个元素，代表初始化生成的DPF(λ3)数据，xm4(0)为初始化位置向量中的第四个元素，代表初始化生成的DPF(λ4)数据，xm5(0)为初始化位置向量中的第五个元素，代表初始化生成的DPF(λ5)数据；第m个粒子的初始化速度向量表示如下：vm(0)＝(vm1(0),vm2(0),vm3(0),vm4(0),vm5(0))其中，vm(0)为第m个粒子的初始化速度向量；vm1(0)为初始化速度向量中的第一个元素；vm2(0)为初始化速度向量中的第二个元素；vm3(0)为初始化速度向量中的第三个元素；vm4(0)为初始化速度向量中的第四个元素；vm5(0)为初始化速度向量中的第五个元素；步骤五：计算步骤四中生成的初始化粒子群中每个粒子的适应度值；过程为：首先将步骤四中生成的所有粒子的初始化位置向量分别带入到步骤三中得到的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号表达式之中，则每个粒子得到十组采用不同近红外波长的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号；对其中的第m个粒子，利用如下的公式来计算采用不同近红外波长组合下所得到氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值Dij,pqm=Σt=1N(Δ[HbO2]λi,λjm(t)-Δ[HbO2]λp,λqm(t))2Σt=1N(Δ[HbO2]λi,λjm(t))2+Σt=1N(Δ[HHb]λi,λjm(t)-Δ[HHb]λp,λqm(t))2Σt=1N(Δ[HHb]λi,λjm(t))2]]>其中p为数字标号，取值范围是p＝1,2,3,4,5；q为数字标号，取值范围是q＝1,2,3,4,5；λp或λq为五波长光源S发出的近红外光的波长λ1、λ2、λ3、λ4和λ5中的某一个，同时数字标号p和q不能取相同值；表示为当采用第m个粒子产生的位置向量时，波长组合(λi,λj)计算得到氧合血红蛋白浓度变化时间信号；表示为当采用第m个粒子产生的位置向量时，波长组合(λp,λq)计算得到氧合血红蛋白浓度变化时间信号；表示为当采用第m个粒子产生的位置向量时，波长组合(λi,λj)计算得到还原血红蛋白浓度变化时间信号；表示为当采用第m个粒子产生的位置向量时，波长组合(λp,λq)计算得到还原血红蛋白浓度变化时间信号；表示为当采用第m个粒子产生的位置向量时，采用波长组合(λi,λj)和(λp,λq)所得到的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值；然后，利用如下的适应度值计算公式来计算第m个粒子的适应度值：Om=Σi,j,p,q=15Dij,pqm]]>步骤六：计算每个粒子的个体极值点和粒子群的全局极值点；过程为：对于生成的初始化粒子群，每个粒子的个体极值点为当前的位置向量，对于第m个粒子，其个体极值点表示为：pbestm(0)＝(pbestm1(0),pbestm2(0),pbestm3(0),pbestm4(0),pbestm5(0))其中pbestm(0)为第m个粒子的个体极值点，pbestm1(0)为个体极值点中的第一个元素，pbestm2(0)为个体极值点中的第二个元素，pbestm3(0)为个体极值点中的第三个元素，pbestm4(0)为个体极值点中的第四个元素，pbestm5(0)为个体极值点中的第五个元素，此时pbestm(0)＝xm(0)；粒子群的全局极值点是按照步骤五计算每个粒子的适应度值之后，具有最小值的适应度值的粒子的位置向量为整个粒子群中的全局极值点，表示为：gbest(0)＝(gbest1(0),gbest2(0),gbest3(0),gbest4(0),gbest5(0))其中gbest1(0)为全局极值点中的第一个元素，gbest2(0)为全局极值点中的第二个元素，gbest3(0)为全局极值点中的第三个元素，gbest4(0)为全局极值点中的第四个元素，gbest5(0)为全局极值点中的第五个元素；步骤七：步骤四中初始化生成的粒子群中的各个粒子的位置向量与速度向量的迭代更新公式分别表示如下：vmd(k)＝wvmd(k‑1)+c1rand1(k‑1)(pbestmd(k‑1)‑xmd(k‑1))+c2rand2(k‑1)(gbestd(k‑1)‑xmd(k‑1))xmd(k)＝xmd(k‑1)+vmd(k)其中，k表示粒子群的迭代更新次数，k＝1,2,…,K，K为预先设置的迭代总次数，K为正整数；d为位置向量和速度向量中的元素标号，d＝1,2,3,4,5；vmd(k)表示粒子群中第m个粒子在第k次迭代时速度向量vm(k)＝(vm1(k),vm2(k),vm3(k),vm4(k),vm5(k))中的第d个元素，d＝1,2,3,4,5；xmd(k)表示粒子群中第m个粒子在第k次迭代时位置向量xm(k)＝(xm1(k),xm2(k),xm3(k),xm4(k),xm5(k))中的第d个元素，d＝1,2,3,4,5；pbestmd(k‑1)表示粒子群中第m个粒子在第k‑1次迭代时个体极值点中的第d个元素；gbestd(k‑1)表示粒子群在第k‑1次迭代时全局极值点中的第d个元素；w为惯性因子；c1和c2是加速因子；当k＝0时，rand1(0)＝rand2(0)＝0，当k≠0时，rand1(k)和rand2(k)为随机数，rand1(k)取值范围为0＜rand1(k)＜1，rand2(k)取值范围为0＜rand2(k)＜1；步骤八：将步骤七中得到的第k次迭代更新后的各粒子的位置向量按照步骤五中的公式分别计算各粒子的适应度值；对第m个粒子，如果第k次迭代得到的位置向量xm(k)计算得到的适应度值小于前一次迭代得到的个体极值点计算得到的适应度值，则将位置向量xm(k)设置为第k次迭代得到个体极值点pbestm(k)；否则维持个体极值点不变，即pbestm(k)＝pbestm(k‑1)；如果所有粒子在第k次迭代得到个体极值点计算得到适应度值中的最小值小于前一次迭代得到全局极值点计算得到适应度值，则将该粒子的位置向量记为全局极值点，并更新gbest(k)；否则维持全局极值点不变，即gbest(k)＝gbest(k‑1)；步骤九：判断迭代次数k是否达到K，如果没有达到则重复步骤七和步骤八，不断更新各粒子的个体极值点和整个粒子群的全局极值点；当迭代次数k达到K时，则得到的全局极值点表示为：gbest(K)＝(gbest1(K),gbest2(K),gbest3(K),gbest4(K),gbest5(K))其中gbest1(K)为近红外光波长的差分路径因子DPF(λ1)的最优估计值，gbest2(K)为近红外光波长的差分路径因子DPF(λ2)的最优估计值，gbest3(K)为近红外光波长的差分路径因子DPF(λ3)的最优估计值，gbest4(K)为近红外光波长的差分路径因子DPF(λ4)的最优估计值，gbest5(K)为近红外光波长的差分路径因子DPF(λ5)的最优估计值；步骤十，利用步骤九中得到的各个波长下的差分路径因子最优估计值构建如下的方程组来求取氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号；步骤十一：将步骤十中的方程组改写为如下的矩阵形式：Bz＝0其中，B为增广矩阵，z为待求解矩阵；步骤十二：对步骤十一中的增广矩阵B进行奇异值分解，即将增广矩阵B分解，表示为：B＝UΣVH其中，H为共轭转置；矩阵为增广矩阵B的左奇异向量矩阵，数据uij为矩阵U第i行第j列对应元素，i,j＝1,2,3,4,5；矩阵为增广矩阵B的右奇异向量矩阵，数据vkl为矩阵V第k行第l列对应元素，k,l＝1,2,3；矩阵为对角矩阵，对角元素σ1、σ2和σ3为增广矩阵B的奇异值；步骤十三：利用步骤十二中求取的矩阵V中第三列的对应元素，得到检测器D处的氧合血红蛋白浓度变化时间信号和检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解。