一种基于分层重建和多波长图像融合的脑功能成像方法

本发明涉及近红外脑功能成像，特别是涉及一种基于分层重建和多波长图像融合的脑功能成像方法。

背景技术：

1、近红外脑功能成像利用近红外光入射大脑组织，探测分析经由大脑组织吸收散射后再次出射的扩散光，计算大脑内部光学参数变化信息，并据此获得大脑兴奋分布图像。作为一种主动测量技术，近红外脑功能成像具有无创、非电离、合理时间/空间分辨率等优点，已逐渐成为传统脑功能成像技术(如核磁共振成像、脑电图等)的重要补充，越来越多的应用于脑功能临床研究[1]。

2、扩散光学层析成像(diffuse optical tomography,dot)作为近红外脑功能成像技术的重要发展，具备三维(three dimension,3d)成像能力，可以为脑功能成像提供更加丰富的信息。由于大脑组织的高散射特性和测量数据量的局限性，dot计算过程具有严重的不适定性，很难获得准确解。为了缓解dot技术的不适定性，可以假设在特定层的深度方向吸收系数保持不变，将3d问题变为二维(two dimension,2d)问题，大幅度减少未知量的个数，该技术称为半三维(semi three dimension,s3d)dot重建。此外，由于采用反射测量模式，dot技术对表层生理干扰更加敏感，这些生理干扰和测量过程带来的随机噪声在某些情况下可能“湮没”大脑兴奋信号。

3、为了提高近红外脑功能成像灵敏度，实际中常将带通滤波和多周期平均相结合来有效抑制生理干扰和随机噪声。然而，带通滤波频带选择往往具有主观性，当生理干扰和大脑兴奋频带重叠时将带来严重误差。虽然多周期平均是应对随机噪声的有效手段，但却带来冗长的测量时间。一些研究者提出了利用小光源-探测器距离(source-detector-separation,sds)通道测量信号或者辅助测量信号作为参考信号，完成自适应滤波过程实现干扰和噪声的有效抑制[2-4]。但是，小sds通道测量信号中可能包含深层兴奋信息，以此作为干扰和噪声的估计信号进行滤波将引入误差，且如何为各个测量通道提供合理的参考信号也是实践中的一个难点。此外，辅助测量会增加系统的复杂性，其调试配准过程将产生额外的时间成本。

4、当前，临床上对高速高灵敏脑功能成像的需求越来越大，这对近红外脑功能成像中的干扰和噪声滤除算法提出了更高的要求：无需多周期重复和辅助测量即可有效抑制生理干扰和随机噪声。

技术实现思路

1、本发明的目的是针对现有技术中存在的技术缺陷，而提供一种基于分层重建和多波长图像融合的脑功能成像方法。该方法可以同时抑制生理干扰和随机噪声，提高重建图像质量。在较少甚至单一周期下获得与传统多周期平均方案相当的重建质量，有效减少测量时间。

2、为实现本发明的目的所采用的技术方案是：

3、一种基于分层重建和多波长图像融合的脑功能成像方法，包括以下步骤：

4、步骤1：对成像区域进行离散化，并计算其3d敏感度矩阵；建立各通道光强变化与各离散节点吸收系数变化的线性关系；

5、步骤2：将成像区域沿着深度方向划分，并对各层进行不同空间密度的剖分，得到3d立方体网格；

6、步骤3：将步骤1获得的3d敏感度矩阵插值到步骤2划分的3d立方体网格中，并沿着深度方向在各层内进行积分，得到各层的s3d敏感度矩阵；建立各通道光强变化与各层吸收系数变化的线性关系；

7、步骤4：特定波长及特定时间点下，将表层的吸收系数变化插值到与深层相同的空间密度，利用插值后的表层吸收系数变化对深层吸收系数变化进行空间对应的自适应滤波，获得滤波后的深层吸收系数变化，即特定波长下深层吸收图像；

8、步骤5：对各个波长下深层吸收图像进行小波分解，获得各个波长下深层吸收图像的特征分量，然后对所得特征分量进行融合，最后利用小波反变换获得修正模板；

