《基于正则化迭代的有限角度光声成像重建方法及装置.pdf》由会员分享，可在线阅读，更多相关《基于正则化迭代的有限角度光声成像重建方法及装置.pdf（10页完整版）》请在专利查询网上搜索。
1、(10)申请公布号 CN 103310472 A (43)申请公布日 2013.09.18 CN 103310472 A *CN103310472A* (21)申请号 201310251165.1 (22)申请日 2013.06.21 G06T 11/00(2006.01) A61B 5/00(2006.01) A61B 8/00(2006.01) (71)申请人 中国科学院自动化研究所 地址 100190 北京市海淀区中关村东路 95 号 (72)发明人 马喜波 田捷 刘学彦 杨鑫 (74)专利代理机构 中科专利商标代理有限责任 公司 11021 代理人 宋焰琴 (54) 发明名称 基于正则。
2、化迭代的有限角度光声成像重建方 法及装置 (57) 摘要 本发明公开了一种基于正则化迭代的有限角 度扫描光声成像的重建方法及装置。基于有限角 度的光声信号， 本发明在每次迭代过程实现残差 更新图像和正则化更新图像， 计算采集信号和重 建图像的计算机模拟信号的残差， 通过将残差信 号反投影得到残差图像， 将残差图像叠加到上一 次重建图像得到更新图像， 然后利用局部正则化 更新重建图像， 结合信号残差反投影法和正则化 方法获得重建图像， 本发明能在有限角度扫描情 形下， 快速精确地重建光声图像， 对减少重建时 间、 降低设备硬件成本有一定意义。 (51)Int.Cl. 权利要求书 3 页 说明书 。
3、5 页 附图 1 页 (19)中华人民共和国国家知识产权局 (12)发明专利申请 权利要求书3页 说明书5页 附图1页 (10)申请公布号 CN 103310472 A CN 103310472 A *CN103310472A* 1/3 页 2 1. 一种基于正则化迭代的有限角度扫描光声成像重建方法， 其包括 ： 步骤 1、 利用有限角度滤波反投影法得到初始的重建光声图像 A0(r) ； 步骤 2、 利用上轮迭代获得的重建光声图像计算投影加权系数， 对于第一轮迭代， 上轮 迭代获得的重建光声图像为初始的重建光声图像 A0(r) ； 步骤 3、 利用上轮迭代获得的重建光声图像和投影加权系数获取计。
4、算机模拟投影信 号 ； 步骤 4、 计算所述计算机模拟投影信号与采集信号之间的信号残差， 并根据所述信号残 差获取残差图像并修正重建光声图像， 得到更新后的重建光声图像 ； 步骤 5、 利用正则化运算对更新后的重建光声图像进行二次优化， 得到本轮迭代的重建 光声图像， 并转步骤 2 进行下轮迭代， 直至迭代完成。 2. 如权利要求 1 所述的光声成像重建方法， 其特征在于， 所述初始的重建光声图像如 下获得 ： 其中， A(r) 表示生物组织内部的光吸收分布， p(r0， t) 表示信号采集位置 r0处 t 时刻 的光声信号， c 表示生物组织中声波传播的速度， 表示声压膨胀系数， Cp为比热。
5、系数， 1 为探测位置的最小接收角度， 2为探测位置的最大接收角度， e为有效扫描角度， 对于图 像中的每个像素， 其定义为像素与最大接受位置点之间线段和像素与最小接受位置点之间 线段的夹角。 3. 如权利要求 1 所述的光声成像重建方法， 其特征在于， 所述投影加权系数如下计算 ： 其中， i表示投影加权系数， Ai-1(r) 表示上次迭代得到的生物组织内部信号产生位置 r 处的重建光声图像， p(r0， t) 表示信号采集位置 r0处 t 时刻的光声信号。 4. 如权利要求 1 所述的光声成像重建方法， 其特征在于， 所述计算机模拟投影信号如 下计算 ： iPAi-1(r) 其中， i表示。
6、投影加权系数， Ai-1(r) 表示上次迭代得到的生物组织内部信号产生位置 r 处的重建光声图像， P 表示离散化投影矩阵。 5. 如权利要求 4 所述的光声成像重建方法， 其特征在于， 投影矩阵 P 如下计算 ： 其中， Pmn 表示投影矩阵 P 中第 m 行的第 n 个元素， 投影矩阵 P 的大小为 MN2， M 为光 声信号向量元素的个数， 大小为采样位置个数和每个位置信号长度k的乘积， N2为重建光声 图像上所有像素点个数 ；(xn， yn) 是光声图像像素坐标， (x0， y0) 权 利 要 求 书 CN 103310472 A 2 2/3 页 3 是信号采集位置坐标， c 是声速，。
7、 dt 是信号采集的时间间隔， m 是光声信号索引指标， n 是图 像像素索引指标， Smn表示第 n 个像素与第 m 条投影弧线相交面积的近似值， mod 是指取 余数运算。 6.权利要求1所述的光声成像重建方法， 其特征在于， 步骤4中所述信号残差通过实际 采集到的信号 p(r0， t) 减去所述计算机模拟投影信号 iPAi-1(r) 获得。 7. 如权利要求 1 所述的光声成像重建方法， 其特征在于， 步骤 4 中， 更新后的重建光声 图像如下表示 ： Ai(r) Ai-1(r)+Ai(r) Ai-1(r)+R(p(r0， t)-iPAi-1(r) 其中， Ai(r) 表示本次迭代得到的。
8、更新后的重建光声图像， Ai-1(r) 表示上次迭代得到的 重建光声图像， (0， 1) 为残差图像的校正系数， Ai(r) 表示残差图像， R 表示离散化 反投影矩阵， p(r0， t) 表示信号采集位置 r0处 t 时刻的光声信号， i表示投影加权系数， P 表示离散化投影矩阵。 8.如权利要求7所述的光声成像重建方法， 其特征在于， 反投影矩阵R的大小为N2M， 其如下获得 ： 其中， N2为重建光声图像上所有像素点个数， M 为光声信号向量元素的个数， 大小为 采样位置个数和每个位置信号长度 k 的乘积 ； Rnm表示反投影矩阵 R 中第 n 行第 m 个元素 ； (xn， yn) 是。
9、重建光声图像的像素坐标， (x0， y0) 是信号采样位置坐 标， c 是声速， dt 是信号采集的时间间隔， m 是个数信号索引指标， n 是图像像素索引指标， mod 表示取余数运算， round 是四舍五入运算。 9.如权利要求1所述的光声成像重建方法， 其特征在于， 对步骤7进行二次优化得到本 轮迭代的重建光声图像如下表示 ： Ai(r) Ai-1(r)+R(p(r0， t)-PAi-1(r)+KAi-1(r) 其中， Ai(r) 表示本次迭代得到的更新后的重建光声图像， Ai-1(r) 表示上次迭代得到的 重建光声图像， (0， 1) 为残差图像的校正系数， R 表示离散化反投影矩阵。
10、， p(r0， t) 表 示信号采集位置 r0处 t 时刻的光声信号， i表示投影加权系数， P 表示离散化投影矩阵， 表示正则化参数， K 为正则化算子， 如下表示 ： 其中， 重建光声图像上一定大小邻域内 dmn为像素 A(m)和像素 A(n)之间距离的指数， 二 维空间中该邻域大小为 33， em为 RN的标准正交基， N 表示标准正交基空间的维数。 10. 一种基于正则化迭代的有限角度扫描光声成像的重建装置， 其包括 ： 初始重建光声图像获取模块， 其利用有限角度滤波反投影法得到初始的重建光声图像 A0(r) ； 投影加权系数计算模块， 其利用上轮迭代获得的重建光声图像计算投影加权系数。
11、， 对 权 利 要 求 书 CN 103310472 A 3 3/3 页 4 于第一轮迭代， 上轮迭代获得的重建光声图像为初始的重建光声图像 A0(r) ； 计算机模拟投影信号获取模块， 其利用上轮迭代获得的重建光声图像和投影加权系数 获取计算机模拟投影信号 ； 重建光声图像更新模块， 其计算所述计算机模拟投影信号与采集信号之间的信号残 差， 并根据所述信号残差获取残差图像并修正重建光声图像， 得到更新后的重建光声图 像 ； 优化模块， 其利用正则化运算对更新后的重建光声图像进行二次优化， 得到本轮迭代 的重建光声图像， 并由投影加权系数计算模块进行下轮迭代， 直至迭代完成。 权 利 要 求 。
12、书 CN 103310472 A 4 1/5 页 5 基于正则化迭代的有限角度光声成像重建方法及装置 技术领域 0001 本发明涉及一种光声成像 (Photoacoustic Tomography， 简称 PAT) 技术领域， 具 体涉及一种基于正则化迭代的有限角度扫描光声成像重建方法及装置。 背景技术 0002 光声成像是新发展起来的一种高分辨率和高对比度的生物医学成像技术， 近年来 获得迅速发展并已达到预临床阶段。光声图像传递的信息量大， 可以提供组织结构形态及 功能信息， 光声成像技术在生物组织成像中有广泛应用， 如血管成像、 肿瘤检测、 脑功能成 像等相关研究。 0003 光声信号的产。
13、生是复合媒介间能量转换的过程， 短脉冲激光照射到成像物体， 组 织局部吸收光能而产生热膨胀向四周辐射超声波， 通过超声换能器在不同位置采集的光声 信号， 使用重建算法算出生物组织中的吸收强度分布。成像方法是光声图像重建的关键技 术， 现有的重建算法有滤波反投影方法， 反卷积方法， 时域重建方法， 延迟求和方法等。 0004 上述图像重建算法都需要采集全方位的完备数据， 数据采集速度慢。本发明提出 的图像重建算法， 能基于有限角度扫描的欠采样信号精确地重建出光声图像。 发明内容 0005 本发明目的在于克服现有图像重建技术的缺点， 基于有限角度扫描的欠采样数据 和正则化迭代的光声成像重建方法， 。
14、提供一种生物组织吸收分布的快速精确成像方法， 解 决了生物组织无法全方位扫描时的问题。 0006 本发明提出的一种基于正则化迭代的有限角度扫描光声成像的重建方法， 其包 括 ： 0007 步骤 1、 利用有限角度滤波反投影法得到初始的重建光声图像 A0(r) ； 0008 步骤 2、 利用上轮迭代获得的重建光声图像计算投影加权系数， 对于第一轮迭代， 上轮迭代获得的重建光声图像为初始的重建光声图像 A0(r) ； 0009 步骤 3、 利用上轮迭代获得的重建光声图像和投影加权系数获取计算机模拟投影 信号 ； 0010 步骤 4、 计算所述计算机模拟投影信号与采集信号之间的信号残差， 并根据所述。
15、信 号残差获取残差图像并修正重建光声图像， 得到更新后的重建光声图像 ； 0011 步骤 5、 利用正则化运算对更新后的重建光声图像进行二次优化， 得到本轮迭代的 重建光声图像， 并转步骤 2 进行下轮迭代， 直至迭代完成。 0012 本发明提出的一种基于正则化迭代的有限角度扫描光声成像的重建装置， 其包 括 ： 0013 初始重建光声图像获取模块， 其利用有限角度滤波反投影法得到初始的重建光声 图像 A0(r) ； 0014 投影加权系数计算模块， 其利用上轮迭代获得的重建光声图像计算投影加权系 说 明 书 CN 103310472 A 5 2/5 页 6 数， 对于第一轮迭代， 上轮迭代获。
16、得的重建光声图像为初始的重建光声图像 A0(r) ； 0015 计算机模拟投影信号获取模块， 其利用上轮迭代获得的重建光声图像和投影加权 系数获取计算机模拟投影信号 ； 0016 重建光声图像更模块， 其计算所述计算机模拟投影信号与采集信号之间的信号 残差， 并根据所述信号残差获取残差图像并修正重建光声图像， 得到更新后的重建光声图 像 ； 0017 优化模块， 其利用正则化运算对更新后的重建光声图像进行二次优化， 得到本轮 迭代的重建光声图像， 并由投影加权系数计算模块进行下轮迭代， 直至迭代完成。 0018 本发明提出的上述方案， 基于有限角度的光声信号， 在每次迭代过程实现残差更 新图像。
17、和正则化更新图像， 计算采集信号和重建图像的计算机模拟信号的残差， 通过将残 差信号反投影得到残差图像， 将残差图像叠加到上一次重建图像得到更新图像， 然后利用 局部正则化更新重建图像， 结合信号残差反投影法和正则化方法获得重建图像， 本发明能 在有限角度扫描情形下， 快速精确地重建光声图像， 对减少重建时间、 降低设备硬件成本有 一定意义。 附图说明 0019 图 1 是本发明中基于有限角度数据光声成像重建的方法流程图。 0020 图 2 是本发明中计算机仿真成像光吸收分布和有限角度重建的结果示意图。 0021 图 3 是本发明中琼脂仿体的光吸收分布和有限角度重建的结果示意图。 具体实施方式。
18、 0022 为使本发明的目的、 技术方案和优点更加清楚明白， 以下结合具体实施例， 并参照 附图， 对本发明作进一步的详细说明。 0023 图 1 示出了本发明提出的基于正则化迭代的有限角度扫描光声成像的重建方法 流程图。如图 1 所示， 该方法的具体步骤如下 ： 0024 步骤 1、 利用有限角度滤波反投影成像方法得到迭代重建初始图像。 0025 用短脉冲激光均匀照射生物组织， 超声探头或者阵列在成像平面内进行扫描接收 光声信号， 根据热传导方程和超声的波动方程， 光声信号的产生服从的光声成像方程为 ： 0026 0027 其中， 代表汉密尔顿算子， p(r0， t) 表示信号采集位置 r0。
19、处 t 时刻的光声信号， c 表示生物组织中声波传播的速度， 表示声压膨胀系数， Cp为比热系数， A(r) 是组织内部 信号产生位置 r 处的光能量吸收分布的图像， I(t) 是激光脉冲能量函数。光声图像重建是 典型的逆问题， 即如何由 p(r0， t) 求出 A(r)。 0028 方程的 Green 函数解可以表示为 0029 说 明 书 CN 103310472 A 6 3/5 页 7 0030 (2) 式是表示光声信号和生物组织光吸收分布的关系， 将图像和采集信号离散化， 重建光吸收分布图像 A(r)， 大小为 NN， 将图像 A(r) 重排为 N2的向量， 采集信号 p(r0， t)。
20、， 长度为 M， 简记为 p， (2) 式转换为矩阵形式 ： 0031 p PA 0032 其中 p， P， A 的大小分别为 M1， MN2， N21。这里 P 表示投影矩阵， 对于每一个 采样信号 pm， 都先计算投影向量 Pm， 大小为 1N2， 对于任意的 1 n N2有 0033 0034 其中， Pmn表示投影矩阵 P 中第 m 行的第 n 个元素， 投影矩阵 P 的大小为 MN2， M 为光声信号向量元素的个数， 大小为采样位置个数和每个位置信号长度 k 的乘积， N 为重建 光声图像行和列的像素个数 ；(xn， yn)是光声图像像素坐标， (x0， y0) 是信号采集位置坐标，。
21、 c 是声速， dt 是信号采集的时间间隔， m 是光声信号索引指标， n 是图像像素索引指标， Smn表示第 n 个像素与第 m 条投影弧线相交面积的近似值， mod 是 指取余数运算。 0035 对于二维有限角度扫描， 通过引入有效扫描角的概念可以改进传统的滤波反投影 算法， 有限角度滤波反投影方法可近似表示为 ： 0036 0037 其中， 1为所述超声探头或者阵列探测位置的最小接收角度， 2为探测位置的 最大接收角度， e为有效扫描角度， 对于图像中的每个像素， 其定义为像素与最大接受位 置点之间线段和像素与最小接受位置点之间线段的夹角。式 (3) 表示采集信号的反投影过 程， 将其写。
22、为矩阵形式 ： 0038 A(r) Rp 0039 这里 R 表示反投影矩阵， 大小为 N2M， 对于每一个像素 A(n)， 计算光声信号 p 对其 的贡献， 对于任意的 1 m M 有 ： 0040 0041 其中， N2为重建光声图像上所有像素点个数， M 为光声信号向量元素的个数， 大小 为采样位置个数和每个位置信号长度 k 的乘积 ； Rnm表示反投影矩阵 R 中第 n 行的第 m 个元 素 ；(xn， yn) 是重建光声图像的像素坐标， (x0， y0) 是信号采样位 置坐标， c 是声速， dt 是信号采集的时间间隔， m 是个数信号索引指标， n 是图像像素索引指 说 明 书 C。
23、N 103310472 A 7 4/5 页 8 标， mod 表示取余数运算， round 是四舍五入运算。基于采集的光声信号， 利用式 (3) 计算得 到初始重建图像 A0(r)。 0042 步骤 2、 计算投影加权系数 i。 0043 根据上次迭代计算得到的迭代重建图像和采集的光声信号， 投影信号的加权系数 可计算如下 ： 0044 0045 其中， i 表示迭代次数， 对于第一轮迭代， 使用步骤 1 中计算得到的初始重建图像 作为上次迭代重建图像。 0046 步骤 3、 基于迭代图像获取计算机模拟投影信号。 0047 根据上轮迭代得到的迭代重建图像和步骤 2 计算所得的投影信号加权系数，。
24、 利用 离散化投影矩阵 P 获取计算机模拟投影信号 0048 iPAi-1(r) (5) 0049 步骤 4、 计算模拟信号与采集信号的残差。 0050 用实际采集信号减去步骤 3 获取的计算机模拟投影信号可得信号残差 ： 0051 p(r0， t)-iPAi-1(r) (6) 0052 步骤 5、 根据信号残差进而获取残差图像并修正重建图像。 0053 将步骤 4 得到的残差信号进行反投影 ： 0054 Ai(r) R(p(r0， t)-iPAi-1(r) (7) 0055 将残差信号反投影得到残差图像按一定比例叠加到初始重建图像， 可以得到更新 图像 ： 0056 Ai(r) Ai-1(r。
25、)+Ai(r) Ai-1(r)+R(p(r0， t)-iPAi-1(r) (8) 0057 其中， (0， 1) 为残差图像的校正系数。 0058 步骤 6、 对迭代更新图像实施正则化运算， 完成重建图像的二次优化。 0059 为了进一步对重建结果优化， 对重建过程加入正则化项。下面给出一个正则化改 进的重建算法的目标函数 ： 0060 0061 其中， 为正则化参数， dmn为一定大小的邻域内像素 A(m)和像素 A(n)之间距离的 指数， 该邻域大小优选为 33， 其余 dmn 0， N 表示空间的维数， 即标准正交基的维数 ； 利 用二次正则化的最小二乘理论， 该目标函数对像素求导并使得。
26、结果等于 0， 可以得到如下方 程 ： 0062 0063 其中， em为 RN的标准正交基， (10) 式右侧最后一项显然是重建图像的正则化项， 将正则化算子简化记为 K， 应用最速下降法， 可以得到目标函数的最优解为 ： 0064 Ai(r) Ai-1(r)+PT(p(r0， t)-PAi-1(r)+4KAi-1(r) (11) 说 明 书 CN 103310472 A 8 5/5 页 9 0065 将 (11) 式右侧的正则化项代入式 (8) 中， 迭代重建图像的过程更新为 ： 0066 Ai(r) Ai-1(r)+R(p(r0， t)-PAi-1(r)+KAi-1(r) (12) 00。
27、67 其中， K 为正则化算子， 如下表示 ： 0068 0069 其中， 重建光声图像上一定大小邻域内dmn为像素A(m)和像素A(n)之间距离的指数， 二维空间中该邻域大小为 33， em为 RN的标准正交基， N 表示标准正交基空间的维数。 0070 步骤 7、 若达到了迭代误差或最大迭代次数， 则输出重建图像， 否则转步骤 2， 重新 进行下一次迭代 ； 其中， 迭代次数或者迭代误差为预先设置的终止迭代条件。 通过上述迭代 可以得到重建光声图像。 0071 本发明中， 校正系数乘以投影算子计算校正后的模拟信号， 通过计算模拟信号和 采集信号的残差， 将残差信号反投影得到残差图像， 将残。
28、差图像按一定比例叠加到初始重 建图像， 可以得到更新图像。 0072 在计算机上进行仿真实验时， 根据光声信号产生方程， 设定已知的光吸收分布图， 并在 180和 90弧形扫描情况下采集光声信号， 扫描半径为 45mm， 扫描步长为 3， 用仿 真光声信号重建光吸收分布情况。 0073 在进行琼脂仿体的实验时， 用 1的脂肪乳剂， 6琼脂粉和 93的水混合加热到 70摄氏度， 然后注入直径30mm的圆柱内冷却制成琼脂仿体， 然后在其中加入直径0.5mm， 长 度为 10mm 和 5mm 的碳棒， 用来模拟生物组织的吸收分布。 0074 用波长为 532nm， 脉冲宽度为 6.5ns， 重复频率。
29、为 10Hz 倍频的 Q-Switched Nd:YAG 脉冲激光器作为激发光源， 脉冲激光经扩束、 均匀后照射到琼脂仿体。采用灵敏度 950mv/ Pa， 频率带宽 200KHz-15MHz， 直径为 1mm 的英国 Precision Acoustics 公司生产的 HP1 型 超声探头作为信号采集装置， 用 350MHz 的带宽， 取样速率最高可达 2.5GS/s 的泰克公司的 MSO4034数字示波器作为信号后处理装置。 探头的旋转实现由重复定位精度为0.005度， 分 辨率为 0.00125 度的 ERSP100 旋转台带动实现。 0075 图 2(a) 显示了计算机仿真的吸收分布图。
30、， 图像尺寸大小为 21mm21mm， 像素个数 为 256256。图 2(b) 显示了基于 180范围内 60 个采样信号经过 5 次迭代重建的光声图 像。 0076 图3(a)显示了琼脂仿体的照片， 图3(b)显示了基于90范围内30个采样信号经 过5次迭代重建的光声图像， 图3(c)显示了基于180范围内60个采样信号经过5次迭代 重建的光声图像。 0077 由琼脂仿体实验和计算机仿真实验的结果可以看出， 本发明的重建方法基于 180信号重建图像和原始图像基本一致， 说明本发明利用有限角度信号以少量迭代次数 精确地重建了光声图像， 对简化成像系统构造具有一定的实际意义。 0078 以上所述的具体实施例， 对本发明的目的、 技术方案和有益效果进行了进一步详 细说明， 应理解的是， 以上所述仅为本发明的具体实施例而已， 并不用于限制本发明， 凡在 本发明的精神和原则之内， 所做的任何修改、 等同替换、 改进等， 均应包含在本发明的保护 范围之内。 说 明 书 CN 103310472 A 9 1/1 页 10 图 1 图 2 图 3 说 明 书 附 图 CN 103310472 A 10 。