小动物活体多模态分子成像系统和成像方法

摘要

本发明公开了一种小动物活体多模态分子成像系统和成像方法，该成像系统包括切伦科夫发光成像装置、核素成像装置和MicroCT装置；这三个装置分别获取小动物体表的切伦科夫光学信号、核素信号和X光信号。采用锥束反投影算法对X光信号进行重建，获取小动物的解剖结构信息；基于混合光传输理论，对切伦科夫光学信号进行重建，获取切伦科夫靶向目标的三维分布信息；将获取的切伦科夫靶向目标信息与获取的解剖结构信息进行融合并显示，完成多模态成像。本发明实现了一种集核素成像的高特异性、光学成像的高灵敏性和X光成像的高分辨率为一体的单标记探针/多模态成像方式，使获得的小动物在体多模影像信息更加完整、全面和精确。

说明书

技术领域

本发明属于图像处理技术领域，涉及基于放射线核素的核素成像和切伦科夫发光 成像，MicroCT成像，尤其是多模态分子成像装置和光学成像三维成像方法，可用于 小动物在体分子成像。

背景技术

分子影像技术在过去十年中取得了飞速的发展，已经广泛应用于肿瘤检测和药物 研发等小动物实验和预临床实验。目前其常用成像技术主要包括核素成像、MRI和光 学成像等，这些模态各自不同的特点，在系统特异性、灵敏度、设备成本等方法有着 各自的优势和不足。核素成像，如PET、SPECT成像灵敏度高、不受探测深度的限制， 但其设备昂贵，空间分辨率低；MRI空间分辨率和组织分辨力很高，不受组织深度的 影响，但其敏感性很低；光学成像，如生物自发光成像、激发荧光成像具有灵敏度高、 成本低、成像速度快的特点，其缺点是穿透深度低，空间分辨率低，单一成像模态无 法为疾病诊断提供足够全面的信息。为了实现对重大疾病的早期诊断与精确诊断，结 合不同成像模式的优势获得更灵敏、更准确、更全面的生理病理信息的多模态分子成 像成为分子影像发展的热点和发展趋势。

多模态靶向分子探针是多模态分子成像的重要基础，靶向分子与靶标结合的特异 性和亲和力是决定分子成像探针成败的首要因素。传统多模态分子探针制备工艺复 杂，多次进行化学合成或标记，影响分子探针与靶分子的亲和力、特异性及稳定性等 生物学行为等，制约了多模态分子影像技术的快速发展。目前一种新的成像方式—— 切伦科夫发光成像CLI，为多模态分子影像技术的发展提供了新的研究思路。切伦科 夫效应是指带电粒子以超光速在介质运行时，带电粒子可将其能量的一部分转化为可 见光。多种核素的高能射线均可通过切伦科夫效应产生光学信号继而被光学成像系统 探测，因而切伦科夫发光成像在2009年首先报道后便受到极大关注，参见：R.Robertson， M.S.Germanos，C.Li，G.S.Mitchell，S.R.Cherry，M.D.Silva，“Optical imaging of  Cerenkov light generation from positron-emitting radiotracers”，Physics in Medicine and  Biology，54：N355-65(2009).切伦科夫发光成像实现了放射性核素的光学成像系统检 测，使得光学成像的研究对象从传统光学探针拓展到放射线核素探针，拓宽了光学分 子成像的研究方向，为其临床转化提供了可能。同时，核素-切伦科夫发光多模态分 子成像，实现了集核素显像的高特异性、光学成像的高灵敏性为一体的单标记探针/ 多模态成像，将是小动物在体分子成像的重要发展方向。

基于核素成像和光学成像的多模态成像系统，美国Carestream Health公司在其专 利申请文件“Apparatus and method for multi-modal imaging”(多模态成像的装置和方法) (美国专利号773425)中提出了一种基于光学/核素/X射线成像的多模态成像装置。 所有模态的数据采集全部由同一个CCD探测器完成，用X光磷屏和同位素磷屏完成X 光或核素信号向可见光的转变。该专利第一次实现了核素成像、光学成像和X光成像 的多模态融合，但该系统的成像目标始终固定平放在系统的样品台上，不能获得多角 度的成像信息。因此仅仅是在二维上对数据进行采集、后处理和分析，不能获取靶向 目标的准确定位信息，不适合分子影像领域的发展需求。清华大学在其专利申请文件 “旋转平台式小动物在体多模成像检测系统”(公开号CN1994229，公开日2007.07.11， 授权号ZL 200610144284.7，授权日2009.01.14)中提出了一种在体多模成像系统，包 括红外光、可见光和紫外光成像系统，以及核素成像装置。但该专利是针对不同模态 的分子探针进行成像，不同探针标记的位点和浓度存在很大差异，在数据后处理与分 析上会带来不可预知的误差。并且该专利中的系统皆使用普通相机获取的小动物体多 角度白光图像进行表面轮廓成像，不能获得小动物体的解剖结构信息，采集的多个角 度图像也只是进行了多模配准和叠加，很难实现对靶向目标的定位和定量成像。

核素和切伦科夫发光多模态成像方法中，核素成像的三维重建算法已经非常成 熟，切伦科夫发光成像的三维重建理论和算法属于光学断层成像范畴。切伦科夫发光 信号的光谱范围是400-900nm，能量更多地集中在高频谱段。研究人员提出了将宽谱 信号进行分谱段扫描成像的方法。可以参见：S.Antonello，K.Chaincy，R.Brad， C.Riccardo，M.Pasquina，S.Andrea，B.Federico，Multispectral Cerenkov luminescence  tomography for small animal optical imaging，Optics Express，19(13)：12605-18(2011). 该方法利用一组6块带宽为20nm的滤波片获取多谱段信息，对三维成像精度有一定的 改善。但滤波片太少不能覆盖所有的波段，提供的检测信息不够完整和准确，太多会 使得采集信号所需时间过长会导致多光谱信号的不一致。因此要考虑一种能够实时进 行多光谱探测的方法，才能保证既不丢失多光谱信号，也能保证不会因为时间过长而 导致信号不准确。

切伦科夫发光成像的三维成像方法中需要建立一个能精确描述光在任意复杂生 物体中传输的数学模型。研究人员基于辐射传输方程的低阶扩散近似模型开展了切伦 科夫发光成像的三维重建，参见：Z.Hu，J.Liang，W.Yang，J.Tian，Experimental  Cerenkov luminescence tomography of the mouse model with SPECT imaging validation. Optics Express.18(24)：24441-24450(2010).但扩散近似只能准确描述光在高散射生 物组织中的传输过程，对于低散射特性组织具有很低的求解精度。辐射传输方程的高 阶近似可以对同时存在高散射和低散射特性组织的小动物体进行准确成像，参见：J. Zhong，J.Tian，X.Yang，C.Qin，Whole-Body Cerenkov Luminescence Tomography with  the Finite Element SP3 Method，Annals of Biomedical Engineering，39：1728-1735(2011). 该方法能够保证一定的求解精度，但高阶近似方程求解难度非常大，对于具有复杂解 剖结构的非匀质光传输模型，该方法计算复杂度将远远超过实际承受能力。

切伦科夫发光信号在生物组织中的传输包括大量的散射和吸收作用，光能量会随 着传输深度的增加而快速的减小，从而体表的测量信号很弱，不能准确、灵敏地反映 深层组织中的目标信息。

发明内容

本发明的目的之一是针对上述成像系统的问题，提出一种小动物活体多模态分子 成像系统，以单一标记探针完成三维多模态成像，实现对切伦科夫发光成像、核素成 像和X光成像的信号的同时采集。

本发明的目的之二是针对上述成像方法的问题，提出一种小动物活体多模态分子 成像方法，以完成实时多光谱信号探测，并且考虑了光能量随深度信息的变化情况， 获得了更加完整、全面和精确的小动物在体多模态影像信息。

为实现上述目的，本发明提出小动物活体多模态分子成像系统，包括信号采集子 系统和计算机子系统，其特征在于：

信号采集子系统，包括切伦科夫发光成像装置、核素成像装置、MicroCT装置、 小动物支架和旋转平移台；该切伦科夫发光成像装置包括高灵敏度的制冷CCD探测 器、扫描式光谱仪装置和滤波片组，滤波片组位于高灵敏度的制冷CCD探测器的正前 方，且两者与扫描式光谱仪装置以旋转台为中心在空间呈一条直线安装在实验台上； 该核素成像装置包括一对PET探测器，它们以旋转台为中心在空间呈一条直线安装在 实验台上；该MicroCT装置包括X射线发射器和X射线平板探测器，两者以旋转台为 中心在空间呈一条直线安装在实验台上；所述的切伦科夫发光成像装置、核素成像装 置和MicroCT装置在旋转台四周对应均匀分布；

计算机子系统，包括系统控制模块、信号处理模块和数据重建模块；该系统控模 块，用于对切伦科夫发光成像装置、核素成像装置和MicroCT装置进行控制，以完成 各个装置与计算机之间的信号传输；该信号处理模块，用于对MicroCT探测到的X光 信号和PET探测到的核素信号进行增强和去噪处理；该数据重建模块，用于对MicroCT 探测到的X光信号和PET探测到的核素信号进行重建，得到小动物体的解剖结构信息 和小动物体内核素靶向目标分布信息，并且对切伦科夫发光成像装置探测到的光学信 号进行重建，得到小动物体内的切伦科夫靶向目标的分布信息。

所述的旋转平移台，包括一个旋转台、一对相互垂直的平移台和一个旋转平移控 制箱，旋转台正上方固定小动物支架。

所述的数据重建模块，包括：切伦科夫光学信号重建子模块，用于将切伦科夫发 光成像装置探测到的伦科夫光学信号，转化为小动物体内的切伦科夫靶向目标的分布 信息；核素信号重建子模块，用于将核素成像装置探测到的核素信号，转化为小动物 体内的核素靶向目标的分布信息；X光信号重建子模块，用于将MicroCT装置探测到 的X光信号，转化为小动物的解剖结构信息。

