k) = D(k, i)，即
[0042] 其中
其中，第i个节点的 坐标为（Xi, yi, Zi)，第k个节点的坐标为（xk, yk, Zk)，Xii为第i个节点的能量，X \为第k个 节点的能量；
[0043] 步骤6. 2 :记初始参考度矩阵P为0,初始适选矩阵A为0,初始代表矩阵R为0,设 定聚类迭代次数；
[0044] 步骤6. 3 :计算代表矩阵R°kl和适选矩阵A °kl，计算公式如下：
[0046] 其中，R°kl(i，k)表示点k适合作为点i的类代表点的合适程度，A° kl(i，k)表示点 i适合作为点k的类代表点的合适程度，k'满足k'辛k，i'满足i'辛k ;
[0047] 步骤6. 4 :计算更新后的代表矩阵R和适选矩阵A，交替更新过程如下：
[0049] 其中，Rn+1表示第n+1次迭代得到的代表矩阵Rn+1° kl更新后的代表矩阵；Rn表示第 η次迭代得到的代表矩阵Rn°kl更新后的代表矩阵;A n+1表示第n+1次迭代得到的代表矩阵 An+1°kl更新后的代表矩阵;An表示第η次迭代得到的代表矩阵A n°kl更新后的代表矩阵；λ表 示阻尼系数；
[0050] 步骤6. 5 :判断聚类迭代次数是否达到设定的聚类迭代次数，如果达到，终止迭 代，输出代表矩阵R和适选矩阵Α，否则，转至步骤6. 3 ;
[0051] 步骤6. 6 :根据代表矩阵R和适选矩阵A对各个节点进行分类，并确定每一类的聚 类中心，具体方法如下：
[0052] 对于任意节点i，如果满足argkmax{A(i, k)+R(i, k)}，则节点i和节点k属于同一 类；对于每一类内的任意节点m，计算A(m, m)+R(m, m)，选取A(m, m)+R(m, m)中的最大值对应 的节点m，作为此类的聚类中心。
[0053] 具体地，所述步骤7的过程包括：
[0054] 步骤7. 1 :根据步骤3得到的网格节点和四面体之间的对应关系，计算所有四面体 的能量值，选取能量值非〇的四面体作为聚类对象，计算此聚类对象的几何中心；
[0055] 步骤7. 2 :依次计算每个能量值非0的四面体的几何中心与各初始聚类中心之间 的距离ΙΛ，…IJ，其中，c为初始聚类中心的数目，根据最小距离I 1将该四面体划分为第i 类；
[0056] 步骤7. 3 :重新计算每一类的聚类中心C :
[0058] 其中p (i)、a(i)、c (i)为每一类第i个非0四面体的能量、体积以及几何中心。
[0059] 与现有技术相比，本发明具有以下技术效果：
[0060] 1、本发明利用光学特性参数与解剖结构信息作为先验知识，结合多光谱荧光数 据，提高了 BLT光源重建的准确性以及重建图像的分辨率。
[0061] 2、本发明采用自适应可行区域迭代收缩策略，减少了 BLT重建问题的病态性，降 低了重建结果对于重建算法的依赖性。
[0062] 3、本发明将重建光源中心的确定问题转换为聚类分析问题，采用AP聚类和 k-means聚类算法，有效地利用了 AP聚类无需事先指定聚类数目的优点，以及k-means聚类 算法的相对可伸缩性，同时又避免了 AP聚类的不稳定性以及k-means聚类必须已知簇均 值的局限性；在混合聚类过程不仅直接使用了网格节点的能量分布信息，而且充分考虑了 节点所属四面体的重建能量、中心坐标、体积体等信息，在最大程度上提升了聚类算法的稳 定性，减少了计算时间复杂度。
【附图说明】
[0063] 图1为本发明的方法流程图；
[0064] 图2 (a)为仿真实验模型及中心相距IOmm的光源分布图；
[0065] 图2 (b)为仿真实验模型及中心相距5mm的光源分布图；
[0066] 图3 (a)为中心相距IOmm的探测器测量生物表面光强分布图；
[0067] 图3 (b)为中心相距5mm的探测器测量生物表面光强分布图；
[0068] 图4(a)为本发明在光源中心相距为IOmm情况下的重建光源分布图；
[0069] 图4(b)为本发明在光源中心相距为IOmm情况下的光源聚类中心分布图；
[0070] 图5(a)为本发明在光源中心相距为5mm情况下的重建光源分布图；
[0071] 图5(b)为本发明在光源中心相距为5_情况下的光源聚类中心分布图；
[0072] 下面结合附图和【具体实施方式】对本发明的方案做进一步详细地解释和说明。
【具体实施方式】
[0073] 实施例
[0074] 遵从上述技术方案，参见图1，本实施例的基于多光源分辨的生物发光断层成像重 建方法，具体包括以下步骤：
[0075] 步骤1 :获得成像目标的解剖结构信息、光学特异性参数以及表面多光谱荧光数 据
[0076] 针对皮下病变部位进行研究，成像目标定义为不考虑边界，且具有一定厚度的皮 肤组织，参照文献【1】，生成半径为42_，高度为24mm的圆柱体进行模拟。
[0077] 对成像目标进行CT扫描，得到其解剖结构信息，如形状、尺寸、位置等信息。
[0078] 将上述成像目标固定在电控旋转台上，将光学检测仪器放置在成像目标的上 方，成像目标内部生成的光源发出荧光，其波长由生物探针材料本身决定，其波长范围为 590~650nm，上述荧光在皮肤中的传输的光学特异性参数，参照文献【1】予以设置。
[0079] 荧光穿透成像目标，光学检测仪器检测穿透成像目标的多光谱荧光数据，上述荧 光数据并未采集成像目标侧面及底面的荧光数据，因此称为单角度数据采集。
[0080] 所述的光学检测仪器采用高灵敏度的C⑶相机。
[0081] 步骤2 :对步骤1得到的表面多光谱荧光数据进行归一化处理
[0082] 由于不同光谱的荧光吸收和散射性质不同，因此透射到成像目标表面的荧光光强 有很大的差异，为了消除这种差异对重建结果的影响，对每一波段的表面多光谱荧光数据 进行归一化处理：
[0084] 其中，Φ ( λ )为波长为λ的表面焚光数据，
i为归一化后的波长为λ的表面 荧光数据。
[0085] 步骤3 :基于光传输模型和有限元方法构建表面荧光数据与成像目标内部光源分 布的线性关系。具体实现方式如下：
[0086] 步骤3. 1 :根据步骤1得到的解剖结构信息和光学特异性参数得到光传输模型，所 述的光传输模型为扩散近似方程，利用扩散近似方程来描述光在成像目标内的传输过程。 [0087] 步骤3. 2 :利用nirfast工具对成像目标进行网格离散，获得成像目标的离散网格 数据，包括网格节点和四面体，得到网格节点和四面体之间的对应关系。
[0088] 步骤3. 3 :根据离散网格数据和有限元理论，对步骤3. 1中的扩散近似方程进行离 散，进而构建表面焚光数据与成像目标内部光源分布的线性关系：
[0089] Φ ( λ ) = M ( λ ) X
[0090] 其中，Μ(λ)是系统矩阵，表示成像目标内所有单元对光的吸收和散射作用，X是 要求解的成像目标内部光源能量分布，是非负的。
[0091] 步骤3. 4 :由于步骤2对不同光谱下的表面荧光数据做了归一化处理，因此对步骤 3. 3中的系统矩阵M做如下处理：
[0093] 则步骤3. 3中的线性关系也变为
为处理后的系统矩 阵。
[0094] 步骤4 :将步骤3. 4得到的线性关系转化为最优化问题：
[0096] 步骤5 :对步骤4中最优化问题