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

Abstract

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

Description

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

技术领域

本发明涉及一种光声成像（Photoacoustic Tomography，简称PAT）技术，具体涉及迭代自适应加权的有限视角光声图像重建方法。

背景技术

光声成像是近十年来新兴的生物医学成像技术，它结合了光学成像高对比度和超声成像的高分辨率的特点，近年来获得飞速的发展。光声图像不仅能够提供组织的结构形态，还能反映生理代谢等功能信息，光声成像技术在生物组织成像中已经得到广泛的应用，如肿瘤检测、血管成像、脑功能成像等相关研究。

短脉冲激光照射到成像样品，样品组织局部吸收光能而产生热膨胀向四周辐射超声波，通过超声探测器在不同位置的采集光声信号，使用图像重建方法算出组织的吸收强度分布。图像重建是由信号转化成图像的关键，而成像方法又是图像重建的核心技术，目前针对规则扫描形状已经提出了重建方法，包括滤波反投影法，时域重建法，频域重建法，延迟求和法等。

上述图像重建方法都在实际应用中，都需要全方位扫描获取完备数据，数据采集速度慢，对系统要求高。本发明提出的图像重建方法，对有限视角扫描的欠采样信息也能精确地重建出光声图像。

发明内容

（一）发明目的

本发明目的在于克服现有图像重建技术数据冗余的缺点，提供一种针对有限视角信息的光声成像方法，基于有限角度少量扫描的欠采样数据和迭代自适应加权方法，精确重建生物组织吸收强度分布图像，解决了生物组织无法全方位扫描的成像问题。

（二）技术方案

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

根据本发明的实施例，提供了一种基于迭代自适应加权的有限视角光声成像重建方法，包括以下步骤：

步骤1，重建图像初始化；

步骤2，计算加权系数；

步骤3，计算模拟信号；

步骤4，计算信号残差；

步骤5，计算重建图像；

步骤6，计算模型误差，设置迭代终止条件。

根据本发明的一个实施例，重建图像初始化包括采集光声信号的预处理，重建图像的离散化，以及计算信号与图像像素之间的关系。

根据本发明的一个实施例，计算信号与图像像素之间的关系包括投影矩阵和反投影矩阵的计算。

附图说明

图1为本发明的基于迭代自适应加权的有限视角光声成像重建方法流程图；

图2为计算机模拟投影信号的几何分布；

图3a、3b和图3c为基于本发明方法的计算机仿真数据重建结果，其中；其中图3a是用90°范围内30个探测信号重建的光声图像；图3b是用135°范围内45个探测信号重建的光声图像；图3c是用180°范围内60个探测信号重建的光声图像；

图4a和图4b为基于本发明方法的离体生物组织的重建结果，其中图4a为离体组织仿体照片，图4b本发明方法重建图像。

具体实施方式

用短脉冲激光均匀照射成像样品，超声探测器在成像平面内扫描，接收样品产生的光声信号，光声信号的产生服从的光声成像方程为：

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

其中

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

基于迭代自适应加权的有限视角光声成像重建方法如下：

第一步，重建图像初始化，包括采集光声信号的预处理和重建图像的离散化。采集光声信号的预处理是指采样数据的选取和滤波，根据成像区域的大小和位置在采集信号p_in中选取对应的信号，对选取信号进行Ramp滤波得到光声信号p(r₀,t)，完成对采集信号中高频噪声的预处理。

重建图像离散化即在重建区域上叠加一个N×N方格网，每一个网格内像素值相同，用A⁽ⁿ⁾表示第n像素的像素值，像素总数为N²；用p^(j)表示第j个计算机模拟光声信号，计算机模拟光声信号的总数与预处理后的采集光声信号相同，表示为M=m×k，它是投影位置个数m和每个位置数据长度k的乘积。如图2所示，A⁽ⁿ⁾和p^(j)的关系可以表示如下

$\underset{n=1}{\overset{N^2}{\mathrm{\Σ}}}{w}_{\mathrm{jn}}{A}^{\left(n\right)}=p^{\left(j\right)}j=\mathrm{1,2},\·\·\·,M---\left(2\right)$

式中w_jn表示第n像素值A⁽ⁿ⁾对第j个计算机模拟光声信号p^(j)的贡献，这里w_jn可以按以下方式计算：把投影弧线看成宽度为τ，间隔为0的弧线束，w_jn表示第n像素与第j条弧线相交的面积的近似值，将w_jn按信号采集顺序排列组成离散投影矩阵P。

对于二维的光声成像，通过引入有效扫描角的概念，改进传统的滤波反投影算法，得到有限视角的光声成像图像：

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

其中θ₁为探测器的最小接收角度，θ₂为探测器的最大接收角度，θ_e为有效扫描角度，其定义是每个像素点与信号采集位置最大接收角度和最小接收角度的夹角，为计算每个光声信号对图像像素的贡献，我们给出如下定义：

A(r)=Qp

其中A(r)表示重建图像，大小为N²×1，p表示光声信号，大小为M×1，Q是表示信号对图像像素的贡献值的反投影矩阵，大小为N²×M，对于每一个像素值A⁽ⁿ⁾，计算模拟光声信号p对其的贡献，对于任意的1≤j≤M有：

其中，Q_nj表示反投影矩阵Q中第n行的第j个元素，j是信号个数索引指标，n是图像像素索引指标，

(x_n,y_n)是重建光声图像的像素坐标，(x₀,y₀)是信号采样位置坐标，c是声速，f是信号采集的频率，mod表示取余数运算。

根据有限视角采集的光声信号，给定信号采集范围最大接受角度θ₂和最小接受角度θ₁，计算每个像素值对应的有效视角，通过傅里叶变换求得光声信号的求导运算值，然后根据每个光声信号对图像每个像素的贡献值，利用式（3）对非零贡献值进行积分运算，可以得到迭代重建光声图像的初始值A₀(r)。

第二步，加权系数的计算，用重建光声图像初始值A₀(r)和处理信号p(r₀,t)计算校正系数：

${\mathrm{\α}}_1=\frac{\max \left(p(r_0,t)\right)-\min \left(p(r_0,t)\right)}{\max \left(A_0\left(r\right)\right)-\min \left(A_0\left(r\right)\right)}---\left(4\right)$

第三步，模拟信号的计算，基于离散化的投影矩阵P和初始光声图像A₀(r)，我们可以计算计算机模拟光声信号p₁，

p₁=PA₀(r)

第四步，信号残差的计算，为了校正计算机模拟信号p₁和采集光声信号p(r₀,t)的差异，利用用校正系数α₁和模拟光声信号p₁，得到和采集光声信号采集信号同一量级的残差信号

Δp₁=p(r₀,t)-α₁p₁    (5)

第五步，重建图像的计算，将残差信号Δp₁反投影可以得到残差校正图像：

ΔA₁(r)=QΔp₁    (6)将残差校正图像ΔA₁(r)按一定比例加到初始重建图像A₀(r)，可以得到更新图像，即第一步迭代重建图像A₁(r)

A₁(r)=A₀(r)+λΔA₁(r)    (7)其中λ∈(0,1)为残差图像的校正系数，常根据重建图像效果的经验值设定，目的是防止过度校正。

重复上述步骤二至四及步骤五的上述步骤，可以得到迭代重建图像

A_i(r)=A_i-1(r)+λΔA_i(r)=A_i-1(r)+λQ(p(r₀,t)-α_ip_i)i=1,2,…  (8)

第六步，模型误差的计算，因为不知道初始图像，我们选择相对的重建误差e_i作为判定迭代终止的标准：

e_i=||A_i(r)-A_i-1(r)||₂/||A_i-1(r)||₂<ε(i=1,2,…n)    (9)其中ε为设定好的最小迭代误差。

