CN103310472A - 基于正则化迭代的有限角度光声成像重建方法及装置 - Google Patents

Abstract

本发明公开了一种基于正则化迭代的有限角度扫描光声成像的重建方法及装置。基于有限角度的光声信号，本发明在每次迭代过程实现残差更新图像和正则化更新图像，计算采集信号和重建图像的计算机模拟信号的残差，通过将残差信号反投影得到残差图像，将残差图像叠加到上一次重建图像得到更新图像，然后利用局部正则化更新重建图像，结合信号残差反投影法和正则化方法获得重建图像，本发明能在有限角度扫描情形下，快速精确地重建光声图像，对减少重建时间、降低设备硬件成本有一定意义。

Description

基于正则化迭代的有限角度光声成像重建方法及装置

技术领域

本发明涉及一种光声成像(Photoacoustic Tomography，简称PAT)技术领域，具体涉及一种基于正则化迭代的有限角度扫描光声成像重建方法及装置。

背景技术

光声成像是新发展起来的一种高分辨率和高对比度的生物医学成像技术，近年来获得迅速发展并已达到预临床阶段。光声图像传递的信息量大，可以提供组织结构形态及功能信息，光声成像技术在生物组织成像中有广泛应用，如血管成像、肿瘤检测、脑功能成像等相关研究。

光声信号的产生是复合媒介间能量转换的过程，短脉冲激光照射到成像物体，组织局部吸收光能而产生热膨胀向四周辐射超声波，通过超声换能器在不同位置采集的光声信号，使用重建算法算出生物组织中的吸收强度分布。成像方法是光声图像重建的关键技术，现有的重建算法有滤波反投影方法，反卷积方法，时域重建方法，延迟求和方法等。

上述图像重建算法都需要采集全方位的完备数据，数据采集速度慢。本发明提出的图像重建算法，能基于有限角度扫描的欠采样信号精确地重建出光声图像。

发明内容

本发明目的在于克服现有图像重建技术的缺点，基于有限角度扫描的欠采样数据和正则化迭代的光声成像重建方法，提供一种生物组织吸收分布的快速精确成像方法，解决了生物组织无法全方位扫描时的问题。

本发明提出的一种基于正则化迭代的有限角度扫描光声成像的重建方法，其包括：

步骤1、利用有限角度滤波反投影法得到初始的重建光声图像A₀(r)；

步骤2、利用上轮迭代获得的重建光声图像计算投影加权系数，对于第一轮迭代，上轮迭代获得的重建光声图像为初始的重建光声图像A₀(r)；

步骤3、利用上轮迭代获得的重建光声图像和投影加权系数获取计算机模拟投影信号；

步骤4、计算所述计算机模拟投影信号与采集信号之间的信号残差，并根据所述信号残差获取残差图像并修正重建光声图像，得到更新后的重建光声图像；

步骤5、利用正则化运算对更新后的重建光声图像进行二次优化，得到本轮迭代的重建光声图像，并转步骤2进行下轮迭代，直至迭代完成。

本发明提出的一种基于正则化迭代的有限角度扫描光声成像的重建装置，其包括：

初始重建光声图像获取模块，其利用有限角度滤波反投影法得到初始的重建光声图像A₀(r)；

投影加权系数计算模块，其利用上轮迭代获得的重建光声图像计算投影加权系数，对于第一轮迭代，上轮迭代获得的重建光声图像为初始的重建光声图像A₀(r)；

计算机模拟投影信号获取模块，其利用上轮迭代获得的重建光声图像和投影加权系数获取计算机模拟投影信号；

重建光声图像更模块，其计算所述计算机模拟投影信号与采集信号之间的信号残差，并根据所述信号残差获取残差图像并修正重建光声图像，得到更新后的重建光声图像；

优化模块，其利用正则化运算对更新后的重建光声图像进行二次优化，得到本轮迭代的重建光声图像，并由投影加权系数计算模块进行下轮迭代，直至迭代完成。

本发明提出的上述方案，基于有限角度的光声信号，在每次迭代过程实现残差更新图像和正则化更新图像，计算采集信号和重建图像的计算机模拟信号的残差，通过将残差信号反投影得到残差图像，将残差图像叠加到上一次重建图像得到更新图像，然后利用局部正则化更新重建图像，结合信号残差反投影法和正则化方法获得重建图像，本发明能在有限角度扫描情形下，快速精确地重建光声图像，对减少重建时间、降低设备硬件成本有一定意义。