为实现上述目的，本发明提出的小动物活体多模态分子成像方法，包括如下步骤：

(1)多光谱数据的获取

在高灵敏制冷CCD探测器前放置不同波段的带通滤光片组，利用CCD探测器采集 不同光谱段的小动物体表切伦科夫发光光通量密度，同时采集带有支架标记点信息的 白光图像以及背景图像；利用光谱仪装置对体表切伦科夫光学信号进行扫描探测，得 到该波段发光信号的光谱分布；将不同光谱段的切伦科夫发光光通量密度离散成更窄 的波段上的光通量密度，最终得到小动物表面测得的切伦科夫发光光通量密度为Φ^m。

(2)数据预处理及配准

利用采集到的背景图像对切伦科夫光学信号进行统计分析和去噪预处理；利用锥 束反投影算法对MicroCT装置探测的X光信号进行重建，获取小动物的三维体数据； 利用小动物支架上的球形标记点，对切伦科夫光学信号和小动物的三维体数据进行坐 标匹配和能量映射；然后经过体数据分割，获得小动物的三维结构信息；

(3)将多光谱信息与混合多级自适应有限元相结合，重建切伦科夫靶向目标：

(3a)根据辐射传输方程RTE高低阶近似的混合光传输理论，用扩散方程DE近似 描述切伦科夫光学信号在高散射特性组织区域的传输过程，用简化球谐函数SP_N描述 光学信号在非高散射特性组织区域的传输过程；

(3b)划定先验光源可行区，利用混合多级自适应有限元网格细分策略，在第k 级离散网格上，结合扩散方程DE和简化球谐函数SP_N在不同特性组织区域的光传输模 型的近似，使用罗宾边界条件进行耦合，最终得到第k级离散网格上内部靶向目标能 量密度S_k和表面测得的切伦科夫信号光通量密度之间的线性关系方程；

$A_k{S}_k={\Phi }_k^m;$

其中A_k为第k级离散网格上的系统矩阵；

(3c)根据深度的变化将系统矩阵A_k的子矩阵顺序排列，每个子矩阵与某一深 度值对应，求得各子矩阵的最大特征值；通过逆向排列各子矩阵的最大特征值，建立 其深度补偿矩阵C；

(3d)校正深度变化对第k级离散网格上的系统矩阵A_k的影响，得到深度补偿后 的第k级网格上的系统矩阵

$A_k^{*}=A_{k}C;$

(3e)利用深度补偿后系统矩阵建立内部靶向目标能量密度S_k和表面实际 测得的光通量密度之间的新的线性关系方程；

$A_k^{*}{S}_k=A_k{\mathrm{CS}}_k={\Phi }_k^m;$

(3f)构造第k级网格上的内部靶向目标能量密度S_k的优化求解表达式，选取L₁范数正则化进行约束；

其中λ_k为正则化参数，为矩阵论中的L₂范数的平方，||·||₁为矩阵论中的L₁范数；

(3g)利用指数自适应迭代软阈值方法，求解(3f)中的优化表达式，得到内部 靶向目标能量密度S_k；利用所得内部靶向目标能量密度S_k，通过步骤(3e)中的线 性关系方程计算第k级离散网格上的光通量密度，记为判断所得到的是否满 足以下条件：其中ε为设定的优化停止条件，通常为一个极小的正数； 如果满足转步骤(3i)，上述所得的S_k即为内部靶向目标能量密度值，否则执行步骤 (3h)；

(3h)根据第k级网格上重建出的内部靶向目标能量密度S_k，采用移动网格以 及自适应网格细分策略，对k+1级的网格进行调整；然后在k+1级离散网格上确定 新的光源可行区，转步骤(3a)；

(3i)显示结果。

所述步骤(1)中将不同光谱段的切伦科夫发光光通量密度离散成更窄的波段上 的光通量密度，按如下步骤进行：

(1a)将CCD探测器获取的小动物体表N个连续波段的切伦科夫光学信号分为N 个波段τ₁，τ₂，...τ_N，其中τ_n＝[λ_n-1，λ_n]，n＝1，2，...，N，λ₀＝0，其中λ_n为波长，每个宽 波段τ_n测量的切伦科夫发光光通量密度为N为波段个数；

(1b)利用扫描光谱仪装置获取的光谱分布信息，将切伦科夫光学信号离散成比 τ_n更窄的M个波段κ₁，κ₂，...κ_M，同时保证每个窄波段κ_l只属于某一τ_n区间内，每个 τ_n区间内有多个窄波段κ_l，其中κ_l＝[λ_l-1，λ_l]，，l＝1，2，...，M，M为窄波段个数；

(1c)计算每一窄波段κ_l在各波段τ_n上所占的发光能量比重ω(κ_l)；

$\omega \left({\kappa }_l\right)=\frac{F\left({\kappa }_l\right)}{F\left({\tau }_n\right)};$

其中，$F\left({\tau }_n\right)={\int }_{{\lambda }_{n-1}}^{{\lambda }_n}f\left(\lambda \right)\mathrm{d\lambda },$$F\left({\kappa }_l\right)={\int }_{{\lambda }_{l-1}}^{{\lambda }_l}f\left(\lambda \right)\mathrm{d\lambda },$f(λ)为切伦科夫光学信号在波 长λ上的概率密度函数；

