一种基于迭代自适应加权的有限视角光声成像重建方法.pdf

《一种基于迭代自适应加权的有限视角光声成像重建方法.pdf》由会员分享，可在线阅读，更多相关《一种基于迭代自适应加权的有限视角光声成像重建方法.pdf（10页完整版）》请在专利查询网上搜索。

1、(10)申请公布号 CN 103345770 A (43)申请公布日 2013.10.09 CN 103345770 A *CN103345770A* (21)申请号 201310302494.4 (22)申请日 2013.07.18 G06T 11/00(2006.01) A61B 5/00(2006.01) A61B 8/00(2006.01) (71)申请人 中国科学院自动化研究所 地址 100190 北京市海淀区中关村东路 95 号 (72)发明人 田捷 刘振宇 彭冬 马喜波 董迪 徐敏 (74)专利代理机构 中科专利商标代理有限责任 公司 11021 代理人 宋焰琴 (54) 发明名。

2、称 一种基于迭代自适应加权的有限视角光声成 像重建方法 (57) 摘要 本发明提供一种基于迭代自适应加权的有限 视角光声成像重建方法。基于有限视角采集的光 声信号， 为了校正计算投影信号和采集信号的差 异， 基于采集信号与上一步迭代重建图像的比值 对投影数据加权， 为了补偿重建图像和实际图像 之间的差异， 计算采集信号和计算机投影信号的 残差， 将残差信号反投影得到残差图像， 从而保证 迭代过程中信号残差逐步变小， 进而保证重建图 像的收敛性。 本发明能在有限视角扫描的情形下， 解决由于信号缺失引起的重建病态问题， 减少重 建伪迹， 提高重建精度， 对减少信号采集时间具有 一定的指导意义。 (。

3、51)Int.Cl. 权利要求书 2 页 说明书 5 页 附图 2 页 (19)中华人民共和国国家知识产权局 (12)发明专利申请 权利要求书2页 说明书5页 附图2页 (10)申请公布号 CN 103345770 A CN 103345770 A *CN103345770A* 1/2 页 2 1. 一种基于迭代自适应加权的有限视角光声成像重建方法， 其特征在于包括以下步 骤 ： 步骤 1， 重建图像初始化 ； 步骤 2， 计算加权系数 ； 步骤 3， 计算模拟信号 ； 步骤 4， 计算信号残差 ； 步骤 5， 计算重建图像 ； 步骤 6， 计算模型误差， 设置迭代终止条件。 2.根据权利要求。

4、1所述的方法， 所述步骤1包括采集光声信号的预处理， 重建图像的离 散化， 以及计算信号与图像像素之间的关系。 3. 根据权利要求 2 所述的方法， 其中采集光声信号的预处理包括对采样数据的选取和 滤波， 预处理后的采集光声信号 p(r0,t) 的总数表示为 M=mk， m 为投影位置个数， k 为每 个位置的数据长度， 滤波完成对采集信号中高频噪声的预处理。 4. 根据权利要求 2 所述的方法， 其中重建图像的离散化具体为 ： 在重建区域上叠加一个 NN 方格网， 每一个网格内像素值相同， 用 A(n)表示第 n 像素 的像素值。 5. 根据权利要求 2 所述的方法， 其中计算信号与图像像素。

5、之间的关系包括投影矩阵和 反投影矩阵的计算 ： 投影矩阵用来表示图像像素值A(n)对计算机模拟光声信号p的贡献， 用p(j)表示第j个 计算机模拟光声信号， 计算关系如下 ： 其中 wjn表示第 n 像素值 A(n)对第 j 个计算机模拟光声信号 p(j)的贡献值， wjn按以下 方式计算 ： 把投影弧线看成宽度为 ， 间隔为 0 的弧线束， wjn表示第 n 像素与第 j 条弧线 相交的面积的近似值， 将 wjn按信号采集顺序排列组成离散投影矩阵 P ； 反投影矩阵用来表示计算机模拟信号 p 对重建图像像素值 A(n)的贡献值， 计算关系如 下， 对于任意的 1 j M 有 ： 其中， Qn。

