[发明专利]基于五层脑结构模型的OT导引脑功能半三维DOT方法

摘要

本发明公开了一种基于五层脑结构模型的OT导引脑功能半三维DOT方法，主要包括设计一种五层脑结构模型，用辐射传输方程和扩散方程分别对光子在非散射区域和散射区域中的传输进行建模且将他们相互耦合以准确描述光子在整个五层脑结构模型中的传输。利用同一组测量数据，先用OT粗略重建提供先验位置信息，再对DOT逆问题进行正则化以改善其病态性，提高重建图像的质量。假设吸收系数的变化量在灰质层内不随深度变化，将三维逆问题简化成一个二维图像重建问题，即半三维图像重建。本发明充分结合了OT鲁棒性强和DOT空间分辨率高、定量性强的优势，简化了测量过程和图像重建过程，改善了成像质量。

主权项

一种基于五层脑结构模型的OT导引脑功能半三维DOT方法，其包括以下步骤：步骤1、构建一个采用交叠采样方式的测量系统，包括：设计一个五层脑结构模型，即：所述五层脑结构模型所占的空间是在三维坐标系中，X、Y和Z方向分别为0‑100mm，0‑120mm和0‑32mm；所述五层脑结构模型由自外而内依次布置的头皮层、头骨层、脑脊液层、灰质层和白质层构成；在所述五层脑结构模型的头皮层上布配20个光源点，20个光源点按照4行5列以20mm距离等间距呈网格状布置，并逐行自左向右对光源点进行编号，光源点编号i＝1，2，3，…，20；在每4个相邻光源点的几何中心布置一个探测器，从而在头皮层上形成了12个探测点；每个光源点上设有1至3个选自600～1000nm范围内的不同波长的光源，波长分别表示为λj，j＝1，2，3；针对每个光源点，将与该光源点距离分别为L1＝14.14mm、L2＝31.62mm、L3＝42.40mm和L4＝50.99mm的所有探测器与该光源点之间形成的通道依次定义为第一近邻通道、第二近邻通道、第三近邻通道和第四近邻通道；其中包括有44个第二近邻通道，在所有第二近邻通道的中点处均设有一采样点；步骤2、根据射线‑扩散理论建立针对所述五层脑结构模型的光传输模型：利用如下偏微分方程和边界条件对光子在五层脑结构模型中的传输进行描述，-D(r)▿2Φ(r)+μa(r)Φ(r)=q0(r),r∈Ω,(a)Φ(r)+2D(r)∂Φ(r)∂n^=12π∫∂Ω+Φ(r′)G(r-r′)dr′,r′∈∂Ω+,r∈∂Ω-,(b)Φ(r)+2D(r)∂Φ(r)∂n^=12π∫∂Ω-Φ(r′)G(r-r′)dr′,r′∈∂Ω-,r∈∂Ω+,(c)Φ(r)+2D(r)∂Φ(r)∂n^=0,r∈∂Ω/(∂Ω-∪∂Ω+(d)---(1)]]>式(1)中，(a)表示光子在头皮层、头骨层、灰质层和白质层中传输时需要满足的扩散方程；(b)、(c)和(d)分别表示边界条件；r表示空间某位置；μa(r)和μ′s(r)分别表示r位置处的吸收系数和约化散射系数；和Φ(r)分别表示r位置处的扩散系数和光子密度；q0(r)表示光源；Ω表示整个待考察的散射区域，包括所述五层脑结构模型中的头皮层、头骨层、灰质层和白质层；表示待考察散射区域的所有边界，和分别表示所述五层脑结构模型中脑脊液层的上下界面；表示界面和处指向非散射区域的法向量；步骤3、由光学拓扑成像导引脑功能半三维扩散光学层析成像，包括：3‑1)利用光学拓扑成像获得一幅反映吸收系数变化量的二维粗略图像：用光源照射各个光源点，先后在静息状态和任务状态下测量所有第二近邻通道中的探测器位置处出射光强度，并分别记作和其中，λj表示编号为j的波长；i表示第i个第二近邻通道的编号；R和T分别表示静息状态和任务状态；用式(2)计算在第i个第二近邻通道采样点的吸收系数变化量△μa；Δμa≈-1Blln(IT,iλjIR,iλj)---(2)]]>式(2)中，B表示差分路径因子，l表示第i个第二近邻通道的长度；然后将它们插值在五层脑结构模型水平截面之内有限元剖分栅格上，从而获得一幅反映吸收系数变化量的二维粗略图像N表示五层脑结构模型的水平截面之内有限元剖分栅格节点的个数；3‑2)对光学拓扑成像结果进行图像分割获取吸收系数变化量的先验位置信息：以吸收系数变化量最大值的30％作为阈值，将上述二维粗略图像中吸收系数变化量大于该阈值的有限元剖分栅格节点赋值为1，其余节点赋值为0，获取一幅二值图像；赋值为1的区域表示感兴趣区域，记作ROI，赋值为0的区域表示非感兴趣区域，记作non‑ROI；3‑3)将三维扩散光学层析成像简化成半三维扩散光学层析成像：通过式(3)将三维扩散光学层析成像逆问题简化成二维图像重建问题ln[IT(ξd,ζs)IR(ξd,ζs)]=-1I(ξd,ζs)∫Rxy{[∫ZgG(ξd,r)Φ(r,ζs)dz]Δμag(x,y)}dxdy---(3)]]>式(3)中，ξd(d＝1…D,D＝12)和ξs(s＝1…S,S＝20)分别代表头皮层上12个探测点位置和20个光源点位置；IT(ξd,ζs)和IR(ξd,ζs)分别表示在任务状态和静息状态下在ξd处探测所得的由ξs处光源激励所引起的出射光强度；I(ξd,ζs)表示静息状态下根据式(1)正向模拟所得的当光源处于ξs，在ξd处测量得到的出射光强度；G(ξd,r)和Φ(r,ζs)分别表示当光源放在ξd和ζs处时，所对应的光子密度的空间分布；△μag(x,y)表示灰质层中(x，y)处吸收系数的变化量；Zg表示灰质层在Z‑轴方向的范围；Rxy表示灰质区在X‑Y平面的投影；扩散光学层析成像图像重建问题最终表示为如下矩阵等式M＝W△μa                           (4)式(4)中，△μa表示在灰质层内吸收系数变化量的二维分布，W是边界光流量变化量相对于各有限元剖分栅格节点吸收系数变化量的权重矩阵，其元素如下：W[(s-1)D+d,n]=∫Rxy{[∫ZgG(ξd,r)Φ(r,ζs)dz]Δμag(x,y)}dxdy---(5)]]>式(5)中，s,d和n分别表示光源、探测器和有限元剖分网格节点的序号；3‑4)在半三维扩散光学层析成像中引入先验位置信息，最终获得反映灰质层内各位置处吸收系数变化量的二维图像：采用软先验正则化策略通过式(6)将步骤3‑2)获得的先验位置信息引入到扩散光学层析成像逆问题的求解中Δμ^a=argmin{||M-WΔμa||22+α||LΔμa||22}---(6)]]>式(6)中，L是依赖于ROI的拉普拉斯算子型正则化矩阵Lij=0ifΔμa(b)(xi,yi)≠Δμa(b)(xj,yj)-1/Nbifi≠jandΔμa(b)(xi,yi)=Δμa(b)(xj,yj)=b1ifi=j---(7)]]>式(7)中，Nb是ROI(b＝1)和non‑ROI(b＝0)内像素的个数；α是根据选择所得的吉洪诺夫正则化参数，其中，{γi|i＝1…N}，{μi|i＝1…N}分别是矩阵W和L的奇异值；式(6)所示的最优化问题简化为如下的矩阵等式WαLΔμa=M0---(8)]]>使用代数重建技术求解式(8)所示的矩阵方程即可得到高空间分辨、高定量性的反映灰质层内各位置处吸收系数变化量的二维图像。