设置最小迭代误差和/或最大迭代次数作为迭代终止条件，具体为：先判断重建误差是否小于设定的最小迭代误差，如果达到设定的最小迭代误差，则终止迭代；若没有达到设定的最小迭代误差，再判断迭代次数是否达到最大迭代次数，若达到最大迭代次数，则迭代终止。满足迭代终止条件则输出重建结果，通过上述迭代可以重建光声图像。发明人为验证本发明方法的有效性进行了如下实验：

计算机仿真实验：给定已知的光吸收分布图，在180°弧形扫描情况下采集光声信号，步长为3°，根据有限视角光声方程计算仿真信号，用仿真光声信号重建光吸收分布情况。

离体组织实验：在准备好的脂肪上凿三个深1mm、直径1mm的小孔，小孔之间的距离约为2mm，然后将准备好的猪肝样品填到脂肪中，最后再覆上一层厚度约为3mm的脂肪，成为一体后进行光声信号采集。

实验中采用的光源是Nd：YAG脉冲激光器（LS-2134，LOTIS），激光的波长是532nm，脉宽是7ns，重复频率是10Hz，脉冲激光经凹透镜扩束、毛玻璃均匀后照射到样品上，单脉冲能量控制在20mJ/cm²以下。采用中心频率5MHz，直径为13mm的美国OLYMPUS公司生产的超声探头采集信号，超声探测器和样品都放在耦合液里，目的是更好地耦合超声信号，探测器在一个位置采集完信号后，探头的旋转实现由重复定位精度为0.005度，分辨率为0.00125度的ERSP100旋转台带动实现，旋转半径为45mm，扫描间隔是3°，探测器采集信号经过前置放大器预处理，包括超低噪声时间增益补偿、模数转换等，信号通过一个采样频率为100MHz的数字示波器(MSO4000B,Tektronix)进行采集。探测器在一个位置采集50个信号后旋转到下一个位置，一共需要旋转60次，信号传输到计算机后，用MATLAB(version7.6,Mathworks)对信号处理，然后用基于迭代自适应加权的有限视角光声成像重建方法进行图像重建。

图3a、3b和3c显示了基于本发明方法的计算机仿真数据的重建结果，像素个数为256×256。其中图3a显示了用90°范围内30个探测信号重建的光声图像，图3b显示了用135°范围内45个探测信号重建的光声图像，图3c显示了用180°范围内60个探测信号重建的光声图像。

图4（a）显示了离体生物组织的照片，图4（b）显示了用180°范围内60个探测信号重建的图像。

由上述2组实验结果表明，本发明的重建方法在180°信号下重建图像和原始图像基本一致，说明本发明能在有限视角扫描信号下，精确地重建出组织的光吸收分布，对减少信号采集时间具有一定的实际意义。

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

Claims (11)

1.一种基于迭代自适应加权的有限视角光声成像重建方法，其特征在于包括以下步骤：

步骤1，重建图像初始化；

步骤2，计算加权系数；

步骤3，计算模拟信号；

步骤4，计算信号残差；

步骤5，计算重建图像；

步骤6，计算模型误差，设置迭代终止条件。

2.根据权利要求1所述的方法，所述步骤1包括采集光声信号的预处理，重建图像的离散化，以及计算信号与图像像素之间的关系。

3.根据权利要求2所述的方法，其中采集光声信号的预处理包括对采样数据的选取和滤波，预处理后的采集光声信号p(r₀,t)的总数表示为M=m×k，m为投影位置个数，k为每个位置的数据长度，滤波完成对采集信号中高频噪声的预处理。

4.根据权利要求2所述的方法，其中重建图像的离散化具体为：

在重建区域上叠加一个N×N方格网，每一个网格内像素值相同，用A⁽ⁿ⁾表示第n像素的像素值。

5.根据权利要求2所述的方法，其中计算信号与图像像素之间的关系包括投影矩阵和反投影矩阵的计算：

投影矩阵用来表示图像像素值A⁽ⁿ⁾对计算机模拟光声信号p的贡献，用p^(j)表示第j个计算机模拟光声信号，计算关系如下：

$\underset{n=1}{\overset{N^2}{\mathrm{\&Sigma;}}}{w}_{\mathrm{jn}}{A}^{\left(n\right)}=p^{\left(j\right)}j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,M$

其中w_jn表示第n像素值A⁽ⁿ⁾对第j个计算机模拟光声信号p^(j)的贡献值，w_jn按以下方式计算：把投影弧线看成宽度为τ，间隔为0的弧线束，w_jn表示第n像素与第j条弧线相交的面积的近似值，将w_jn按信号采集顺序排列组成离散投影矩阵P；

反投影矩阵用来表示计算机模拟信号p对重建图像像素值A⁽ⁿ⁾的贡献值，计算关系如下，对于任意的1≤j≤M有：

其中，Q_nj表示反投影矩阵Q中第n行的第j个元素，j是信号个数索引指标，n是图像像素索引指标，

(x_n,y_n)是重建光声图像的像素坐标，(x₀,y₀)是信号采样位置坐标，c是声速，f是信号采集的频率，mod表示取余数运算。

6.根据权利要求5所述的方法，对于二维的有限视角光声成像，重建图像A(r)和光声信号p(r₀,t)的关系近似为：

$A\left(r\right)\&ap;-\frac{r_0^2{C}_p}{2{\mathrm{\&pi;\&beta;c}}^4}{\&Integral;}_{{\mathrm{\&theta;}}_1}^{{\mathrm{\&theta;}}_2}\frac{1}{{\mathrm{\&theta;}}_{e}t}\frac{\&PartialD;p(r_0,t)}{\&PartialD;t}\mathrm{d\&theta;}{|}_{t=|r_0-r|/c}$

其中θ₁为探测器的最小接收角度，θ₂为探测器的最大接收角度，θ_e为有效扫描角度，表示每个像素点与信号采集位置最大接收角度和最小接收角度的夹角；

根据采集信号p(r₀,t)对图像每个像素的贡献值，利用上式计算得到重建光声图像的初始值A₀(r)。

7.根据权利要求1所述的方法，所述步骤2为采用重建光声图像初始值A₀(r)和处理信号p(r₀,t)利用下式计算校正系数：

${\mathrm{\&alpha;}}_1=\frac{\max \left(p(r_0,t)\right)-\min \left(p(r_0,t)\right)}{\max \left(A_0\left(r\right)\right)-\min \left(A_0\left(r\right)\right)}.$

8.根据权利要求1所述的方法，所述步骤3为基于离散化的投影矩阵P和初始光声图像A₀(r)，计算计算机模拟光声信号p₁，

p₁=PA₀(r)。

9.根据权利要求1所述的方法，所述步骤4为利用校正系数α₁和模拟光声信号p₁，得到和采集光声信号p(r₀,t)同一量级的残差信号

Δp₁=p(r₀,t)-α₁p₁。

10.根据权利要求1所述的方法，所述步骤5为将残差信号Δp₁反投影可以得到残差校正图像：

ΔA₁(r)=QΔp₁

将残差校正图像ΔA₁(r)按一定比例λ∈(0,1)加到初始重建图像A₀(r)，得到更新图像，即第一步迭代重建图像A₁(r)

A₁(r)=A₀(r)+λΔA₁(r)，

其中λ为残差图像的校正系数；

重复步骤2-4及步骤5中的上述步骤，得到迭代重建图像

A_i(r)=A_i-1(r)+λΔA_i(r)=A_i-1(r)+λQ(p(r₀,t)-α_ip_i)i=1,2,…。

11.根据权利要求1所述的方法，所述步骤6具体为计算重建误差e_i

e_i=||A_i(r)-A_i-1(r)||₂/||A_i-1(r)||₂；

设置最小迭代误差和/或最大迭代次数作为迭代终止条件，具体为：先判断重建误差是否小于设定的最小迭代误差，如果达到设定的最小迭代误差，则终止迭代；若没有达到设定的最小迭代误差，再判断迭代次数是否达到最大迭代次数，若达到最大迭代次数，则迭代终止。

Images