(1d)结合CCD探测器获取每个波段τ_n上的切伦科夫发光光通量密度，得到每 一窄波段κ_l的光通量密度

${\Phi }_k^m\left({\kappa }_l\right)={\Phi }_k^m\left({\tau }_n\right)·\omega \left({\kappa }_l\right);$

(1f)组合每一窄波段κ_l上的光通量密度，得到第k级离散网格上总的光通量密 度

${\Phi }_m^k=\left(\begin{array}{c}{\Phi }_k^m\left({\kappa }_1\right)\\ {\Phi }_k^m\left({\kappa }_2\right)\\ ···\\ {\Phi }_k^m\left({\kappa }_M\right)\end{array}\right).$

所述步骤(2)中的切伦科夫光学信号和小动物的三维体数据进行坐标匹配和能 量映射，是依据小动物支架上的球形标记点进行的，具体步骤如下：

(2a)将CCD探测到的带有标记点的二维切伦科夫光学信号与X光信号重建后得 到的带有标记点的三维体数据信息建立对应关系，经过坐标变换使两者位于同一三维 坐标系中，完成坐标匹配；

(2b)应用非接触式光学断层成像方法中的生物体表面三维能量映射技术，将坐 标匹配后的二维切伦科夫光学信号投影到小动物的三维表面上，获取小动物体表面的 三维光学数据分布信息；

本发明与现有技术相比具有如下优点：

第一，本发明将核素成像、光学成像和X光成像的优势进行了有效整合，开发了 一种集核素成像的高特异性、光学成像的高灵敏性和X光成像的高分辨率为一体的小 动物活体多模态分子成像系统，并实现了对单一标记探针的多模态同步成像，使获得 的小动物在体多模影像信息更加完整、全面和精确，解决了单一模态或多标记探针成 像不能避免的技术难题。

第二，本发明借助扫描式光谱仪装置获取切伦科夫光学信号在整个发光波长上的 光谱分布信息，有效计算了每一窄波段在各自宽波段上所占的发光能量密度的比重， 得到了每一窄波段上的光通量密度，为靶向目标的三维成像提供更多的信息；成像过 程中深入考虑小动物体内的靶向目标在不同深度情况下对体表测量信号的影响，构建 深度影响补偿矩阵，建立更切合实际的系统方程，能够实现对靶向目标更精确的成像。

第三，本发明根据生物体在解剖结构和组织光学参数方面的差异，建立基于RTE 高低阶近似理论的混合光传输数学模型，建立使用L1范数正则化约束的优化目标函 数，并使用指数自适应迭代软阈值方法进行求解，实现对复杂生物体内部靶向目标的 准确、快速成像。

附图说明

图1为本发明中小动物活体多模态分子成像系统的组成框图；

图2是本发明中小动物活体多模态分子成像系统示意图；

图3是本发明的小动物活体多模态分子成像主流程图；

图4是本发明对动物结构体数据分割后的示意图；

图5是本发明中小动物活体多模态分子成像方法的切伦科夫靶向目标三维重建子 流程图。

具体实施方式

下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。所描述的实例 仅旨在便于对本发明的理解。

参照图1和图2，本发明的小动物活体多模态分子成像系统包括：信号采集子系统 1和计算机子系统2，其中：

信号采集子系统1，包括切伦科夫发光成像装置3、核素成像装置4和MicroCT装 置5，旋转平移台11、电动平移及旋转台控制箱12、小动物支架13和暗箱14。该切伦 科夫发光成像装置3，包括CCD探测器31、CCD探测器电源箱32、滤波片组33和扫描 式光谱仪装置34；该核素成像装置4，包括两个PET探测41和探测器电源箱42；该 MicroCT装置5，包括X射线发射器51、发射器电源52、X射线CMOS平板探测器53。

上述各部件的结构关系如下：

切伦科夫发光成像装置3、核素成像装置4、MicroCT装置5以旋转台为中心安装 在实验台上，三套装置各组成部分在空间按照60度等间距均匀对称分布排列。切伦科 夫发光成像装置3中，滤波片组33位于CCD探测器31的正前方；其中CCD探测器31、 CCD探测器电源箱32和滤波片组33，三者与扫描式光谱仪34以旋转台为中心在空间呈 一条直线安装在实验台上。核素成像装置4中，一对PET探测器41及其探测器电源箱 42以旋转台为中心在空间呈一条直线安装在实验台上。MicroCT装置5中，X射线发射 器51和发射器电源52，与X射线平板探测器53以旋转台为中心在空间呈一条直线安装 在实验台上。暗箱14遮挡住三套成像装置的各组成部分，计算机位于暗箱的外部。暗 箱14根据成像系统的尺寸进行设计加工，尺寸为1800mm×1200mm×800mm，由黑 色吸光毡布、铅板、防氧化铝合金三层组成。暗箱14可从侧边打开，方便实验操作， 其功能是防止X射线对周围环境的污染，同时也起到保护装置部件和防尘的作用。

计算机子系统2，包括系统控制模块6、信号处理模块7以及数据重建模块8。该系 统控制模块6，实现对信号采集子系统1中各个成像系统的主要器件的控制，包括对电 动平移及旋转台11的平移和旋转控制，对CCD探测器31和PET探测器41的温度监控和 数据采集指令控制，对X射线发射器51的电压和电流的控制，对X射线平板探测器53 的数据采集指令控制，对滤光片组33的滤波轮旋转控制，对扫描式光谱仪装置34参数 设定和数据扫描的指令控制。该信号处理模块7，对信号采集子系统1中获取到的带有 各种噪声的信号，包括切伦科夫光学信号、核素信号和X光信号进行预处理，消除 信号采集过程CCD暗电流、直流偏置各种原因产生的噪声，并获取数据重建时的必要 参数，如球形标记点的坐标信息。该数据重建模块8，首先基于GPU加速的锥束反投 影算法，对预处理后的多角度核素信号和X光信号进行重建。将切伦科夫光学信号与 X光学信号重建后得到的三维体数据进行坐标匹配和能量映射，然后采用自适应双阈 值分割方法分割小动物的骨骼，采用区域生长法分割小动物的肝脏和肾脏，采用融合 先验知识的方法分割心脏、肝脏、胃等软组织器官，最终获得小动物的三维结构信息； 结合小动物三维结构信息和多光谱光学参数这些先验信息确定光源可行区域，进行切 伦科夫靶向目标的三维重建。

CCD探测器31和CCD探测器电源箱32之间由CCD电源电缆连接，CCD探测器31 与计算机子系统2之间由CCD数据电缆连接，滤波片组33与计算机子系统2之间由通过 滤波轮上USB控制电缆连接。CCD探测器31采用2048×2048的CCD阵列的科学级高 灵敏度相机，最低温度能达到-70℃，可以捕获生物体内微弱的生物发光信号，在450nm 至800nm的波长范围内量子效率高达75％以上。扫描式光谱仪装置34用于实现对目标 从可见到近红外光波段的光谱探测。PET探测41与计算机子系统2之间用数据电缆连 接，PET探测器41直径为11.8cm，分辨率为1.8mm。X射线发射器51与发射器电源52 之间由电源电缆连接，发射器电源52与计算机子系统2之间由电源控制电缆连接。X 射线发射器51的焦斑为35微米，X射线CMOS平板探测器53的成像区域尺寸为 12cm×12cm，像素大小为50μm，平板探测器53和计算机子系统2之间用数据电缆连 接。电动平移及旋转台11与电动平移及旋转台控制箱12之间通过串口控制电缆连接， 电动平移及旋转台控制箱12与计算机子系统2之间通过USB控制电缆连接。旋转台通 过步进电机驱动，最小分辨率为0.33″，最大中心负载为4kg，自重为3kg。平移台用 于调整固定在支架13上的活体小动物在暗箱14中的位置，最大行程为200mm，分辨 率为3.2μm，中心负载为30kg，自重5.5kg。

参照图3，本发明的小动物活体多模态分子成像方法，其实现步骤如下：

步骤1：数据采集。

在高灵敏制冷CCD探测器31前放置不同波段的带通滤光片组33，利用CCD探测器 采集连续光谱段的小动物体表切伦科夫光学信号，本例设连续光谱段为4个，但不局 限于4个，并采集带有球形标记点信息的白光图像以及背景图像；

对于多光谱信号的采集：使用CCD探测器获取到小动物体表N个连续波段的切伦 科夫发光信号；N个宽带连续滤波片将整个小动物体的发光信号分为N个波段 τ₁，τ₂，...τ_N，其中τ_n＝[λ_n-1，λ_n]，n＝1，2，...，N，λ₀＝0，其中λ_n为波长，每个宽波段τ_n测量的切伦科夫发光光通量密度为

利用光谱仪装置34对体表切伦科夫光学信号在500-800nm波段进行扫描探测，获 取其光谱分布信息；借助扫描式光谱仪装置获得的切伦科夫发光信号在整个发光波长 上的光谱分布信息，将光谱分布离散成比τ_n更窄的M个波段κ₁，κ₂，...κ_M，其中 κ_l＝[λ_l-1，λ_l]，l＝1，2，...，M；同时保证每个窄波段κ_l只属于某一τ_n区间内，每个宽τ_n区 间内可以有多个窄波段κ_l。

采用注射GE healthcare公司生产的医用碘海醇造影剂的方法来提高小动物体内 部软组织的对比度，并利用MircoCT装置5采集小动物活体多个角度的X光二维投影信 号，本例为360个角度，但不局限于360。

步骤2：数据预处理及配准。

(2a)利用计算机子系统2的信号处理模块7对采集到的切伦科夫光学信号和X光 信号进行如下预处理：

首先，对于切伦科夫光学信号而言，由于放射性同位素的韧致辐射信号对CCD 相机的噪声干扰很大，因而需要利用中值滤波方法对信号进行多次处理，以有效地降 低辐射噪声的影响；

其次，对于CCD暗电流、直流偏置各种原因产生的噪声，利用采集的背景信号对 目标信号进行统计分析，然后进行去噪；

最后，对探测器采集的二维投影信号进行暗场去除、坏点补偿以及亮场校正预处 理，以解决X射线探测器53中集成电路的制作工艺问题，

(2b)数据配准：

数据的坐标匹配主要是依据小动物支架13上的球形标记点进行的；本例中球形标 记点有4个，但不局限于4个；具体配准策略如下：

采用基于GPU加速的锥束反投影算法，利用计算机子系统2数据重建模块8，对360 个角度的X光二维投影信号进行重建，获取小动物体的表面轮廓和带有标记点信息的 三维体数据；

提取球形标记点在二维白光图像中的坐标信息，进而得到球形标记点在二维切伦 科夫光学图像中的坐标信息；由于切伦科夫光学信号与白光图像都是由CCD探测器31 拍摄到的二维图像，所以两者具有相同的坐标信息；

将带有球形标记点信息的二维切伦科夫光学信号与X光信号重建后得到的带有 标记点的三维体数据建立对应关系，经过坐标变换使两者位于同一三维坐标系中，完 成坐标匹配；

利用非接触式光学断层成像方法见ZL200910024292.1专利中描述的生物体表面 三维能量重建技术，将坐标匹配后的二维切伦科夫光学信号投影到小动物的三维表面 上，获取小动物体表面的三维切伦科夫发光信号分布。

步骤3：解剖结构数据获取。

数据采用前，采用注射医用碘海醇造影剂的方法来提高小动物体内部软组织的对 比度，通过锥束反投影算法重建后得到的CT体数据已经具备了一定的解剖结构信息， 但是其边界信息依然模糊，因此需要进行体数据分割。本发明将自动分割与手动分割 相结合，然后利用不同的分割方法完成对体数据的分割，获得全面的小动物解剖结构 信息。

对于小动物骨骼而言，其灰度范围与其他器官的灰度范围存在较大差异，因此采 集自适应双阈值分割方法：

对于肺而言，其内部与相邻的器官有较为明显的边界，可以将区域为处理对象， 采用区域生长法进行分割，再将分割后的肺部数据与高斯滤波器求卷积，然后取阈值 进行二值化处理；

对于小动物体的脂肪和肾脏也使用区域生长法来分割；

对于心脏、肝脏、胃这些软组织器官，其灰度值与相邻组织的差异较小，需要融 合待分割目标的先验知识进行分割：采用统计形状模型的方法充分利用特定器官的形 状和纹理的先验知识，以训练和统计分析为基础，采用主成分分析方法建立起关于目 标的平均形状和纹理及其变化模式的先验模型，在待分割数据中进行迭代搜索匹配， 使先验模型上的标记点逐步拟合到待分割器官的实际边界，同时手动对某些单张切片 的分割结果进行调整，最终准确分割出特定器官。

图4为本发明对动物结构体数据分割后的示意图，本例一共分割出7种器官：骨骼 401、肺402、脂肪403、肾脏404、心脏405、肝脏406、胃407。

步骤4：光学参数确定。

通过对小动物体CT体数据的分割，得到一个包括多种器官的非匀质生物组织模 型，利用信号处理模块7对该模型的表面和内部器官进行网格剖分，其中网格网格的 类型为四面体网格。假设每个器官内部的所有四面体单元在组织光学特征参数的分布 上都是一致的，通过参考G.Alexandrakis，F.R.Rannou，and A.F.Chatziioannou， “Tomographic bioluminescence imaging by use of a combined optical(OPET)system：a  computer simulation feasibility study，”Phys.Med.Biol.50，4225-4241，2005，可以确 定每种生物组织器官在500-800nm波长范围内的光学参数分布。

表1列出了上述生物组织器官在常用波长下的光学特性参数，其中，胃属于囊状 腔体，光在其内部传输时只发生光的吸收，几乎不发生散射；μ′_s是生物组织的约化 散射系数；μ_a是吸收系数，单位为mm^-1。

表1不同生物组织器官在不同波长下的组织光学特性参数

步骤5：辐射传输方程RTE的高低阶理论近似。

由于生物自发光信号和切伦科夫发光信号的光谱分布范围为可见光和近红外光 波段，故可使用辐射传输方程精确对其进行描述。但辐射传输方程又是一个复杂的积 分-微分方程，对于复杂模型很难直接求解，时间代价很高。以往的光学分子成像研 究针对绝大多数生物组织的高散射特性，采用辐射传输方程的一阶近似，即扩散方程 DE来描述光在生物组织中的传输过程。但是在小动物体中，不同的生物组织在不同 的波长下的传输特性不一样，不同生物组织在一定波段上会表现出高散射或者低散射 特性。对于切伦科夫发光成像，其光学信号具有较宽的光谱分布，并且在可见光波段 的能量分布较高。因此仅适用于高散射生物组织光传输特性的扩散方程在研究切伦科 夫发光成像中具有明显的局限性。从表1可以看出，发光光谱对应的不同器官的光学 特征参数波动范围很大，某些器官在某一波段上会表现出高吸收特性。

