CN107928631B - 基于差分路径因子估计的近红外脑功能信号处理方法 - Google Patents

基于差分路径因子估计的近红外脑功能信号处理方法 Download PDF

Info

Publication number
CN107928631B
CN107928631B CN201711392782.8A CN201711392782A CN107928631B CN 107928631 B CN107928631 B CN 107928631B CN 201711392782 A CN201711392782 A CN 201711392782A CN 107928631 B CN107928631 B CN 107928631B
Authority
CN
China
Prior art keywords
wavelength
particle
lambda
gbest
near infrared
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201711392782.8A
Other languages
English (en)
Other versions
CN107928631A (zh
Inventor
刘昕
张瞫
张岩
刘丹
孙金玮
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN201711392782.8A priority Critical patent/CN107928631B/zh
Publication of CN107928631A publication Critical patent/CN107928631A/zh
Application granted granted Critical
Publication of CN107928631B publication Critical patent/CN107928631B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0075—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence by spectroscopy, i.e. measuring spectra, e.g. Raman spectroscopy, infrared absorption spectroscopy
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0033—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
    • A61B5/004—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for image acquisition of a particular organ or body part
    • A61B5/0042—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for image acquisition of a particular organ or body part for the brain
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/145—Measuring characteristics of blood in vivo, e.g. gas concentration, pH value; Measuring characteristics of body fluids or tissues, e.g. interstitial fluid, cerebral tissue
    • A61B5/14546—Measuring characteristics of blood in vivo, e.g. gas concentration, pH value; Measuring characteristics of body fluids or tissues, e.g. interstitial fluid, cerebral tissue for measuring analytes not otherwise provided for, e.g. ions, cytochromes
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/145—Measuring characteristics of blood in vivo, e.g. gas concentration, pH value; Measuring characteristics of body fluids or tissues, e.g. interstitial fluid, cerebral tissue
    • A61B5/1455—Measuring characteristics of blood in vivo, e.g. gas concentration, pH value; Measuring characteristics of body fluids or tissues, e.g. interstitial fluid, cerebral tissue using optical sensors, e.g. spectral photometrical oximeters
    • A61B5/14551—Measuring characteristics of blood in vivo, e.g. gas concentration, pH value; Measuring characteristics of body fluids or tissues, e.g. interstitial fluid, cerebral tissue using optical sensors, e.g. spectral photometrical oximeters for measuring blood gases
    • A61B5/14553—Measuring characteristics of blood in vivo, e.g. gas concentration, pH value; Measuring characteristics of body fluids or tissues, e.g. interstitial fluid, cerebral tissue using optical sensors, e.g. spectral photometrical oximeters for measuring blood gases specially adapted for cerebral tissue
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Animal Behavior & Ethology (AREA)
  • Signal Processing (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Psychiatry (AREA)
  • Neurology (AREA)
  • Physiology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Optics & Photonics (AREA)
  • Medicines That Contain Protein Lipid Enzymes And Other Medicines (AREA)
  • Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

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

Description

基于差分路径因子估计的近红外脑功能信号处理方法
技术领域
本发明涉及近红外脑功能信号处理方法。
背景技术
连续波近红外光谱技术可检测脑组织中的氧合血红蛋白及还原血红蛋白浓度的变化信息，提供脑功能活动过程中的脑组织血氧变化信息，用于脑功能活动检测。与传统脑功能检测方法如功能性磁共振成像、正电子放射断层扫描等相比，基于近红外光谱技术的脑功能检测方法具有低成本、易实施、非侵入、安全性好等诸多优点。
当采用连续波近红外光谱检测技术对脑功能活动进行测量时，需要使用修正郎伯比尔定律对脑部安装的光源检测器所采集到的光密度变化量的时间序列信号进行信号处理来获取还原血红蛋白浓度变化时间信号和氧合血红蛋白浓度变化时间信号。然而，在修正郎伯比尔定律中需要使用到的差分路径因子一般为文献中的参考值，由于不同的实际测量对象存在较大的个体差异如不同测量对象的脑部头皮、颅骨、脑脊髓液、脑灰质、脑白质等各层脑组织的厚度因人而异，因而导致文献中的差分路径因子参考值与实际测量对象的真实差分路径因子之间通常会存在较大的差异，同时光源检测器所采集到的光密度变化量的时间序列信号中也存在着测量误差干扰，从而导致了利用修正郎伯比尔定律处理得到的还原血红蛋白浓度变化时间信号和氧合血红蛋白浓度变化时间信号中存在误差干扰，进而导致后续对连续波近红外脑功能活动响应信号测量提取精度降低。
发明内容
本发明目的是为了解决现有技术中修正郎伯比尔定律所使用的差分路径因子参考值与实际测量对象的真实差分路径因子存在较大的差异，同时光源检测器所采集到的光密度变化量的时间序列信号中也存在着测量误差干扰，导致对连续波近红外脑功能活动响应信号测量提取精度低的问题，从而提出了基于差分路径因子估计的近红外脑功能信号处理方法。
上述的发明目的是通过以下技术方案实现的：
步骤一：在待测脑组织头皮表面放置一个由五波长光源S与检测器D所构成的近红外探头，五波长光源S与检测器D之间的直线距离为R，五波长光源S发出的近红外光的波长分别为λ1、λ2、λ3、λ4和λ5，检测器D用于获取大脑安静状态下的漫反射光强和大脑诱发激励状态下的漫反射光强，从而获得五个不同波长的近红外光在距检测器D相同距离R下的光密度变化量的时间信号：
Figure BDA0001517884240000021
和
Figure BDA0001517884240000022
其中，t为采样时刻，t＝1,2,…,N，N为正整数(此处表示t的取值范围是从1到N)；
Figure BDA0001517884240000023
为五波长光源S发出近红外光的波长为λ1时，检测器D获得的光密度变化量时间信号；
Figure BDA0001517884240000024
为五波长光源S发出近红外光的波长为λ2时，检测器D获得的光密度变化量时间信号；
Figure BDA0001517884240000025
为五波长光源S发出近红外光的波长为λ3时，检测器D获得的光密度变化量时间信号；
Figure BDA0001517884240000026
为五波长光源S发出近红外光的波长为λ4时，检测器D获得的光密度变化量时间信号；
Figure BDA0001517884240000027
为五波长光源S发出近红外光的波长为λ5时，检测器D获得的光密度变化量时间信号；
步骤二：对步骤一获得的五个不同波长下的近红外光在距检测器D相同距离R下的光密度变化量时间信号中每两组采用修正郎伯比尔定律构建一个方程组，总计构建十个方程组，具体的方程组表示为：
Figure BDA0001517884240000028
此公式中
Figure BDA0001517884240000029
和
Figure BDA00015178842400000210
为未知量，需要求解。εHHb(λi)、εHHb(λj)、
Figure BDA00015178842400000211
为已知参数；DPF(λi)、DPF(λj)也为已知参数(可以使文献中直接查找到的λi和λj波长下的差分路径因子；也可以是后面粒子群算法中不断迭代产生的差分路径因子值)。
其中，i为数字标号，取值范围是i＝1,2,3,4,5；j为数字标号，取值范围是j＝1,2,3,4,5；λi或λj为五波长光源S发出的近红外光的波长λ1、λ2、λ3、λ4和λ5中的某一个，同时方程组中的数字标号i和j不能取相同值；
Figure BDA00015178842400000212
和
Figure BDA00015178842400000213
表示五波长光源S发出近红外光的波长为λi或λj时，检测器D获得的光密度变化量时间信号；εHHb(λi)或εHHb(λj)为五波长光源S发出近红外光的波长为λi或λj时的还原血红蛋白消光系数；
Figure BDA0001517884240000031
或
Figure BDA0001517884240000032
为五波长光源S发出近红外光的波长为λi或λj时的氧合血红蛋白消光系数；DPF(λi)或DPF(λj)为五波长光源S发出近红外光的波长为λi或λj时的差分路径因子；
Figure BDA0001517884240000033
为采用近红外光波长组合(λi,λj)时，检测器D处的氧合血红蛋白浓度变化时间信号；
Figure BDA0001517884240000034
为采用近红外光波长组合(λi,λj)时，检测器D处的还原血红蛋白浓度变化时间信号。
步骤三：对步骤二中得到的十个方程组分别进行求解，得到十组采用不同近红外波长的氧合血红蛋白浓度变化时间信号
Figure BDA0001517884240000035
和还原血红蛋白浓度变化时间信号
Figure BDA0001517884240000036
分别表示如下：
Figure BDA0001517884240000037
Figure BDA0001517884240000038
步骤四：生成初始化的粒子群，其中第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个粒子，利用如下的公式来计算采用不同近红外波长组合下所得到氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值
Figure BDA0001517884240000041
Figure BDA0001517884240000042
其中p为数字标号，取值范围是p＝1,2,3,4,5；q为数字标号，取值范围是q＝1,2,3,4,5；λp或λq为五波长光源S发出的近红外光的波长λ1、λ2、λ3、λ4和λ5中的某一个，同时数字标号p和q不能取相同值；
Figure BDA0001517884240000043
表示为当采用第m个粒子产生的位置向量时，波长组合(λi,λj)计算得到氧合血红蛋白浓度变化时间信号；
Figure BDA0001517884240000044
表示为当采用第m个粒子产生的位置向量时，波长组合(λp,λq)计算得到氧合血红蛋白浓度变化时间信号；
Figure BDA0001517884240000045
表示为当采用第m个粒子产生的位置向量时，波长组合(λi,λj)计算得到还原血红蛋白浓度变化时间信号；
Figure BDA0001517884240000046
表示为当采用第m个粒子产生的位置向量时，波长组合(λp,λq)计算得到还原血红蛋白浓度变化时间信号；
Figure BDA0001517884240000047
表示为当采用第m个粒子产生的位置向量时，采用波长组合(λi,λj)和(λp,λq)所得到的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值；
然后，利用如下的适应度值计算公式来计算第m个粒子的适应度值：
Figure BDA0001517884240000048
步骤六：计算每个粒子的个体极值点和粒子群的全局极值点；过程为：
对于生成的初始化粒子群，每个粒子的个体极值点为当前的位置向量，对于第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次迭代更新后的各粒子的位置向量按照步骤五中的公式分别计算各粒子的适应度值；
从步骤四开始根据位置向量得到k次迭代的DPF(λi)，在根据DPF(λi)得到步骤三中的
Figure BDA0001517884240000061
和
Figure BDA0001517884240000062
然后根据步骤五中第一个公式得到
Figure BDA0001517884240000063
最终得到适应度值；
对第m个粒子，如果第k次迭代得到的位置向量xm(k)计算得到的适应度值小于前一次迭代得到的个体极值点计算得到的适应度值，则将位置向量xm(k)设置为第k次迭代得到个体极值点pbestm(k)；否则维持个体极值点不变，即pbestm(k)＝pbestm(k-1)；
如果所有粒子在第k次迭代得到个体极值点计算得到适应度值中的最小值小于前一次迭代得到全局极值点计算得到适应度值，则将该粒子(第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
其中，
Figure BDA0001517884240000071
B为增广矩阵，z为待求解矩阵；
步骤十二：对步骤十一中的增广矩阵B进行奇异值分解，即将增广矩阵B分解，表示为：
B＝UΣVH
其中，H为共轭转置；矩阵
Figure BDA0001517884240000072
为增广矩阵B的左奇异向量矩阵，数据uij为矩阵U第i行第j列对应元素，i,j＝1,2,3,4,5；矩阵
Figure BDA0001517884240000073
为增广矩阵B的右奇异向量矩阵，数据vkl为矩阵V第k行第l列对应元素，k,l＝1,2,3；矩阵
Figure BDA0001517884240000074
为对角矩阵，对角元素σ1、σ2和σ3为增广矩阵B的奇异值；
步骤十三：利用步骤十二中求取的矩阵V中第三列的对应元素，得到检测器D处的氧合血红蛋白浓度变化时间信号和检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解。
发明效果
本发明首先通过采用五波长近红外光源构建多组不同波长组合下的基于修正郎伯比尔定律的求解方程组，从而得到多组不同波长组合下的近红外波长的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号表达式，当这五个不同波长的差分路径因子均采用真实值时，不同波长组合下的基于修正郎伯比尔定律的方程组所得到的近红外波长的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号基本相同，此时不同近红外波长组合下所得到氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值之和应为最小，将粒子群算法用于五个不同波长的差分路径因子真实值的迭代优化估计，可实现对不同波长的差分路径因子真实值的最优化估计，从而改善因所采用的差分路径因子文献值与实际真实值差异导致的近红外脑功能活动响应信号检测提取精度降低的问题，实现对真实脑功能活动信号的高质量检测。
粒子群算法是一种模拟鸟群觅食的仿生算法，类似遗传算法，粒子群算法也是一种基于群智能迭代的优化算法，但粒子群算法中不包含变异算子与交叉算子，因此粒子群算法的优化速度较快，同时粒子群算法还具有对种群数量大小不敏感的优点，即种群数量减小时优化迭代性能下降不大，因此可有效降低算法的计算复杂度，使其适用于近红外脑功能信号处理中的差分路径因子真实值最优化估计。
其次，针对差分路径因子最优化估计值与真实值之间依然存在较小偏差以及测量光强数据中存在着测量误差干扰的问题，通过采用五组不同的近红外波长所得到的光强信号来构建超定方程组，并利用总体最小二乘法来求解该方程，总体最小二乘法同时综合考虑了测量数据误差干扰以及差分路径因子优化估计参数偏差对解算结果的影响，可有效抑制误差干扰，实现对误差干扰下的真实近红外脑功能活动信号的高精度检测，更适用于实际应用中的近红外脑功能活动信号处理分析。
附图说明
图1为本发明中所采用的单距离近红外脑功能活动检测探头与五层脑组织模型结构示意图，其中S表示五波长光源，D表示光源检测器，R表示五波长光源S与近端检测器D之间的直线距离，L1为头皮，L2为颅骨，L3为脑脊髓液，L4为脑灰质，L5为脑白质；
图2为本发明流程图。
具体实施方式
具体实施方式一：结合图1、图2说明本实施方式，本实施方式的基于差分路径因子估计的近红外脑功能信号处理方法，具体是按照以下步骤实现的：
步骤一：在待测脑组织头皮表面放置一个由五波长光源S与检测器D所构成的近红外探头，五波长光源S与检测器D之间的直线距离为R，五波长光源S发出的近红外光的波长分别为λ1、λ2、λ3、λ4和λ5，检测器D用于获取大脑安静状态下的漫反射光强和大脑诱发激励状态下的漫反射光强，从而获得五个不同波长的近红外光在距检测器D相同距离R下的光密度变化量的时间信号：
Figure BDA0001517884240000091
和
Figure BDA0001517884240000092
其中，t为采样时刻，t＝1,2,…,N，N为正整数(此处表示t的取值范围是从1到N)；
Figure BDA0001517884240000093
为五波长光源S发出近红外光的波长为λ1时，检测器D获得的光密度变化量时间信号；
Figure BDA0001517884240000094
为五波长光源S发出近红外光的波长为λ2时，检测器D获得的光密度变化量时间信号；
Figure BDA0001517884240000095
为五波长光源S发出近红外光的波长为λ3时，检测器D获得的光密度变化量时间信号；
Figure BDA0001517884240000096
为五波长光源S发出近红外光的波长为λ4时，检测器D获得的光密度变化量时间信号；
Figure BDA0001517884240000097
为五波长光源S发出近红外光的波长为λ5时，检测器D获得的光密度变化量时间信号；
步骤二：对步骤一获得的五个不同波长下的近红外光在距检测器D相同距离R下的光密度变化量时间信号中每两组采用修正郎伯比尔定律构建一个方程组，总计构建十个方程组，具体的方程组表示为：
Figure BDA0001517884240000098
此公式中
Figure BDA0001517884240000099
和
Figure BDA00015178842400000910
为未知量，需要求解。εHHb(λi)、εHHb(λj)、
Figure BDA00015178842400000911
为已知参数；DPF(λi)、DPF(λj)也为已知参数(可以使文献中直接查找到的λi和λj波长下的差分路径因子；也可以是后面粒子群算法中不断迭代产生的差分路径因子值)。
其中，i为数字标号，取值范围是i＝1,2,3,4,5；j为数字标号，取值范围是j＝1,2,3,4,5；λi或λj为五波长光源S发出的近红外光的波长λ1、λ2、λ3、λ4和λ5中的某一个，同时方程组中的数字标号i和j不能取相同值；
Figure BDA00015178842400000912
和
Figure BDA00015178842400000913
表示五波长光源S发出近红外光的波长为λi或λj时，检测器D获得的光密度变化量时间信号；εHHb(λi)或εHHb(λj)为五波长光源S发出近红外光的波长为λi或λj时的还原血红蛋白消光系数；
Figure BDA0001517884240000101
或
Figure BDA0001517884240000102
为五波长光源S发出近红外光的波长为λi或λj时的氧合血红蛋白消光系数；DPF(λi)或DPF(λj)为五波长光源S发出近红外光的波长为λi或λj时的差分路径因子；
Figure BDA0001517884240000103
为采用近红外光波长组合(λi,λj)时，检测器D处的氧合血红蛋白浓度变化时间信号；
Figure BDA0001517884240000104
为采用近红外光波长组合(λi,λj)时，检测器D处的还原血红蛋白浓度变化时间信号。
步骤三：对步骤二中得到的十个方程组分别进行求解，得到十组采用不同近红外波长的氧合血红蛋白浓度变化时间信号
Figure BDA0001517884240000105
和还原血红蛋白浓度变化时间信号
Figure BDA0001517884240000106
分别表示如下：
Figure BDA0001517884240000107
Figure BDA0001517884240000108
步骤四：生成初始化的粒子群，其中第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个粒子，利用如下的公式来计算采用不同近红外波长组合下所得到氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值
Figure BDA0001517884240000111
Figure BDA0001517884240000112
其中p为数字标号，取值范围是p＝1,2,3,4,5；q为数字标号，取值范围是q＝1,2,3,4,5；λp或λq为五波长光源S发出的近红外光的波长λ1、λ2、λ3、λ4和λ5中的某一个，同时数字标号p和q不能取相同值；
Figure BDA0001517884240000113
表示为当采用第m个粒子产生的位置向量时，波长组合(λi,λj)计算得到氧合血红蛋白浓度变化时间信号；
Figure BDA0001517884240000114
表示为当采用第m个粒子产生的位置向量时，波长组合(λp,λq)计算得到氧合血红蛋白浓度变化时间信号；
Figure BDA0001517884240000115
表示为当采用第m个粒子产生的位置向量时，波长组合(λi,λj)计算得到还原血红蛋白浓度变化时间信号；
Figure BDA0001517884240000116
表示为当采用第m个粒子产生的位置向量时，波长组合(λp,λq)计算得到还原血红蛋白浓度变化时间信号；
Figure BDA0001517884240000117
表示为当采用第m个粒子产生的位置向量时，采用波长组合(λi,λj)和(λp,λq)所得到的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值；
然后，利用如下的适应度值计算公式来计算第m个粒子的适应度值：
Figure BDA0001517884240000118
步骤六：计算每个粒子的个体极值点和粒子群的全局极值点；过程为：
对于生成的初始化粒子群，每个粒子的个体极值点为当前的位置向量，对于第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次迭代更新后的各粒子的位置向量按照步骤五中的公式分别计算各粒子的适应度值；
从步骤四开始根据位置向量得到k次迭代的DPF(λi)，在根据DPF(λi)得到步骤三中的
Figure BDA0001517884240000131
和
Figure BDA0001517884240000132
然后根据步骤五中第一个公式得到
Figure BDA0001517884240000133
最终得到适应度值；
对第m个粒子，如果第k次迭代得到的位置向量xm(k)计算得到的适应度值小于前一次迭代得到的个体极值点计算得到的适应度值，则将位置向量xm(k)设置为第k次迭代得到个体极值点pbestm(k)；否则维持个体极值点不变，即pbestm(k)＝pbestm(k-1)；
如果所有粒子在第k次迭代得到个体极值点计算得到适应度值中的最小值小于前一次迭代得到全局极值点计算得到适应度值，则将该粒子(第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
其中，
Figure BDA0001517884240000141
B为增广矩阵，z为待求解矩阵；
步骤十二：对步骤十一中的增广矩阵B进行奇异值分解，即将增广矩阵B分解，表示为：
B＝UΣVH
其中，H为共轭转置；矩阵
Figure BDA0001517884240000142
为增广矩阵B的左奇异向量矩阵，数据uij为矩阵U第i行第j列对应元素，i,j＝1,2,3,4,5；矩阵
Figure BDA0001517884240000143
为增广矩阵B的右奇异向量矩阵，数据vkl为矩阵V第k行第l列对应元素，k,l＝1,2,3；矩阵
Figure BDA0001517884240000144
为对角矩阵，对角元素σ1、σ2和σ3为增广矩阵B的奇异值；
步骤十三：利用步骤十二中求取的矩阵V中第三列的对应元素，得到检测器D处的氧合血红蛋白浓度变化时间信号和检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解。
具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤十中利用步骤九中得到的各个波长下的差分路径因子最优估计值构建如下的方程组来求取氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号；
具体方程组可表示如下：
Figure BDA0001517884240000151
其中，εHHb(λ1)为光源S发出近红外光的波长为λ1时的还原血红蛋白消光系数；εHHb(λ2)为光源S发出近红外光的波长为λ2时的还原血红蛋白消光系数；εHHb(λ3)为光源S发出近红外光的波长为λ3时的还原血红蛋白消光系数；εHHb(λ4)为光源S发出近红外光的波长为λ4时的还原血红蛋白消光系数；εHHb(λ5)为光源S发出近红外光的波长为λ5时的还原血红蛋白消光系数；
Figure BDA0001517884240000152
为光源S发出近红外光的波长为λ1时的氧合血红蛋白消光系数；
Figure BDA0001517884240000153
为光源S发出近红外光的波长为λ2时的氧合血红蛋白消光系数；
Figure BDA0001517884240000154
为光源S发出近红外光的波长为λ3时的氧合血红蛋白消光系数；
Figure BDA0001517884240000155
为光源S发出近红外光的波长为λ4时的氧合血红蛋白消光系数；
Figure BDA0001517884240000156
为光源S发出近红外光的波长为λ5时的氧合血红蛋白消光系数；Δ[HbO2](t)为检测器D处的氧合血红蛋白浓度变化时间信号；Δ[HHb](t)为检测器D处的还原血红蛋白浓度变化时间信号。
其它步骤及参数与具体实施方式一相同。
具体实施方式三：本实施方式与具体实施方式一至二之一不同的是：所述步骤十三中利用步骤十二中求取的矩阵V中第三列的对应元素，得到检测器D处的氧合血红蛋白浓度变化时间信号和检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解；分别表示为：
Figure BDA0001517884240000157
Figure BDA0001517884240000158
其中，Δ[HbO2]TLS(t)为检测器D处的氧合血红蛋白浓度变化时间信号的总体最小二乘解；Δ[HHb]TLS(t)为检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解。
具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤一中5mm≤R≤40mm；
其中，R表示五波长光源S与近端检测器D之间的直线距离。
其它步骤及参数与具体实施方式一相同。
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤一中五波长光源S所发出的五种光，其各波长分别是λ1为670nm、λ2为770nm、λ3为810nm、λ4为850nm、λ5为950nm。
其它步骤及参数与具体实施方式一至四之一相同。
具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：预先设置的迭代总次数K取值为2000。
其它步骤及参数与具体实施方式一至五之一相同。
具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：所述加速因子c1设置为0.5≤c1≤5，加速因子c2设置为0.5≤c2≤5。
其它步骤及参数与具体实施方式一至六之一相同。
具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：所述加速因子典型取值可设置为c1＝c2＝2。
其它步骤及参数与具体实施方式一至七之一相同。
具体实施方式九：本实施方式与具体实施方式一至八之一不同的是：所述惯性因子w取值范围为0.4≤w≤0.9。
其它步骤及参数与具体实施方式一至八之一相同。

Claims (8)

1.基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：所述方法具体过程为：
步骤一：在待测脑组织头皮表面放置一个由五波长光源S与检测器D所构成的近红外探头，五波长光源S与检测器D之间的直线距离为R，五波长光源S发出的近红外光的波长分别为λ1、λ2、λ3、λ4和λ5，检测器D用于获取大脑安静状态下的漫反射光强和大脑诱发激励状态下的漫反射光强，从而获得五个不同波长的近红外光在距检测器D相同距离R下的光密度变化量的时间信号：
Figure FDA0002555461580000011
和
Figure FDA0002555461580000012
其中，t为采样时刻，t＝1,2,…,N，N为正整数；
Figure FDA0002555461580000013
为五波长光源S发出近红外光的波长为λ1时，检测器D获得的光密度变化量时间信号；
Figure FDA0002555461580000014
为五波长光源S发出近红外光的波长为λ2时，检测器D获得的光密度变化量时间信号；
Figure FDA0002555461580000015
为五波长光源S发出近红外光的波长为λ3时，检测器D获得的光密度变化量时间信号；
Figure FDA0002555461580000016
为五波长光源S发出近红外光的波长为λ4时，检测器D获得的光密度变化量时间信号；
Figure FDA0002555461580000017
为五波长光源S发出近红外光的波长为λ5时，检测器D获得的光密度变化量时间信号；
步骤二：对步骤一获得的五个不同波长下的近红外光在距检测器D相同距离R下的光密度变化量时间信号中每两组采用修正郎伯比尔定律构建一个方程组，总计构建十个方程组，具体的方程组表示为：
Figure FDA0002555461580000018
其中，i为数字标号，取值范围是i＝1,2,3,4,5；j为数字标号，取值范围是j＝1,2,3,4,5；λi或λj为五波长光源S发出的近红外光的波长λ1、λ2、λ3、λ4和λ5中的某一个，同时方程组中的数字标号i和j不能取相同值；
Figure FDA0002555461580000019
和
Figure FDA00025554615800000110
表示五波长光源S发出近红外光的波长为λi或λj时，检测器D获得的光密度变化量时间信号；εHHb(λi)或εHHb(λj)为五波长光源S发出近红外光的波长为λi或λj时的还原血红蛋白消光系数；
Figure FDA0002555461580000021
或
Figure FDA0002555461580000022
为五波长光源S发出近红外光的波长为λi或λj时的氧合血红蛋白消光系数；DPF(λi)或DPF(λj)为五波长光源S发出近红外光的波长为λi或λj时的差分路径因子；
Figure FDA0002555461580000023
为采用近红外光波长组合(λi,λj)时，检测器D处的氧合血红蛋白浓度变化时间信号；
Figure FDA0002555461580000024
为采用近红外光波长组合(λi,λj)时，检测器D处的还原血红蛋白浓度变化时间信号；
步骤三：对步骤二中得到的十个方程组分别进行求解，得到十组采用不同近红外波长的氧合血红蛋白浓度变化时间信号
Figure FDA0002555461580000025
和还原血红蛋白浓度变化时间信号
Figure FDA0002555461580000026
分别表示如下：
Figure FDA0002555461580000027
Figure FDA0002555461580000028
步骤四：生成初始化的粒子群，其中第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个粒子，利用如下的公式来计算采用不同近红外波长组合下所得到氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值
Figure FDA0002555461580000031
Figure FDA0002555461580000032
其中p为数字标号，取值范围是p＝1,2,3,4,5；q为数字标号，取值范围是q＝1,2,3,4,5；λp或λq为五波长光源S发出的近红外光的波长λ1、λ2、λ3、λ4和λ5中的某一个，同时数字标号p和q不能取相同值；
Figure FDA0002555461580000033
表示为当采用第m个粒子产生的位置向量时，波长组合(λi,λj)计算得到氧合血红蛋白浓度变化时间信号；
Figure FDA0002555461580000034
表示为当采用第m个粒子产生的位置向量时，波长组合(λp,λq)计算得到氧合血红蛋白浓度变化时间信号；
Figure FDA0002555461580000035
表示为当采用第m个粒子产生的位置向量时，波长组合(λi,λj)计算得到还原血红蛋白浓度变化时间信号；
Figure FDA0002555461580000036
表示为当采用第m个粒子产生的位置向量时，波长组合(λp,λq)计算得到还原血红蛋白浓度变化时间信号；
Figure FDA0002555461580000037
表示为当采用第m个粒子产生的位置向量时，采用波长组合(λi,λj)和(λp,λq)所得到的氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号之间的差异值；
然后，利用如下的适应度值计算公式来计算第m个粒子的适应度值：
Figure FDA0002555461580000041
步骤六：计算每个粒子的个体极值点和粒子群的全局极值点；过程为：
对于生成的初始化粒子群，每个粒子的个体极值点为当前的位置向量，对于第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)的最优估计值；
步骤十，利用步骤九中得到的各个波长下的差分路径因子最优估计值构建如下的方程组来求取氧合血红蛋白浓度变化时间信号和还原血红蛋白浓度变化时间信号；
具体方程组可表示如下：
Figure FDA0002555461580000061
其中，εHHb(λ1)为光源S发出近红外光的波长为λ1时的还原血红蛋白消光系数；εHHb(λ2)为光源S发出近红外光的波长为λ2时的还原血红蛋白消光系数；εHHb(λ3)为光源S发出近红外光的波长为λ3时的还原血红蛋白消光系数；εHHb(λ4)为光源S发出近红外光的波长为λ4时的还原血红蛋白消光系数；εHHb(λ5)为光源S发出近红外光的波长为λ5时的还原血红蛋白消光系数；
Figure FDA0002555461580000062
为光源S发出近红外光的波长为λ1时的氧合血红蛋白消光系数；
Figure FDA0002555461580000063
为光源S发出近红外光的波长为λ2时的氧合血红蛋白消光系数；
Figure FDA0002555461580000064
为光源S发出近红外光的波长为λ3时的氧合血红蛋白消光系数；
Figure FDA0002555461580000065
为光源S发出近红外光的波长为λ4时的氧合血红蛋白消光系数；
Figure FDA0002555461580000066
为光源S发出近红外光的波长为λ5时的氧合血红蛋白消光系数；Δ[HbO2](t)为检测器D处的氧合血红蛋白浓度变化时间信号；Δ[HHb](t)为检测器D处的还原血红蛋白浓度变化时间信号；
步骤十一：将步骤十中的方程组改写为如下的矩阵形式：
Bz＝0
其中，
Figure FDA0002555461580000067
B为增广矩阵，z为待求解矩阵；
步骤十二：对步骤十一中的增广矩阵B进行奇异值分解，即将增广矩阵B分解，表示为：
B＝UΣVH
其中，H为共轭转置；矩阵
Figure FDA0002555461580000071
为增广矩阵B的左奇异向量矩阵，数据uij为矩阵U第i行第j列对应元素，i,j＝1,2,3,4,5；矩阵
Figure FDA0002555461580000072
为增广矩阵B的右奇异向量矩阵，数据vkl为矩阵V第k行第l列对应元素，k,l＝1,2,3；矩阵
Figure FDA0002555461580000073
为对角矩阵，对角元素σ1、σ2和σ3为增广矩阵B的奇异值；
步骤十三：利用步骤十二中求取的矩阵V中第三列的对应元素，得到检测器D处的氧合血红蛋白浓度变化时间信号和检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解。
2.根据权利要求1所述基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：所述步骤十三中利用步骤十二中求取的矩阵V中第三列的对应元素，得到检测器D处的氧合血红蛋白浓度变化时间信号和检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解；分别表示为：
Figure FDA0002555461580000074
Figure FDA0002555461580000075
其中，Δ[HbO2]TLS(t)为检测器D处的氧合血红蛋白浓度变化时间信号的总体最小二乘解；Δ[HHb]TLS(t)为检测器D处的还原血红蛋白浓度变化时间信号的总体最小二乘解。
3.根据权利要求2所述基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：所述步骤一中5mm≤R≤40mm；
其中，R表示五波长光源S与近端检测器D之间的直线距离。
4.根据权利要求3所述基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：所述步骤一中五波长光源S所发出的五种光，其各波长分别是λ1为670nm、λ2为770nm、λ3为810nm、λ4为850nm、λ5为950nm。
5.根据权利要求4所述基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：预先设置的迭代总次数K取值为2000。
6.根据权利要求5所述基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：所述加速因子c1设置为0.5≤c1≤5，加速因子c2设置为0.5≤c2≤5。
7.根据权利要求6所述基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：所述加速因子c1＝c2＝2。
8.根据权利要求7所述基于差分路径因子估计的近红外脑功能信号处理方法，其特征在于：所述惯性因子w取值范围为0.4≤w≤0.9。
CN201711392782.8A 2017-12-21 2017-12-21 基于差分路径因子估计的近红外脑功能信号处理方法 Active CN107928631B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711392782.8A CN107928631B (zh) 2017-12-21 2017-12-21 基于差分路径因子估计的近红外脑功能信号处理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711392782.8A CN107928631B (zh) 2017-12-21 2017-12-21 基于差分路径因子估计的近红外脑功能信号处理方法

Publications (2)

Publication Number Publication Date
CN107928631A CN107928631A (zh) 2018-04-20
CN107928631B true CN107928631B (zh) 2021-05-07

Family

ID=61941617

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711392782.8A Active CN107928631B (zh) 2017-12-21 2017-12-21 基于差分路径因子估计的近红外脑功能信号处理方法

Country Status (1)

Country Link
CN (1) CN107928631B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109106376B (zh) * 2018-09-20 2022-04-19 京东方科技集团股份有限公司 一种血液中总血红蛋白浓度的检测方法及装置
CN109753797B (zh) * 2018-12-10 2020-11-03 中国科学院计算技术研究所 针对流式图的密集子图检测方法及系统
CN110192866A (zh) * 2019-04-28 2019-09-03 上海爱德赞医疗科技有限公司 无创毛细动脉血液组分浓度的监测方法及设备
CN112990224B (zh) * 2021-02-04 2023-07-11 郑州航空工业管理学院 一种用于从fMRI数据中识别脑功能划分的粒子群方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102499651A (zh) * 2011-10-24 2012-06-20 华南理工大学 用于监护系统的报警方法
CN102973279B (zh) * 2012-12-18 2014-09-17 哈尔滨工业大学 独立成分分析联合最小二乘法的近红外脑机接口的信号检测方法
CN106919948B (zh) * 2015-12-28 2021-04-06 西南交通大学 一种驾驶持续性注意水平的识别方法
CN106371610B (zh) * 2016-09-23 2020-06-09 重庆金瓯科技发展有限责任公司 一种基于脑电信号的驾驶疲劳的检测方法
CN106821355A (zh) * 2017-04-01 2017-06-13 泰康保险集团股份有限公司 血压预测的方法及装置
CN107174204B (zh) * 2017-05-12 2020-07-24 哈尔滨工业大学 基于总体最小二乘法的近红外脑功能信号处理方法

Also Published As

Publication number Publication date
CN107928631A (zh) 2018-04-20

Similar Documents

Publication Publication Date Title
CN107928631B (zh) 基于差分路径因子估计的近红外脑功能信号处理方法
KR20170112404A (ko) 원격탐사 초분광 영상이미지를 이용한 병원성 미생물의 농도 예측방법
ES2960485T3 (es) Procedimiento de medición de la concentración de un objeto de medición biométrica mediante el uso de un aprendizaje profundo de inteligencia artificial
CA2481650A1 (en) Modification of the normalized difference method for real-time optical tomography
Habermehl et al. Optimizing the regularization for image reconstruction of cerebral diffuse optical tomography
CN109993808B (zh) 一种基于dsn的动态双示踪pet重建方法
Garrido et al. Introduction of deep learning in thermographic monitoring of cultural heritage and improvement by automatic thermogram pre-processing algorithms
CN111728590A (zh) 基于动态功能连接的个体认知能力预测方法和系统
CN112798654B (zh) 用于电阻抗层析成像的快速梯度法和自适应雅可比矩阵重构方法
CN110068543A (zh) 一种基于迁移学习的太赫兹光谱识别方法
CN110969614A (zh) 基于三维卷积神经网络的脑龄预测方法及系统
Wang et al. Fluorescence molecular tomography reconstruction of small targets using stacked auto-encoder neural networks
CN108846200B (zh) 一种基于迭代法的准静态桥梁影响线识别方法
CN109799472B (zh) 一种基于深度学习的磁共振涡流补偿方法
CN115015126B (zh) 一种粉末状生物粒子材料活性判定方法和系统
Zhang et al. Unrolled convolutional neural network for full-wave inverse scattering
Oldan et al. Channelized Hotelling and human observer study of optimal smoothing in SPECT MAP reconstruction
CN117116434B (zh) 人脑白质结构连接组的个体差异评估方法、应用及装置
CN112581385A (zh) 基于多先验约束扩散峰度成像张量估计方法、介质和设备
CN117041972A (zh) 基于信道-时空注意自编码的车联网传感器异常检测方法
CN107174204B (zh) 基于总体最小二乘法的近红外脑功能信号处理方法
CN111854960B (zh) 一种基于热成像仪测量人体脸部温度校准方法
CN114255183A (zh) 数据与知识联合驱动的智能计算光学成像方法
CN113870377A (zh) 基于V-ResNet肺部成像方法
CN117454200A (zh) 基于云计算的新生儿黄疸筛查方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant