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

权利要求书

1.基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：所述方法具体过程为：

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

其中，t为采样时刻，t＝1,2,…,N，N为正整数；

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

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

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

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

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

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

其中，i为数字标号，取值范围是i＝1,2,3,4,5；j为数字标号，取值范围是j＝1,2,3,4,5；λ_i或λ_j为五波长光源S发出的近红外光的波长λ₁、λ₂、λ₃、λ₄和λ₅中的某一个，同时方程组中的数字标号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处的还原血红蛋白浓度变化时间信号；

步骤三：对步骤二中得到的十个方程组分别进行求解，得到十组采用不同近红外波长的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号分别表示如下：

步骤四：生成初始化的粒子群，其中第m个粒子的初始化位置向量表示如下：

x_m(0)＝(x_m1(0),x_m2(0),x_m3(0),x_m4(0),x_m5(0))

其中，m为每个粒子的标号，m＝1,2,…,M，M为正整数，代表步骤四中生成的初始化粒子群的粒子数量；x_m1(0)为初始化位置向量中的第一个元素，代表初始化生成的DPF(λ₁)数据；x_m2(0)为初始化位置向量中的第二个元素，代表初始化生成的DPF(λ₂)数据，x_m3(0)为初始化位置向量中的第三个元素，代表初始化生成的DPF(λ₃)数据，x_m4(0)为初始化位置向量中的第四个元素，代表初始化生成的DPF(λ₄)数据，x_m5(0)为初始化位置向量中的第五个元素，代表初始化生成的DPF(λ₅)数据；

第m个粒子的初始化速度向量表示如下：

v_m(0)＝(v_m1(0),v_m2(0),v_m3(0),v_m4(0),v_m5(0))

其中，v_m(0)为第m个粒子的初始化速度向量；v_m1(0)为初始化速度向量中的第一个元素；v_m2(0)为初始化速度向量中的第二个元素；v_m3(0)为初始化速度向量中的第三个元素；v_m4(0)为初始化速度向量中的第四个元素；v_m5(0)为初始化速度向量中的第五个元素；

步骤五：计算步骤四中生成的初始化粒子群中每个粒子的适应度值；过程为：

首先将步骤四中生成的所有粒子的初始化位置向量分别带入到步骤三中得到的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号表达式之中，则每个粒子得到十组采用不同近红外波长的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号；

对其中的第m个粒子，利用如下的公式来计算采用不同近红外波长组合下所得到氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值

其中p为数字标号，取值范围是p＝1,2,3,4,5；q为数字标号，取值范围是q＝1,2,3,4,5；λ_p或λ_q为五波长光源S发出的近红外光的波长λ₁、λ₂、λ₃、λ₄和λ₅中的某一个，同时数字标号p和q不能取相同值；表示为当采用第m个粒子产生的位置向量时，波长组合(λ_i,λ_j)计算得到氧合血红蛋白浓度变化时间信号；表示为当采用第m个粒子产生的位置向量时，波长组合(λ_p,λ_q)计算得到氧合血红蛋白浓度变化时间信号；表示为当采用第m个粒子产生的位置向量时，波长组合(λ_i,λ_j)计算得到还原血红蛋白浓度变化时间信号；表示为当采用第m个粒子产生的位置向量时，波长组合(λ_p,λ_q)计算得到还原血红蛋白浓度变化时间信号；表示为当采用第m个粒子产生的位置向量时，采用波长组合(λ_i,λ_j)和(λ_p,λ_q)所得到的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值；

然后，利用如下的适应度值计算公式来计算第m个粒子的适应度值：

步骤六：计算每个粒子的个体极值点和粒子群的全局极值点；过程为：

对于生成的初始化粒子群，每个粒子的个体极值点为当前的位置向量，对于第m个粒子，其个体极值点表示为：

pbest_m(0)＝(pbest_m1(0),pbest_m2(0),pbest_m3(0),pbest_m4(0),pbest_m5(0))

其中pbest_m(0)为第m个粒子的个体极值点，pbest_m1(0)为个体极值点中的第一个元素，pbest_m2(0)为个体极值点中的第二个元素，pbest_m3(0)为个体极值点中的第三个元素，pbest_m4(0)为个体极值点中的第四个元素，pbest_m5(0)为个体极值点中的第五个元素，此时pbest_m(0)＝x_m(0)；