近年辐射传输方程的一些高阶近似形式，如球谐波近似方程、离散坐标近似方程、 相位近似方程等，能够准确地描述任意光学特性生物组织中的光传输过程，但计算复 杂度将远远超出现有的承受能力。为了平衡精度和求解速度，研究人员又提出了简化 球谐波近似方程SP_N来描述光传输过程。SP_N能够有效的求解光传输方程，同时能够 达到实际应用所要求的求解精度。虽然SP_N能够在保证求解精度的同时，提高求解效 率；但是与扩散方程DE相比，SP_N的求解效率仍是偏低。

为了最优地兼顾求解精度和求解效率，本发明采用联立SP_N和DE，建立基于RTE 高低阶近似理论的光传输混合数学模型，描述切伦科夫发光在复杂生物体中的传输过 程。在光传输混合数学模型的建立过程中，以生物组织光学特征参数随解剖结构的变 化和波长的波动为依据，将小动物体分为高散射特性区域和非高散射特性区域。DE 用于描述光在高散射特性组织区域的传输过程，SP_N用于描述光在非高散射特性组织 区域的传输过程。为了建立统一的光传输混合数学模型，需要在DE和SP_N之间建立边 界耦合条件，本发明用罗宾边界条件进行耦合

步骤6：切伦科夫靶向目标三维重建。

参照图5，本步骤的具体实现如下：

(6a)根据步骤(2b)得到的小动物体表的三维切伦科夫光学信号分布，划定先验 光源可行区；对求解区域进行网格离散化处理并初始化；

(6b)计算每一窄波段κ_l在各自光谱区间τ_n上所占的发光能量比重：

$\omega \left({\kappa }_l\right)=\frac{F\left({\kappa }_l\right)}{F\left({\tau }_n\right)},$l＝1，2，...，M，

其中，$F\left({\tau }_n\right)={\int }_{{\lambda }_{n-1}}^{{\lambda }_n}f\left(\lambda \right)\mathrm{d\lambda }$为波段τ_n上的发光能量，$F\left({\kappa }_l\right)={\int }_{{\lambda }_{l-1}}^{{\lambda }_l}f\left(\lambda \right)\mathrm{d\lambda }$为窄波段 κ_l上的发光能量，f(λ)为切伦科夫发光信号在波长λ上的光谱函数；

(6c)根据发光能量比重计算每一窄波段κ_l对应的切伦科夫发光光通量密度

${\Phi }_k^m\left({\kappa }_l\right)={\Phi }_k^m\left({\tau }_n\right)·\omega \left({\kappa }_l\right),$

其中为第k级离散网格上波段τ_n上的切伦科夫发光光通量密度；

(6d)基于hp混合自适应有限元策略，在第k级离散网格上，结合扩散方程DE和 简化球谐函数SP_N的光传输理论近似，建立在第k级网格某一波段κ_l上内部靶向目标 能量密度S_k与表面光通量密度结之间的线性关系方程：

$A_k\left({\kappa }_l\right)S_k={\Phi }_k^m\left({\kappa }_l\right),$

其中，A_k(κ_l)为S_k与的系统矩阵；

将上述M个窄波段的线性关系方程组合成一个整体，建立第k级网格上内部靶向 目标能量密度S_k和表面测得总的光通量密度之间的线性关系方程：

$A_k{S}_k={\Phi }_k^m,$

其中，$A_k=\left(\begin{array}{c}{A}_k\left({\kappa }_1\right)\\ A_k\left({\kappa }_2\right)\\ ···\\ A_k\left({\kappa }_M\right)\end{array}\right)$为总的系统矩阵，${\Phi }_m^k=\left(\begin{array}{c}{\Phi }_k^m\left({\kappa }_1\right)\\ {\Phi }_k^m\left({\kappa }_2\right)\\ ···\\ {\Phi }_k^m\left({\kappa }_M\right)\end{array}\right)$为表面测得总的光通量 密度；

(6e)构建深度影响补偿矩阵，修正线性关系方程

考虑到小动物体内深层的靶向目标对体表测量信号的影响，构建深度影响补偿矩 阵，校正深度变化对系统矩阵A_k的影响；使用新的系统矩阵重新建立内部靶向目标 能量密度S_k和表面测得总的光通量密度之间的线性关系方程，具体深度补偿策略 如下：

首先，采用动态分层策略将求解域进行圆柱体形式的深度分层，每层代表一个深 度，根据求解域的分层层数对系统矩阵A进行分块，从表面到深层共分成H块子系 统矩阵A₁，A₂，...，A_H，每个子系统矩阵都反映了该层深度值对表面测量的光通量密度 的影响；求得每一子系统矩阵对应的最大特征值s_l，1≤l≤H；

其次，定义第l层的深度补偿子矩阵为：

C_l＝s_H-l+1I；

其中，s_H-l+1为第H-l+1层深度上的子系统矩阵所对应的最大特征值，I为单位 矩阵；

然后，组合各深度补偿子矩阵C_l，获得总的深度影响补偿矩阵C，然后构造深 度补偿后的新系统矩阵