6、j表示反投影矩阵 Q 中第 n 行的第 j 个元素， j 是信号个数索引指标， n 是图像 像素索引指标，(xn,yn) 是重建光声图像的像素坐标， (x0,y0) 是 信号采样位置坐标， c 是声速， f 是信号采集的频率， mod 表示取余数运算。 6. 根据权利要求 5 所述的方法， 对于二维的有限视角光声成像， 重建图像 A(r) 和光声 信号 p(r0,t) 的关系近似为 ： 权 利 要 求 书 CN 103345770 A 2 2/2 页 3 其中 1为探测器的最小接收角度， 2为探测器的最大接收角度， e为有效扫描角 度， 表示每个像素点与信号采集位置最大接收角度和最小接收角度的。

7、夹角 ； 根据采集信号 p(r0,t) 对图像每个像素的贡献值， 利用上式计算得到重建光声图像的 初始值 A0(r)。 7. 根据权利要求 1 所述的方法， 所述步骤 2 为采用重建光声图像初始值 A0(r) 和处理 信号 p(r0,t) 利用下式计算校正系数 ： 8. 根据权利要求 1 所述的方法， 所述步骤 3 为基于离散化的投影矩阵 P 和初始光声图 像 A0(r)， 计算计算机模拟光声信号 p1， p1=PA0(r)。 9. 根据权利要求 1 所述的方法， 所述步骤 4 为利用校正系数 1和模拟光声信号 p1， 得 到和采集光声信号 p(r0,t) 同一量级的残差信号 p1=p(r0,。

8、t)-1p1。 10. 根据权利要求 1 所述的方法， 所述步骤 5 为将残差信号 p1反投影可以得到残差 校正图像 ： A1(r)=Qp1 将残差校正图像 A1(r) 按一定比例 (0,1) 加到初始重建图像 A0(r)， 得到更新 图像， 即第一步迭代重建图像 A1(r) A1(r)=A0(r)+A1(r)， 其中 为残差图像的校正系数 ； 重复步骤 2-4 及步骤 5 中的上述步骤， 得到迭代重建图像 Ai(r)=Ai-1(r)+Ai(r)=Ai-1(r)+Q(p(r0,t)-ipi)i=1,2,。 11. 根据权利要求 1 所述的方法， 所述步骤 6 具体为计算重建误差 ei ei=|。

9、Ai(r)-Ai-1(r)|2/|Ai-1(r)|2； 设置最小迭代误差和 / 或最大迭代次数作为迭代终止条件， 具体为 ： 先判断重建误差 是否小于设定的最小迭代误差， 如果达到设定的最小迭代误差， 则终止迭代 ； 若没有达到设 定的最小迭代误差， 再判断迭代次数是否达到最大迭代次数， 若达到最大迭代次数， 则迭代 终止。 权 利 要 求 书 CN 103345770 A 3 1/5 页 4 一种基于迭代自适应加权的有限视角光声成像重建方法 技术领域 0001 本发明涉及一种光声成像 （Photoacoustic Tomography， 简称PAT） 技术， 具体涉及 迭代自适应加权的有限视。

10、角光声图像重建方法。 背景技术 0002 光声成像是近十年来新兴的生物医学成像技术， 它结合了光学成像高对比度和超 声成像的高分辨率的特点， 近年来获得飞速的发展。光声图像不仅能够提供组织的结构形 态， 还能反映生理代谢等功能信息， 光声成像技术在生物组织成像中已经得到广泛的应用， 如肿瘤检测、 血管成像、 脑功能成像等相关研究。 0003 短脉冲激光照射到成像样品， 样品组织局部吸收光能而产生热膨胀向四周辐射超 声波， 通过超声探测器在不同位置的采集光声信号， 使用图像重建方法算出组织的吸收强 度分布。 图像重建是由信号转化成图像的关键， 而成像方法又是图像重建的核心技术， 目前 针对规则扫。

