，采用自适应可行区域迭代收缩策略进行代数迭代 求解，得到最优解矿。具体实现方法如下：
[0097] 步骤5. 1 :初始光源可行区域R设置为整个成像目标，并且设定区域收缩因子
其中N1为整个成像目标内包含的节点，心为迭代终止时非0解的数目，设 为1，Lnax为区域迭代收缩的次数，设为50。
[0098] 步骤5. 2 :结合可行区域R，步骤4的最优化问题变为
其 中，
为可行区域R下的系统矩阵，具体求解采用代数迭代算法（ART)，得到可行区域 R下的最优解X(R)，其反应可行区域R下的每个节点的能量大小。
[0099] 步骤5. 3 :根据步骤5. 2求得的可行区域R下的最优解X (R)，以及此可 行区域下的系统矩阵
计算成像目标表面的光通量Φ1λ)。计算光学检 测仪器检测到的穿透成像目标的多光谱荧光数据
与光通量之差，即
其中，i表示迭代次数（0 < i < Lmax) 〇
[0100] 步骤5. 4 :将步骤5. 3得到的f⑴与min(f (1)，…，f (i-Ι))进行比较，如果f(i) 小于min (f (1)，…，f (i-1))，则将全局最优解矿(R)用X (R)进行替换，否则，不作替换。。
[0101] 步骤5. 5 :将步骤5. 2求得的最优解X(R)从大到小进行排序，选取前Iength(R)/ β个值对应的节点来构建新的可行区域R。
[0102] 步骤5. 6 :判断迭代次数i是否达到最大值Lniax，如果达到，终止区域迭代收缩，输 出全局最优解Xi(R);如果没有达到，转至步骤5. 2。
[0103] 步骤5.7 :由全局最优解Xi(R)得到整个成像目标内的最优解X%若节点位于可行 区域R内，Xi= Xi(R);若节点位于可行区域R外部，Xi= 0。
[0104] 最优解Xi表示整个成像目标内各个节点的能量分布，其展示在空间分布图上，即 实现成像目标内部生成的光源的重建。
[0105] 为了准确获取成像目标内部生成的光源的中心位置，本发明的方法还包括以下步 骤：
[0106] 步骤6 :将步骤5得到的最优解矿进行AP中心聚类，得到每一类的聚类中心。由 于Xi每个元素对应的是网格节点的能量，在中心聚类时，选取X""中的非0元素对应的节点， 按照其空间分布进行聚类，由于AP聚类可以自动确定聚类数目和聚类中心，因此采用AP聚 类算法在确定光源数目的基础上，初步对各个光源中心定位。
[0107] 具体实现方法如下：
[0108] 步骤6. 1 :选取Xi中的非0元素对应的节点构建相似度矩阵D，其中，D为对称矩 阵 D(i, k) = D(k, i)，即
[oho] 其中
其中，第i个节点的 坐标为（Xi, yi, Zi)，第k个节点的坐标为（xk, yk, Zk)，Xii为第i个节点的能量，X \为第k个 节点的能量。
[0111] 步骤6. 2 :记初始参考度矩阵P为0,初始适选矩阵A为0,初始代表矩阵R为0,聚 类迭代次数为1000。
[0112] 步骤6. 3 :AP聚类算法传递两种类型的消息矩阵：代表矩阵R°kl和适选矩阵A°k1，参 考文献【2】，计算公式如下：
[0114] 其中，R°kl(i，k)表示点k适合作为点i的类代表点的合适程度。
表示点 i适合作为点k的类代表点的合适程度。k'满足k'辛k，i'满足i'辛k。
[0115] 步骤6. 4 :AP聚类算法的迭代过程就是这两个消息矩阵R和A交替更新的过程，设 阻尼系数λ为0.8,计算更新后的代表矩阵R和适选矩阵A，交替更新过程如下：
[0117] 其中，Rn+1表示第η+1次迭代得到的代表矩阵Rn+1° kl更新后的代表矩阵；Rn表示第 η次迭代得到的代表矩阵Rn°kl更新后的代表矩阵;A n+1表示第n+1次迭代得到的代表矩阵 An+1°kl更新后的代表矩阵;A n表示第η次迭代得到的代表矩阵A n°kl更新后的代表矩阵；
[0118] 步骤6.5 :判断聚类迭代次数是否达到设定的聚类迭代次数，如果达到，终止迭 代，输出代表矩阵R和适选矩阵A，否则，转至步骤6. 3。
[0119] 步骤6. 6 :根据代表矩阵R和适选矩阵A确定聚类中心m
[0120] 对于任意节点i，如果满足
则节点i和节点k属于同一 类；对于每一类内的任意节点m，计算A(m, m)+R(m, m)，选取A(m, m)+R(m, m)中的最大值对应 的节点m，作为此类的聚类中心。
[0121] 步骤7 :对步骤6中AP聚类得到的聚类中心，作为k-means聚类算法的初始聚类 中心，进行k-means聚类，得到每一类的聚类中心C。由于最优解的稀疏性导致AP聚类结果 的稳定性较差，因此k-means聚类依据能量非0四面体的几何中心进行聚类。如具体实现 方法如下：
[0122] 步骤7. 1 :根据步骤3得到的网格节点和四面体之间的对应关系，计算所有四面体 的能量值，选取能量值非〇的四面体作为聚类对象，计算此聚类对象的几何中心；上述设计 可充分利用网格节点和四面体之间的相关性，获取更多的非零数据进行聚类，最大程度上 降低了由于最优解矿的稀疏性导致的偶然误差，采用四面体空间分布进行聚类，增加了聚 类数据源的可靠性与稳定性，从而提高了聚类分析的准确度进而降低了试验误差。
[0123] 步骤7. 2 :依次计算每个能量值非0的四面体的几何中心与各初始聚类中心之间 的距离ΙΛ，…IJ，其中，c为初始聚类中心的数目，根据最小距离I 1将该四面体划分为第i 类。
[0124] 步骤7. 3 :重新计算每一类的聚类中心C :
[0126] 其中p (i)、a(i)、c⑴为每一类第i个非0四面体的能量、体积以及几何中心。
[0127] 步骤8 :利用matlab相关工具包，以及nirfast、iso2mesh相关函数将步骤5得到 的最优解，进行几何3D展示，实现成像目标内部光源的重建结果的显示；对步骤7得到的 聚类中心C进行几何3D展示，实现成像目标内部光源中心位置的确定。
[0128] 实验结果分析
[0129] 结合附图2、附图3、附图4对本发明的重建结果和聚类中心做进一步的描述。
[0130] 用于实验的仿体如附图2所示。其为半径为42mm，高为24mm的圆柱，主要模拟深 层皮下组织的病变。在试验中设置了中心相距IOmm和5_的两组仿体实验，其中图2(a) 模拟5个中心相距为IOmm的光源，图2 (b)模拟5个中心相距为5mm的光源。
[0131] 重建过程所需的表面光强分布图如附图3所示。由于模拟皮肤组织，所以图3的表 面光强分布均在仿体上表面采集。图3(a)是5个中心相距为IOmm的光源表面测量数据。 图3(b)是5个中心相距为5mm的光源表面测量数据。对比可知中心相距距离的减少对于 重建难度的影响。
[0132] 最终的多光源重建结果如附图4和附图5所示。其中，附图4为光源中心相距 IOmm情况下，基于本发明的实验重建结果和聚类结果。附图5为光源中心相距5mm情况下， 基于本发明的实验重建结果和聚类结果。其中图4(a)和图5 (a)为重建结果展示图，红色 区域为重建的光源分布，蓝色的圆圈代表真实的光源中心，蓝色的十字代表聚类结果中心。 图4(b)和图5(b)展示了最终的聚类结果图，中心点为最终聚类中心。在临床应用中，由图 4(a)、图5(a)和图4(b)、图5(b)可以清楚的得到病变区域，并准确确定区域中心。基于本 发明的重建方法，重建光源位置误