$A_k^{*}=A_{k}C;$

最后，利用深度补偿后系统矩阵建立内部靶向目标能量密度S_k和表面实际 测得的光通量密度之间的新的线性关系方程：

$A_k^{*}{S}_k=A_k{\mathrm{CS}}_k={\Phi }_k^m;$

(6f)构造第k级离散网格上的最优化目标函数：

其中残差项正则化项g(S_k)＝||S_k||₁，λ_k为正则化参数，即 拉格朗日乘子，函数Θ(S_k)＝f(S_k)+λ_kg(S_k)；为矩阵论中的L₂范数的平方，||·||₁为矩阵论中的L₁范数；

(6g)使用指数自适应迭代软阈值方法求解步骤(6f)中的第k级网格上的优化 目标函数；

该方法的基本思想是首先给定一个较大初始拉格朗日乘子，λ_k→+∞，求得的解 一定非常稀疏；在迭代过程中不断减小λ_k，使重建结果的非零值元素不断增加，进 而得到较好的重建结果；该方法的优点是当λ_k取值很大时函数Θ(S_k)的下降速度很 快，λ_k取值很小时函数Θ(S_k)的下降过程变慢，有利于算法迅速收敛到全局最优解； 其求解步骤如下：

(6g1)初始化：d＝1，λ_k→+∞，其中d为迭代次数，迭代开始

(6g2)采用如下指数函数对阈值进行收缩：

λ_k，d+1＝λ_k，d exp(-ρ_k*d)，

其中exp(·)为指数函数，ρ_k为下降速度因子，λ_k，d和λ_k，d+1分别为第d次和第d+1 次迭代所对应的拉格朗日乘子；

(6g3)利用上式所得λ_k，d+1，对靶向目标能量密度按照如下软阈值方法进行更新：

其中，S_k，d+1和S_k，d分别为第k级网格上第d次和第d+1次迭代所对应的靶向目标能 量密度，$u_{k,d}=S_{k,d}-\frac{1}{{\alpha }_{k,d}}▿f\left(S_{k,d}\right),$为求导函数；

${\alpha }_{k,d}=\frac{{(S_{k,d}-S_{k,d-1})}^T(▿f\left(S_{k,d}\right)-▿f\left(S_{k,d-1}\right))}{{(S_{k,d}-S_{k,d-1})}^T(S_{k,d}-S_{k,d-1})},$(·)^T为转置函数，S_k，d-1为第d-1次迭 代所对应的靶向目标能量密度；

软阈值函数soft(u，a)的表达式为：

$\mathrm{soft}(u,a)=\mathrm{sgn}\left(u\right)\max \left\{\right|u|-a,0\}$

$=\left(\begin{array}{cc}\mathrm{sgn}\left(u\right)\left(\right|u|-a)& \left|u\right|>a\\ 0& \left|u\right|\le a\end{array}\right),$

其中，sgn(·)为符号函数，max{e₁，e₂}为取e₁和e₂两者的最大值函数，|·|为绝对 值函数；

(6g4)利用步骤(6g3)求得的S_k，d+1，根据Θ(S_k，d+1)＝f(S_k，d+1)+λ_kg(S_k，d+1)来 计算函数Θ(S_k，d+1)的值；

(6g5)判断||Θ(S_k，d+1)-Θ(S_k，d)||＜ε_Θ或者d＞D_max是否满足，其中ε_Φ为设定的指 数自适应迭代软阈值方法停止条件，通常为一个极小的正数，D_max为自适应迭代软 阈值方法的最大迭代次数；如果满足，则停止迭代并执行步骤(6h)，否则令d＝d+1 并转步骤(6g2)；

(6h)利用步骤(6g)求解得到的内部靶向目标能量密度S_k，通过步骤(6e) 中的深度补偿后的线性关系方程，计算第k级离散网格上的光通量密度，记为

(6i)判断或者k＞L_max是否满足，其中ε为设定的优化停止条件， 通常为一个极小的正数，L_max为网格的最大级数；如果满足，上述所得的S_k即为内部 靶向目标能量密度，转步骤7，否则执行步骤(6j)；

(6j)根据第k级网格上重建出的光源分布S_k进行移动网格和自适应细分：

移动网格：将求解域中的部分节点移动到重建光源能量较大和能量变化较快的区 域，以减小该区域中解的误差；

自适应hp细分：设s_i为第i个四面体单元的靶向目标的能量密度，s_max为重建的靶 向目标的能量密度的最大值，如果满足计算式：s_i＞βs_max，其中β是一个常数且 0＜β＜1，则对该四面体单元进行p细分，并且对该四面体单元周围的四面体单元采 取h细分以避免产生悬点；

(6k)取k＝k+1，并在新的离散网格上确定光源可行区，然后转步骤(6d)， 在新的网格上进行内部靶向目标能量密度重建。

步骤7：结果显示

利用计算机子系统2的信号处理模块7对步骤6获得的靶向目标信息和步骤3获取 的动物体解剖结构信息进行融合，将重建的靶向目标在动物体中进行三维显示。