附图说明

图1是本发明中基于有限角度数据光声成像重建的方法流程图。

图2是本发明中计算机仿真成像光吸收分布和有限角度重建的结果示意图。

图3是本发明中琼脂仿体的光吸收分布和有限角度重建的结果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

图1示出了本发明提出的基于正则化迭代的有限角度扫描光声成像的重建方法流程图。如图1所示，该方法的具体步骤如下：

步骤1、利用有限角度滤波反投影成像方法得到迭代重建初始图像。

用短脉冲激光均匀照射生物组织，超声探头或者阵列在成像平面内进行扫描接收光声信号，根据热传导方程和超声的波动方程，光声信号的产生服从的光声成像方程为：

${\▿}^{2}p(r,t)-\frac{1}{c^2}\frac{{\∂}^{2}p(r_0,t)}{{\∂t}^2}=\frac{\mathrm{\β}}{C_p}A\left(r\right)\frac{\∂I\left(t\right)}{\∂t}---\left(1\right)$

其中，

代表汉密尔顿算子，p(r₀，t)表示信号采集位置r₀处t时刻的光声信号，c表示生物组织中声波传播的速度，β表示声压膨胀系数，C_p为比热系数，A(r)是组织内部信号产生位置r处的光能量吸收分布的图像，I(t)是激光脉冲能量函数。光声图像重建是典型的逆问题，即如何由p(r₀，t)求出A(r)。

方程的Green函数解可以表示为

(2)式是表示光声信号和生物组织光吸收分布的关系，将图像和采集信号离散化，重建光吸收分布图像A(r)，大小为N×N，将图像A(r)重排为N²的向量，采集信号p(r₀，t)，长度为M，简记为p，(2)式转换为矩阵形式：

p＝PA

其中p，P，A的大小分别为M×1，M×N²，N²×1。这里P表示投影矩阵，对于每一个采样信号p_m，都先计算投影向量P_m，大小为1×N²，对于任意的1≤n≤N²有

其中，P_mn表示投影矩阵P中第m行的第n个元素，投影矩阵P的大小为M×N²，M为光声信号向量元素的个数，大小为采样位置个数和每个位置信号长度k的乘积，N为重建光声图像行和列的像素个数；

(x_n，y_n)是光声图像像素坐标，(x₀，y₀)是信号采集位置坐标，c是声速，dt是信号采集的时间间隔，m是光声信号索引指标，n是图像像素索引指标，ΔS_mn表示第n个像素与第m条投影弧线相交面积的近似值，mod是指取余数运算。

对于二维有限角度扫描，通过引入有效扫描角的概念可以改进传统的滤波反投影算法，有限角度滤波反投影方法可近似表示为：

$A\left(r\right)\≈-\frac{r_0^2{C}_p}{2\mathrm{\π\β}{c}^4}{\∫}_{{\mathrm{\θ}}_1}^{{\mathrm{\θ}}_2}\mathrm{d\θ}\frac{1}{{\text{\θ}}_{e}t}\frac{\∂p(r_0,t)}{\∂t}{|}_{t=|r_0-r|/c}---\left(3\right)$

其中，θ₁为所述超声探头或者阵列探测位置的最小接收角度，θ₂为探测位置的最大接收角度，θ_e为有效扫描角度，对于图像中的每个像素，其定义为像素与最大接受位置点之间线段和像素与最小接受位置点之间线段的夹角。式(3)表示采集信号的反投影过程，将其写为矩阵形式：

A(r)＝Rp

这里R表示反投影矩阵，大小为N²×M，对于每一个像素A⁽ⁿ⁾，计算光声信号p对其的贡献，对于任意的1≤m≤M有：

其中，N²为重建光声图像上所有像素点个数，M为光声信号向量元素的个数，大小为采样位置个数和每个位置信号长度k的乘积；R_nm表示反投影矩阵R中第n行的第m个元素；(x_n，y_n)是重建光声图像的像素坐标，(x₀，y₀)是信号采样位置坐标，c是声速，dt是信号采集的时间间隔，m是个数信号索引指标，n是图像像素索引指标，mod表示取余数运算，round是四舍五入运算。基于采集的光声信号，利用式(3)计算得到初始重建图像A₀(r)。

步骤2、计算投影加权系数λ_i。

根据上次迭代计算得到的迭代重建图像和采集的光声信号，投影信号的加权系数可计算如下：