11、描形状已经提出了重建方法， 包括滤波反投影法， 时域重建法， 频域重建法， 延 迟求和法等。 0004 上述图像重建方法都在实际应用中， 都需要全方位扫描获取完备数据， 数据采集 速度慢， 对系统要求高。 本发明提出的图像重建方法， 对有限视角扫描的欠采样信息也能精 确地重建出光声图像。 发明内容 0005 （一） 发明目的 0006 本发明目的在于克服现有图像重建技术数据冗余的缺点， 提供一种针对有限视角 信息的光声成像方法， 基于有限角度少量扫描的欠采样数据和迭代自适应加权方法， 精确 重建生物组织吸收强度分布图像， 解决了生物组织无法全方位扫描的成像问题。 0007 （二） 技术方案 0。

12、008 本发明提出了一种基于迭代自适应加权的有限视角光声成像重建方法， 具体步骤 如下 ： 首先利用有限视角滤波反投影方法得到迭代初始图像， 基于初始图像和采集光声信 号计算投影信号的加权系数， 利用初始图像和加权系数获取计算机模拟光声信号， 计算模 拟光声信号与采集光声信号的残差， 对残差信号反投影获取残差图像， 将残差图像加到上 一次迭代生成图像获得修正重建图像， 重复上述过程， 设置最大迭代次数和迭代误差作为 迭代结束条件。 0009 根据本发明的实施例， 提供了一种基于迭代自适应加权的有限视角光声成像重建 方法， 包括以下步骤 ： 0010 步骤 1， 重建图像初始化 ； 0011 步。

13、骤 2， 计算加权系数 ； 0012 步骤 3， 计算模拟信号 ； 0013 步骤 4， 计算信号残差 ； 说 明 书 CN 103345770 A 4 2/5 页 5 0014 步骤 5， 计算重建图像 ； 0015 步骤 6， 计算模型误差， 设置迭代终止条件。 0016 根据本发明的一个实施例， 重建图像初始化包括采集光声信号的预处理， 重建图 像的离散化， 以及计算信号与图像像素之间的关系。 0017 根据本发明的一个实施例， 计算信号与图像像素之间的关系包括投影矩阵和反投 影矩阵的计算。 附图说明 0018 图 1 为本发明的基于迭代自适应加权的有限视角光声成像重建方法流程图 ； 0。

14、019 图 2 为计算机模拟投影信号的几何分布 ； 0020 图 3a、 3b 和图 3c 为基于本发明方法的计算机仿真数据重建结果， 其中 ； 其中图 3a 是用 90范围内 30 个探测信号重建的光声图像 ； 图 3b 是用 135范围内 45 个探测信号重 建的光声图像 ； 图 3c 是用 180范围内 60 个探测信号重建的光声图像 ； 0021 图 4a 和图 4b 为基于本发明方法的离体生物组织的重建结果， 其中图 4a 为离体组 织仿体照片， 图 4b 本发明方法重建图像。 具体实施方式 0022 用短脉冲激光均匀照射成像样品， 超声探测器在成像平面内扫描， 接收样品产生 的光声。

15、信号， 光声信号的产生服从的光声成像方程为 ： 0023 0024 其中代表汉密尔顿算子， c 表示生物组织中声波传播的速度， p(r0,t) 表示 r0位 置处t时刻的光声信号， 表示声压膨胀系数， Cp为比热系数， A(r)是组织光能量吸收分布 图像， I(t) 是激光脉冲函数。光声图像重建是典型的逆问题， 即如何由 p(r0,t) 求出 A(r)。 0025 基于迭代自适应加权的有限视角光声成像重建方法如下 ： 0026 第一步， 重建图像初始化， 包括采集光声信号的预处理和重建图像的离散化。 采集 光声信号的预处理是指采样数据的选取和滤波， 根据成像区域的大小和位置在采集信号pin 中。