粒子群的全局极值点是按照步骤五计算每个粒子的适应度值之后，具有最小值的适应度值的粒子的位置向量为整个粒子群中的全局极值点，表示为：

gbest(0)＝(gbest₁(0),gbest₂(0),gbest₃(0),gbest₄(0),gbest₅(0))

其中gbest₁(0)为全局极值点中的第一个元素，gbest₂(0)为全局极值点中的第二个元素，gbest₃(0)为全局极值点中的第三个元素，gbest₄(0)为全局极值点中的第四个元素，gbest₅(0)为全局极值点中的第五个元素；

步骤七：步骤四中初始化生成的粒子群中的各个粒子的位置向量与速度向量的迭代更新公式分别表示如下：

v_md(k)＝wv_md(k-1)+c₁rand₁(k-1)(pbest_md(k-1)-x_md(k-1))+c₂rand₂(k-1)(gbest_d(k-1)-x_md(k-1))

x_md(k)＝x_md(k-1)+v_md(k)

其中，k表示粒子群的迭代更新次数，k＝1,2,…,K，K为预先设置的迭代总次数，K为正整数；d为位置向量和速度向量中的元素标号，d＝1,2,3,4,5；v_md(k)表示粒子群中第m个粒子在第k次迭代时速度向量v_m(k)＝(v_m1(k),v_m2(k),v_m3(k),v_m4(k),v_m5(k))中的第d个元素，d＝1,2,3,4,5；x_md(k)表示粒子群中第m个粒子在第k次迭代时位置向量x_m(k)＝(x_m1(k),x_m2(k),x_m3(k),x_m4(k),x_m5(k))中的第d个元素，d＝1,2,3,4,5；pbest_md(k-1)表示粒子群中第m个粒子在第k-1次迭代时个体极值点中的第d个元素；gbest_d(k-1)表示粒子群在第k-1次迭代时全局极值点中的第d个元素；w为惯性因子；c₁和c₂是加速因子；

当k＝0时，rand₁(0)＝rand₂(0)＝0，当k≠0时，rand₁(k)和rand₂(k)为随机数，rand₁(k)取值范围为0＜rand₁(k)＜1，rand₂(k)取值范围为0＜rand₂(k)＜1；

步骤八：将步骤七中得到的第k次迭代更新后的各粒子的位置向量按照步骤五中的公式分别计算各粒子的适应度值；

对第m个粒子，如果第k次迭代得到的位置向量x_m(k)计算得到的适应度值小于前一次迭代得到的个体极值点计算得到的适应度值，则将位置向量x_m(k)设置为第k次迭代得到个体极值点pbest_m(k)；否则维持个体极值点不变，即pbest_m(k)＝pbest_m(k-1)；

如果所有粒子在第k次迭代得到个体极值点计算得到适应度值中的最小值小于前一次迭代得到全局极值点计算得到适应度值，则将该粒子的位置向量记为全局极值点，并更新gbest(k)；否则维持全局极值点不变，即gbest(k)＝gbest(k-1)；

步骤九：判断迭代次数k是否达到K，如果没有达到则重复步骤七和步骤八，不断更新各粒子的个体极值点和整个粒子群的全局极值点；当迭代次数k达到K时，则得到的全局极值点表示为：

gbest(K)＝(gbest₁(K),gbest₂(K),gbest₃(K),gbest₄(K),gbest₅(K))

其中gbest₁(K)为近红外光波长的差分路径因子DPF(λ₁)的最优估计值，gbest₂(K)为近红外光波长的差分路径因子DPF(λ₂)的最优估计值，gbest₃(K)为近红外光波长的差分路径因子DPF(λ₃)的最优估计值，gbest₄(K)为近红外光波长的差分路径因子DPF(λ₄)的最优估计值，gbest₅(K)为近红外光波长的差分路径因子DPF(λ₅)的最优估计值；

步骤十，利用步骤九中得到的各个波长下的差分路径因子最优估计值构建如下的方程组来求取氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号；

具体方程组可表示如下：

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

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

Bz＝0

其中，

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

步骤十二：对步骤十一中的增广矩阵B进行奇异值分解，即将增广矩阵B分解，表示为：

B＝UΣV^H

其中，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处的还原血红蛋白浓度变化时间信号的总体最小二乘解。