${\mathrm{\λ}}_i=\frac{\max \left(p(r_0,t)\right)-\min \left(p(r_0,t)\right)}{\max \left(A_{i-1}\left(r\right)\right)-\min \left(A_{i-1}\left(r\right)\right)},i=\mathrm{1,2}\·\·\·---\left(4\right)$

其中，i表示迭代次数，对于第一轮迭代，使用步骤1中计算得到的初始重建图像作为上次迭代重建图像。

步骤3、基于迭代图像获取计算机模拟投影信号。

根据上轮迭代得到的迭代重建图像和步骤2计算所得的投影信号加权系数，利用离散化投影矩阵P获取计算机模拟投影信号

λ_iPA_i-1(r)    (5)

步骤4、计算模拟信号与采集信号的残差。

用实际采集信号减去步骤3获取的计算机模拟投影信号可得信号残差：

p(r₀，t)-λ_iPA_i-1(r)    (6)

步骤5、根据信号残差进而获取残差图像并修正重建图像。

将步骤4得到的残差信号进行反投影：

ΔA_i(r)＝R(p(r₀，t)-λ_iPA_i-1(r))    (7)

将残差信号反投影得到残差图像按一定比例叠加到初始重建图像，可以得到更新图像：

A_i(r)＝A_i-1(r)+αΔA_i(r)＝A_i-1(r)+αR(p(r₀，t)-λ_iPA_i-1(r))    (8)

其中，α∈(0，1)为残差图像的校正系数。

步骤6、对迭代更新图像实施正则化运算，完成重建图像的二次优化。

为了进一步对重建结果优化，对重建过程加入正则化项。下面给出一个正则化改进的重建算法的目标函数：

$F(p(r_0,t),A\left(r\right))=\frac{1}{2}{\left|\right|p(r_0,t)-\mathrm{PA}\left(r\right)\left|\right|}^2+\mathrm{\μ}\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{d}_{\mathrm{mn}}{(A^{\left(m\right)}\left(r\right)-A^{\left(n\right)}\left(r\right))}^2---\left(9\right)$

其中，μ为正则化参数，d_mn为一定大小的邻域内像素A^(m)和像素A⁽ⁿ⁾之间距离的指数，该邻域大小优选为3×3，其余d_mn＝0，N表示空间的维数，即标准正交基的维数；利用二次正则化的最小二乘理论，该目标函数对像素求导并使得结果等于0，可以得到如下方程：

$\frac{\∂F(p(r_0,t),A\left(r\right))}{{\∂A}^{\left(m\right)}\left(r\right)}=P^T(p(r_0,t)-\mathrm{PA}\left(r\right))+4\mathrm{\μ}\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\left(\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{d}_{\mathrm{mn}}(A^{\left(m\right)}\left(r\right)-A^{\left(n\right)}\left(r\right))\right)e_m=0,(m=\mathrm{1,2},\·\·\·N)---\left(10\right)$

其中，e_m为R^N的标准正交基，(10)式右侧最后一项显然是重建图像的正则化项，将正则化算子简化记为K，应用最速下降法，可以得到目标函数的最优解为：

A_i(r)＝A_i-1(r)+P^T(p(r₀，t)-PA_i-1(r))+4μKA_i-1(r)    (11)

将(11)式右侧的正则化项代入式(8)中，迭代重建图像的过程更新为：

A_i(r)＝A_i-1(r)+α[R(p(r₀，t)-PA_i-1(r))+μKA_i-1(r)]    (12)

其中，K为正则化算子，如下表示：

$K=\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\left(\underset{n-1}{\overset{N}{\mathrm{\Σ}}}{d}_{\mathrm{mn}}(A^{\left(m\right)}\left(r\right)-A^{\left(n\right)}\left(r\right))\right)e_m$

其中，重建光声图像上一定大小邻域内d_mn为像素A^(m)和像素A⁽ⁿ⁾之间距离的指数，二维空间中该邻域大小为3×3，e_m为R^N的标准正交基，N表示标准正交基空间的维数。

步骤7、若达到了迭代误差或最大迭代次数，则输出重建图像，否则转步骤2，重新进行下一次迭代；其中，迭代次数或者迭代误差为预先设置的终止迭代条件。通过上述迭代可以得到重建光声图像。

本发明中，校正系数乘以投影算子计算校正后的模拟信号，通过计算模拟信号和采集信号的残差，将残差信号反投影得到残差图像，将残差图像按一定比例叠加到初始重建图像，可以得到更新图像。