16、选取对应的信号， 对选取信号进行Ramp滤波得到光声信号p(r0,t)， 完成对采集信号中高 频噪声的预处理。 0027 重建图像离散化即在重建区域上叠加一个 NN 方格网， 每一个网格内像素值相 同， 用 A(n)表示第 n 像素的像素值， 像素总数为 N2； 用 p(j)表示第 j 个计算机模拟光声信号， 计算机模拟光声信号的总数与预处理后的采集光声信号相同， 表示为 M=mk， 它是投影位 置个数 m 和每个位置数据长度 k 的乘积。如图 2 所示， A(n)和 p(j)的关系可以表示如下 0028 0029 式中wjn表示第n像素值A(n)对第j个计算机模拟光声信号p(j)的贡献， 这。

17、里wjn可 以按以下方式计算 ： 把投影弧线看成宽度为 ， 间隔为 0 的弧线束， wjn表示第 n 像素与第 j 条弧线相交的面积的近似值， 将 wjn按信号采集顺序排列组成离散投影矩阵 P。 说 明 书 CN 103345770 A 5 3/5 页 6 0030 对于二维的光声成像， 通过引入有效扫描角的概念， 改进传统的滤波反投影算法， 得到有限视角的光声成像图像 ： 0031 0032 其中 1为探测器的最小接收角度， 2为探测器的最大接收角度， e为有效扫描 角度， 其定义是每个像素点与信号采集位置最大接收角度和最小接收角度的夹角， 为计算 每个光声信号对图像像素的贡献， 我们给出如。

18、下定义 ： 0033 A(r)=Qp 0034 其中 A(r) 表示重建图像， 大小为 N21， p 表示光声信号， 大小为 M1， Q 是表示信 号对图像像素的贡献值的反投影矩阵， 大小为 N2M， 对于每一个像素值 A(n)， 计算模拟光声 信号 p 对其的贡献， 对于任意的 1 j M 有 ： 0035 0036 其中， Qnj表示反投影矩阵Q中第n行的第j个元素， j是信号个数索引指标， n是图 像像素索引指标，(xn,yn) 是重建光声图像的像素坐标， (x0,y0) 是信号采样位置坐标， c 是声速， f 是信号采集的频率， mod 表示取余数运算。 0037 根据有限视角采集的光。

19、声信号， 给定信号采集范围最大接受角度 2和最小接受 角度 1， 计算每个像素值对应的有效视角， 通过傅里叶变换求得光声信号的求导运算值， 然后根据每个光声信号对图像每个像素的贡献值， 利用式 （3） 对非零贡献值进行积分运算， 可以得到迭代重建光声图像的初始值 A0(r)。 0038 第二步， 加权系数的计算， 用重建光声图像初始值 A0(r) 和处理信号 p(r0,t) 计算 校正系数 ： 0039 0040 第三步， 模拟信号的计算， 基于离散化的投影矩阵P和初始光声图像A0(r)， 我们可 以计算计算机模拟光声信号 p1， 0041 p1=PA0(r) 0042 第四步， 信号残差的计。

20、算， 为了校正计算机模拟信号p1和采集光声信号p(r0,t)的 差异， 利用用校正系数 1和模拟光声信号 p1， 得到和采集光声信号采集信号同一量级的残 差信号 0043 p1=p(r0,t)-1p1 (5) 0044 第五步， 重建图像的计算， 将残差信号 p1反投影可以得到残差校正图像 ： 0045 A1(r)=Qp1 (6) 将残差校正图像 A1(r) 按一定比例加到初始重建图像 说 明 书 CN 103345770 A 6 4/5 页 7 A0(r)， 可以得到更新图像， 即第一步迭代重建图像 A1(r) 0046 A1(r)=A0(r)+A1(r) (7) 其中 (0,1) 为残差图。

21、像的校正系数， 常根据 重建图像效果的经验值设定， 目的是防止过度校正。 0047 重复上述步骤二至四及步骤五的上述步骤， 可以得到迭代重建图像 0048 Ai(r)=Ai-1(r)+Ai(r)=Ai-1(r)+Q(p(r0,t)-ipi)i=1,2, (8) 0049 第六步， 模型误差的计算， 因为不知道初始图像， 我们选择相对的重建误差 ei作为 判定迭代终止的标准 ： 0050 ei=|Ai(r)-Ai-1(r)|2/|Ai-1(r)|2(i=1,2,n) (9) 其中 为设定好的 最小迭代误差。 0051 设置最小迭代误差和 / 或最大迭代次数作为迭代终止条件， 具体为 ： 先判断重。