9、步骤6：利用所得修正模板对各个波长下深层吸收图像进行修正并融合，最终获得融合的大脑兴奋图像。

10、在上述技术方案中，步骤1中，对成像区域按照交叠阵列进行离散化；

11、所述交叠阵列包括20个光源点和12个探测点。

12、在上述技术方案中，步骤1中，根据born近似或rytov近似建立各通道光强变化与各个离散节点吸收系数变化的线性方程：

13、m＝j·δμa                                      (1)；

14、其中：

15、δμa＝[δμ1 δμ2 … δμn]t，表示各个离散节点吸收系数变化，n表示离散节点总数；

16、m＝[m(ξ1,ζ1) m(ξ2,ζ1) … m(ξd,ζs)]t，表示各通道光强变化，ζs(s＝1,2，…，s)表示总数s个光源位置，ξd(d＝1,2，…，d)表示总数d个探测器位置；

17、表示3d敏感度矩阵。

18、在上述技术方案中，各通道光强变化m中各元素计算式为：

19、

20、式中，γ表示测量到的光强值；γ0表示参考光强值，通常取基线态测量值的平均。

21、在上述技术方案中，3d敏感度矩阵j各元素计算式为：

22、

23、式中，φ0表示不同位置的光子密度，g0为格林函数，为第n个节点的3d形函数，ν表示光在组织中传输的速率，ω表示重建区域，r′表示重建区域内空间位置信息。

24、在上述技术方案中，步骤2中，将成像区域沿着深度方向划分为sc层、sk层和cc层；所述sc层即表层；所述cc层即深层。

25、在上述技术方案中，所述sc层采用稀疏剖分网格；所述cc层采用密集剖分网格。

26、在上述技术方案中，步骤3中，各通道光强变化与各层吸收系数变化的线性方程：

27、

28、其中，λw代表特定波长；ti代表特定时间点；

29、m1nn(λw,ti)和m2nn(λw,ti)分别表示特定波长及特定时间点下1nn测量通道和2nn测量通道的光强变化向量；

30、和分别代表1nn测量通道/2nn测量通道中sc层/cc层在特定波长下的s3d敏感度矩阵；

31、和分别代表特定波长及特定时间点下，sc层和cc层的吸收系数变化。

32、在上述技术方案中，步骤4中，第p个节点处，滤波后的cc层吸收图像为：

33、

34、其中表示滤波后结果，表示k阶的参考信号，θ(p,ti-1)＝[θ1(p,ti-1),…,θk(pw,ti-1)]t表示权重因子，其初始值为0。

35、在上述技术方案中，利用所得修正模板i(ti)对各个波长下cc层吸收图像进行修正并融合，最终获得融合图像

36、

37、其中⊙代表对应元素相乘。

38、与现有技术相比，本发明的有益效果是：

39、1.本发明提供的脑功能成像方法，利用交叠阵列获得不同sds的测量数据，并基于s3d-dot技术实现分层重建：表层采用稀疏网格以求在一定程度上通过区域平均的方法抑制来自深层兴奋信号的影响，从而准确估计表层干扰信号；底层采用密集网格以求获得大脑兴奋分布的高空间分辨能力；接着，以插值后表层重建结果为参考对深层重建结果进行自适应滤波，有效抑制生理干扰对大脑兴奋分布的影响，明显提高重建结果的空间分辨能力。

40、2.本发明提供的脑功能成像方法，对完成自适应滤波后各个波长的重建图像进行小波分解，然后对各个波长的特定分量进行融合，接着利用小波反变换获得修正模板。最后，利用归一化修正模板对各个波长的滤波图像进行修正，基于图像融合算法实现随机噪声的有效抑制，获得高空间分辨率的大脑兴奋重建图像。

41、3.本方法具有严密的理论模型，可以同时抑制生理干扰和随机噪声对大脑兴奋分布重建结果的影响。过程完全由数据驱动，无需额外辅助测量手段。无需重复测量即可获得良好的图像质量，极大提高脑功能成像效率，在基于近红外脑功能成像的脑机接口领域具有良好的应用前景。