在计算机上进行仿真实验时，根据光声信号产生方程，设定已知的光吸收分布图，并在180°和90°弧形扫描情况下采集光声信号，扫描半径为45mm，扫描步长为3°，用仿真光声信号重建光吸收分布情况。

在进行琼脂仿体的实验时，用1％的脂肪乳剂，6％琼脂粉和93％的水混合加热到70摄氏度，然后注入直径30mm的圆柱内冷却制成琼脂仿体，然后在其中加入直径0.5mm，长度为10mm和5mm的碳棒，用来模拟生物组织的吸收分布。

用波长为532nm，脉冲宽度为6.5ns，重复频率为10Hz倍频的Q-Switched Nd:YAG脉冲激光器作为激发光源，脉冲激光经扩束、均匀后照射到琼脂仿体。采用灵敏度950mv/Pa，频率带宽200KHz-15MHz，直径为1mm的英国Precision Acoustics公司生产的HP1型超声探头作为信号采集装置，用350MHz的带宽，取样速率最高可达2.5GS/s的泰克公司的MSO4034数字示波器作为信号后处理装置。探头的旋转实现由重复定位精度为0.005度，分辨率为0.00125度的ERSP100旋转台带动实现。

图2(a)显示了计算机仿真的吸收分布图，图像尺寸大小为21mm×21mm，像素个数为256×256。图2(b)显示了基于180°范围内60个采样信号经过5次迭代重建的光声图像。

图3(a)显示了琼脂仿体的照片，图3(b)显示了基于90°范围内30个采样信号经过5次迭代重建的光声图像，图3(c)显示了基于180°范围内60个采样信号经过5次迭代重建的光声图像。

由琼脂仿体实验和计算机仿真实验的结果可以看出，本发明的重建方法基于180°信号重建图像和原始图像基本一致，说明本发明利用有限角度信号以少量迭代次数精确地重建了光声图像，对简化成像系统构造具有一定的实际意义。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于正则化迭代的有限角度扫描光声成像重建方法，其包括：

步骤1、利用有限角度滤波反投影法得到初始的重建光声图像A₀(r)；

步骤2、利用上轮迭代获得的重建光声图像计算投影加权系数，对于第一轮迭代，上轮迭代获得的重建光声图像为初始的重建光声图像A₀(r)；

步骤3、利用上轮迭代获得的重建光声图像和投影加权系数获取计算机模拟投影信号；

步骤4、计算所述计算机模拟投影信号与采集信号之间的信号残差，并根据所述信号残差获取残差图像并修正重建光声图像，得到更新后的重建光声图像；

步骤5、利用正则化运算对更新后的重建光声图像进行二次优化，得到本轮迭代的重建光声图像，并转步骤2进行下轮迭代，直至迭代完成。

2.如权利要求1所述的光声成像重建方法，其特征在于，所述初始的重建光声图像如下获得：

$A\left(r\right)\≈-\frac{r_0^2{C}_p}{2\mathrm{\π\β}{c}^4}{\∫}_{{\mathrm{\θ}}_1}^{{\mathrm{\θ}}_2}\mathrm{d\θ}\frac{1}{{\text{\θ}}_{e}t}\frac{\∂p(r_0,t)}{\∂t}{|}_{t=|r_0-r|/c}$

其中，A(r)表示生物组织内部的光吸收分布，p(r₀，t)表示信号采集位置r₀处t时刻的光声信号，c表示生物组织中声波传播的速度，β表示声压膨胀系数，C_p为比热系数，θ₁为探测位置的最小接收角度，θ₂为探测位置的最大接收角度，θ_e为有效扫描角度，对于图像中的每个像素，其定义为像素与最大接受位置点之间线段和像素与最小接受位置点之间线段的夹角。

3.如权利要求1所述的光声成像重建方法，其特征在于，所述投影加权系数如下计算：

${\mathrm{\λ}}_i=\frac{\max \left(p(r_0,t)\right)-\min \left(p(r_0,t)\right)}{\max \left(A_{i-1}\left(r\right)\right)-\min \left(A_{i-1}\left(r\right)\right)},i=\mathrm{1,2}\·\·\·$

其中，λ_i表示投影加权系数，A_i-1(r)表示上次迭代得到的生物组织内部信号产生位置r处的重建光声图像，p(r₀，t)表示信号采集位置r₀处t时刻的光声信号。

4.如权利要求1所述的光声成像重建方法，其特征在于，所述计算机模拟投影信号如下计算：