22、建 误差是否小于设定的最小迭代误差， 如果达到设定的最小迭代误差， 则终止迭代 ； 若没有达 到设定的最小迭代误差， 再判断迭代次数是否达到最大迭代次数， 若达到最大迭代次数， 则 迭代终止。满足迭代终止条件则输出重建结果， 通过上述迭代可以重建光声图像。发明人 为验证本发明方法的有效性进行了如下实验 ： 0052 计算机仿真实验 ： 给定已知的光吸收分布图， 在 180弧形扫描情况下采集光声 信号， 步长为 3， 根据有限视角光声方程计算仿真信号， 用仿真光声信号重建光吸收分布 情况。 0053 离体组织实验 ： 在准备好的脂肪上凿三个深 1mm、 直径 1mm 的小孔， 小孔之间的距 离约。

23、为 2mm， 然后将准备好的猪肝样品填到脂肪中， 最后再覆上一层厚度约为 3mm 的脂肪， 成为一体后进行光声信号采集。 0054 实验中采用的光源是Nd ： YAG脉冲激光器 （LS-2134， LOTIS） ， 激光的波长是532nm， 脉宽是 7ns， 重复频率是 10Hz， 脉冲激光经凹透镜扩束、 毛玻璃均匀后照射到样品上， 单脉 冲能量控制在 20mJ/cm2以下。采用中心频率 5MHz， 直径为 13mm 的美国 OLYMPUS 公司生 产的超声探头采集信号， 超声探测器和样品都放在耦合液里， 目的是更好地耦合超声信 号， 探测器在一个位置采集完信号后， 探头的旋转实现由重复定位精。

24、度为 0.005 度， 分辨 率为 0.00125 度的 ERSP100 旋转台带动实现， 旋转半径为 45mm， 扫描间隔是 3， 探测器 采集信号经过前置放大器预处理， 包括超低噪声时间增益补偿、 模数转换等， 信号通过一 个采样频率为 100MHz 的数字示波器 (MSO4000B,Tektronix) 进行采集。探测器在一个 位置采集 50 个信号后旋转到下一个位置， 一共需要旋转 60 次， 信号传输到计算机后， 用 MATLAB(version7.6,Mathworks) 对信号处理， 然后用基于迭代自适应加权的有限视角光声 成像重建方法进行图像重建。 0055 图3a、 3b和3。

25、c显示了基于本发明方法的计算机仿真数据的重建结果， 像素个数为 256256。其中图 3a 显示了用 90范围内 30 个探测信号重建的光声图像， 图 3b 显示了 用 135范围内 45 个探测信号重建的光声图像， 图 3c 显示了用 180范围内 60 个探测信 号重建的光声图像。 0056 图 4（a） 显示了离体生物组织的照片， 图 4（b） 显示了用 180范围内 60 个探测 信号重建的图像。 0057 由上述 2 组实验结果表明， 本发明的重建方法在 180信号下重建图像和原始图 像基本一致， 说明本发明能在有限视角扫描信号下， 精确地重建出组织的光吸收分布， 对减 说 明 书 。

26、CN 103345770 A 7 5/5 页 8 少信号采集时间具有一定的实际意义。 0058 以上所述的具体实施例， 对本发明的目的、 技术方案和有益效果进行了进一步详 细说明， 所应理解的是， 以上所述仅为本发明的具体实施例而已， 并不用于限制本发明， 凡 在本发明的精神和原则之内， 所做的任何修改、 等同替换、 改进等， 均应包含在本发明的保 护范围之内。 说 明 书 CN 103345770 A 8 1/2 页 9 图 1 说 明 书 附 图 CN 103345770 A 9 2/2 页 10 图 2 图 3a 图 3b图 3c 图 4a 图 4b 说 明 书 附 图 CN 103345770 A 10 。