λ_iPA_i-1(r)

其中，λ_i表示投影加权系数，A_i-1(r)表示上次迭代得到的生物组织内部信号产生位置r处的重建光声图像，P表示离散化投影矩阵。

5.如权利要求4所述的光声成像重建方法，其特征在于，投影矩阵P如下计算：

其中，Pmn表示投影矩阵P中第m行的第n个元素，投影矩阵P的大小为M×N²，M为光声信号向量元素的个数，大小为采样位置个数和每个位置信号长度k的乘积，N²为重建光声图像上所有像素点个数；

(x_n，y_n)是光声图像像素坐标，(x₀，y₀)是信号采集位置坐标，c是声速，dt是信号采集的时间间隔，m是光声信号索引指标，n是图像像素索引指标，ΔS_mn表示第n个像素与第m条投影弧线相交面积的近似值，mod是指取余数运算。

6.权利要求1所述的光声成像重建方法，其特征在于，步骤4中所述信号残差通过实际采集到的信号p(r₀，t)减去所述计算机模拟投影信号λ_iPA_i-1(r)获得。

7.如权利要求1所述的光声成像重建方法，其特征在于，步骤4中，更新后的重建光声图像如下表示：

A_i(r)＝A_i-1(r)+αΔA_i(r)＝A_i-1(r)+αR(p(r₀，t)-λ_iPA_i-1(r))

其中，A_i(r)表示本次迭代得到的更新后的重建光声图像，A_i-1(r)表示上次迭代得到的重建光声图像，α∈(0，1)为残差图像的校正系数，ΔA_i(r)表示残差图像，R表示离散化反投影矩阵，p(r₀，t)表示信号采集位置r₀处t时刻的光声信号，λ_i表示投影加权系数，P表示离散化投影矩阵。

8.如权利要求7所述的光声成像重建方法，其特征在于，反投影矩阵R的大小为N²×M，其如下获得：

其中，N²为重建光声图像上所有像素点个数，M为光声信号向量元素的个数，大小为采样位置个数和每个位置信号长度k的乘积；R_nm表示反投影矩阵R中第n行第m个元素；

(x_n，y_n)是重建光声图像的像素坐标，(x₀，y₀)是信号采样位置坐标，c是声速，dt是信号采集的时间间隔，m是个数信号索引指标，n是图像像素索引指标，mod表示取余数运算，round是四舍五入运算。

9.如权利要求1所述的光声成像重建方法，其特征在于，对步骤7进行二次优化得到本轮迭代的重建光声图像如下表示：

A_i(r)＝A_i-1(r)+α[R(p(r₀，t)-PA_i-1(r))+μKA_i-1(r)]

其中，A_i(r)表示本次迭代得到的更新后的重建光声图像，A_i-1(r)表示上次迭代得到的重建光声图像，α∈(0，1)为残差图像的校正系数，R表示离散化反投影矩阵，p(r₀，t)表示信号采集位置r₀处t时刻的光声信号，λ_i表示投影加权系数，P表示离散化投影矩阵，μ表示正则化参数，K为正则化算子，如下表示：

$K=\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\left(\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{d}_{\mathrm{mn}}(A^{\left(m\right)}\left(r\right)-A^{\left(n\right)}\left(r\right))\right)e_m$

其中，重建光声图像上一定大小邻域内d_mn为像素A^(m)和像素A⁽ⁿ⁾之间距离的指数，二维空间中该邻域大小为3×3，e_m为R^N的标准正交基，N表示标准正交基空间的维数。

10.一种基于正则化迭代的有限角度扫描光声成像的重建装置，其包括：

初始重建光声图像获取模块，其利用有限角度滤波反投影法得到初始的重建光声图像A₀(r)；

投影加权系数计算模块，其利用上轮迭代获得的重建光声图像计算投影加权系数，对于第一轮迭代，上轮迭代获得的重建光声图像为初始的重建光声图像A₀(r)；

计算机模拟投影信号获取模块，其利用上轮迭代获得的重建光声图像和投影加权系数获取计算机模拟投影信号；

重建光声图像更新模块，其计算所述计算机模拟投影信号与采集信号之间的信号残差，并根据所述信号残差获取残差图像并修正重建光声图像，得到更新后的重建光声图像；

优化模块，其利用正则化运算对更新后的重建光声图像进行二次优化，得到本轮迭代的重建光声图像，并由投影加权系数计算模块进行下轮迭代，直至迭代完成。

Images