CN103584835A - 一种基于压缩感知的光声图像重建方法 - Google Patents

Abstract

本发明公开了一种基于压缩感知的光声图像重建方法，包括以下步骤：在目标组织周围放置M个传感器；用宽束短脉冲激光周期性照射组织W次，各传感器采集光声信号；对各传感器采集的W次光声信号求平均，得到M组平均光声信号；对每组平均光声信号做FFT，得到M组频谱数据；设定M个随机向量，根据各随机向量对各组频谱数据进行N点随机抽样；选取正交基矩阵，根据传感器脉冲响应及随机抽样得到的频率值计算测量矩阵；利用压缩感知求解目标生物组织的光吸收分布，得到光声重建图像。本发明结合压缩感知，利用少量传感器采集的光声信号完成目标生物组织内部的三维图像重建，具有传感器数量少、传感器摆放位置相对任意等优点。

Description

一种基于压缩感知的光声图像重建方法

技术领域

本发明涉及图像处理领域，特别是一种基于压缩感知的光声图像重建方法。

背景技术

光声成像以超声波作为信息的载体，反映生物组织的空间光学吸收分布特性，它同时具备了纯光学成像和超声成像二者的优点，在保证一定成像深度的同时又具有高的对比度和空间分辨率，尤其适合于对生物软组织的成像。光声成像可根据待测目标的光学吸收特性选取合适的激励光波长，可定量的测量与光学吸收相关的生理参数，并可在完全无损伤的情况下对生物的高级功能活动进行观察并提供各方向高分辨率和高对比度的组织断层图像。尽管光声成像具备以上优势，但由于使用单传感器进行圆周扫描采集信号周期长或多传感器阵列受制于传感器的带宽、阵元数及间距而制造价格昂贵，该技术若要步入到临床实践阶段，仍需做出大量研究和改善。

压缩感知技术是一种新的采样理论，其认为只要信号是可压缩的或在某个变换域中是稀疏的，就可以利用随机采样直接得到压缩数据，然后通过求解最优化问题从采集到的少量压缩数据中重建出原始信号。利用压缩感知技术可以大大减少光声重建需要采集的数据，从而可以利用少量观测角度采集到的信息或较少数量传感器采集到的信息进行光声重建，大大节省传感器采样时间和制造成本。

发明内容

发明目的：本发明所要解决的技术问题是针对传统光声重建方法信号采样周期长或传感器阵列制造成本高的问题，提供一种基于压缩感知的光声图像重建方法。

为了解决上述技术问题，本发明公开了一种基于压缩感知技术的光声图像重建方法，包括以下步骤：

步骤一，在目标组织周围放置M个传感器；

步骤二，用宽束短脉冲激光周期性照射组织W次，各传感器采集光声信号；

步骤三，对各传感器采集的W次光声信号求平均，得到M组平均光声信号；

步骤四，对每组平均光声信号做FFT，得到M组频谱数据；

步骤五，设定M个随机向量，根据各随机向量对各组频谱数据进行N点随机抽样；

步骤六，选取正交基矩阵，根据传感器脉冲响应及随机抽样得到的频率值计算测量矩阵；

步骤六，利用压缩感知求解目标生物组织的光吸收分布，得到光声重建图像。

本发明中，优选地，所述步骤一中M个传感器可以放置在目标生物组织周围的任何位置，比如放置于组织一侧的矩形网格中的任意网格点上，或放置在包围组织的整个球体网格中的任意网格点上，M的取值与所要重建光声图像经过正交基矩阵变换后的图像的稀疏度有关，图像越稀疏，所需要的传感器个数就越少，通常可选择在30～150之间。

本发明中，优选地，所述步骤二中激光发射器周期性的发射宽束短脉冲激光，照射整个生物组织。每发射一次激光，所有传感器就采集组织发出的光声信号，重复发射W次，接收过程中不需改变传感器位置。

本发明中，步骤三对每个传感器采样得到的W组光声信号求平均，得到M组平均光声信号，以减少随机噪声干扰，其中W的值太小重建图像受噪声影响较大，W的值太大将耗费大量信号采样时间，通常可选择在1～128之间。

本发明中，优选地，所述步骤四中需要对每组平均光声信号进行快速傅里叶变换(FFT)，得到对应的频谱数据。

本发明中，优选地，所述步骤五中针对不同传感器计算得到的频谱数据产生不同的随机向量对该频谱数据进行N点随机抽样，将M组抽样得到的N点频谱数据合并为一个M行N列的频谱抽样矩阵Y，其中N即为重建图像宽度。

本发明中，优选地，所述步骤六中正交基矩阵Φ应根据需要重建图像在该正交基变换后的系数是稀疏的或可压缩的为前提进行选择，通常可以是傅里叶正交基、小波正交基、曲波正交基中的一种；测量矩阵K可根据传感器和声源的相对位置、随机抽样得到的频率值及传感器脉冲响应计算得到。

本发明中，优选地，所述步骤七中利用压缩感知技术，根据随机抽样得到的频谱抽样矩阵Y、用于稀疏信号表示的正交基矩阵Φ以及计算得到的测量矩阵K，通过求解最优化问题，求解得到目标生物组织的空间光吸收分布。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1是本发明方法的流程图。

图2是本发明方法的设备安放示意图。

具体实施方式：

本发明结合压缩感知技术，通过随机放置的少数传感器采集光声信号，并转换到频域进行随机抽样，利用随机抽样得到的频谱数据、计算得到的测量矩阵以及用于稀疏化重建图像的正交基矩阵，通过求解一个最优化过程重建出生物组织的空间光吸收分布。

如图1所示，本发明公开了一种基于压缩感知的光声图像重建方法，包括以下步骤：

步骤一，在目标组织周围放置M个传感器；

步骤二，用宽束短脉冲激光周期性照射组织W次，各传感器采集光声信号；

步骤三，对各传感器采集的W次光声信号求平均，得到M组平均光声信号；

步骤四，对每组平均光声信号做FFT，得到M组频谱数据；

步骤五，设定M个随机向量，根据各随机向量对各组频谱数据进行N点随机抽样；

步骤六，选取正交基矩阵，根据传感器脉冲响应及随机抽样得到的频率值计算测量矩阵；

步骤六，利用压缩感知求解目标生物组织的光吸收分布，得到光声重建图像。

本发明中，步骤一，将M(30～150)个传感器安放在目标生物组织周围的任何位置，比如放置于组织一侧的矩形网格中的任意网格点上，也可以放置在包围组织的整个球体网格中的任意网格点上。一种安放激光发生器、生物组织、球体网格传感器阵列等装置的示意图如图2所示，其中透镜用来将较窄的激光束进行扩束，以覆盖到整个生物组织。

本发明中，步骤二，激光发射器周期性的发射宽束短脉冲激光，照射整个生物组织。每发射一次激光，所有传感器就采集组织发出的光声信号，重复发射W(1～128)次，接收过程中不需改变传感器位置。

本发明中，步骤三对每个传感器采样得到的W组光声信号求平均，得到M组平均光声信号

${\stackrel{\‾}{p}}_m({\stackrel{\→}{r}}_m,t)=\frac{1}{W}\underset{s=1}{\overset{W}{\mathrm{\Σ}}}{p}_{m,s}({\stackrel{\→}{r}}_m,t)---\left(1\right)$

其中，W为每个传感器采集次数，

为第m(m＝1，2，...，M)个传感器第s(s＝1，2，...，W)次采集到的光声信号，t表示时间，

为第m个传感器所在位置，一般来说，光声信号信噪比很低，不到10dB，平均后可以抑制背景噪声。

本发明中，步骤四，

与组织的光吸收分布有如下关系：

$p_{m,s}({\stackrel{\→}{r}}_m,t)=\frac{\mathrm{\β}}{4\mathrm{\πC}}\∫d{\stackrel{\→}{r}}^{\′}A\left({\stackrel{\→}{r}}^{\′}\right)\frac{{\mathrm{\δ}}^{\′}(t-|{\stackrel{\→}{r}}_m-{\stackrel{\→}{r}}^{\′}|/c)}{|{\stackrel{\→}{r}}_m-{\stackrel{\→}{r}}^{\′}|}\⊗g_m\left(t\right)---\left(2\right)$

其中，c为声速，C为比热，β为等压膨胀系数，g_m(t)为第m个传感器的时域脉冲响应函数，

为处组织的光吸收分布，

表示第m个传感器到

处组织的欧氏距离。

公式(2)在傅里叶域的表示形式为：

${\stackrel{\‾}{p}}_m({\stackrel{\→}{r}}_m,k)=-\mathrm{ick}{g}_m\left(k\right)\∫d{\stackrel{\→}{r}}^{\′}A\left({\stackrel{\→}{r}}^{\′}\right)\frac{e^{\mathrm{ik}|{\stackrel{\→}{r}}_m-{\stackrel{\→}{r}}^{\′}|}}{|{\stackrel{\→}{r}}_m-{\stackrel{\→}{r}}^{\′}|}---\left(3\right)$

其中k＝2πf/c，f为频率值，g_m(k)为第m个传感器的频域脉冲响应函数，i为虚数单位。

本发明中，步骤五，对不同传感器计算得到的频谱数据根据设定的N点随机向量进行N点随机抽样。处理完所有频谱数据后，可将M组抽样得到的N点频谱数据合并为一个M×N的频谱抽样矩阵Y：

$Y={\left[\begin{array}{cccc}{\stackrel{\‾}{p}}_1({\stackrel{\→}{r}}_1,k_{\mathrm{1,1}})& {\stackrel{\‾}{p}}_1({\stackrel{\→}{r}}_1,k_{\mathrm{1,2}})& ...& {\stackrel{\‾}{p}}_1({\stackrel{\→}{r}}_1,k_{1,N})\\ {\stackrel{\‾}{p}}_2({\stackrel{\→}{r}}_2,k_{\mathrm{2,1}})& {\stackrel{\‾}{p}}_2({\stackrel{\→}{r}}_2,k_{\mathrm{2,2}})& ...& {\stackrel{\‾}{p}}_2({\stackrel{\→}{r}}_2,k_{2,N})\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\stackrel{\‾}{p}}_M({\stackrel{\→}{r}}_M,k_{M,1})& {\stackrel{\‾}{p}}_M({\stackrel{\→}{r}}_M,k_{M,2})& ...& {\stackrel{\‾}{p}}_M({\stackrel{\→}{r}}_M,k_{M,N})\end{array}\right]}_{M\×N}---\left(4\right)$

其中k_m，n＝2πf_m，n/c，f_m，n为第m(m＝1，...，M)个传感器的第n(n＝1，2，...，N)个频谱抽样值。

本发明中，步骤六，为能够从步骤四中得到的少量采样数据中重建出生物组织的光吸收分布，认为重建得到的图像在某种正交基变换后的系数是稀疏的或可压缩的，从而可利用压缩感知技术从矩阵Y中恢复出组织的光吸收分布

正交基矩阵Φ可以是傅里叶正交基、小波正交基、曲波正交基等，另外，KΦ矩阵中各列向量之间的相关性越小，所得到的重建图像越好。此处以傅里叶正交基矩阵为例：

$\mathrm{\Φ}=\frac{1}{N}{\left[\begin{array}{cccc}{\mathrm{\ω}}^0& {\mathrm{\ω}}^0& ...& {\mathrm{\ω}}^0\\ {\mathrm{\ω}}^0& {\mathrm{\ω}}^1& ...& {\mathrm{\ω}}^{\mathrm{NN}-1}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\ω}}^0& {\mathrm{\ω}}^{\mathrm{NN}-1}& ...& {\mathrm{\ω}}^{(\mathrm{NN}-1)*(\mathrm{NN}-1)}\end{array}\right]}_{\mathrm{NN}\×\mathrm{NN}}---\left(5\right)$

其中ω是1的n次方根的主值，大小为e^-2πi/NN，最终重建图像的大小为N×N。

测量矩阵K可根据公式(3)的离散化结果得到：

$K_{m,n}(s,l)=-\mathrm{ic}{k}_{m,n}\frac{e^{i{k}_{m,n}|{\stackrel{\→}{r}}_m-{\stackrel{\→}{r}}_{s,l}|}}{|{\stackrel{\→}{r}}_m-{\stackrel{\→}{r}}_{s,l}|}{g}_{m,n}---\left(6\right)$

其中(s，l)为需要重建的组织中某断层图像上某像素点的坐标，g_m，n为第m个传感器的脉冲响应在第n个位置处的值。

本发明中，步骤七，根据公式(6)计算得到的测量矩阵K、步骤六中选取的正交基矩阵Φ以及采用公式(4)随机抽样得到的频谱抽样矩阵Y求解压缩感知最优化问题：

$\begin{array}{ccc}\min {\left|\right|\mathrm{\θ}\left|\right|}_{1_1}& s.t.& Y=\mathrm{K\Φ\θ}\end{array}---\left(7\right)$

A＝Φθ    (8)

其中A为目标生物的光吸收分布矩阵，即所需求解变量。θ为矩阵A经正交基Φ变换后的光吸收分布矩阵。

求解公式(7)的方法有很多，以正交匹配跟踪算法(OMP)为例求解公式(7)：

第一步：初始化残差Δ⁽⁰⁾＝Y，其中Y为公式(4)计算出的频谱抽样矩阵，系数集

迭代次数t＝1，最大迭代次数为T；

第二步：令等效矩阵D＝KΦ，其中K由公式(5)求得，Φ为步骤六选取。计算残差与等效矩阵D中各列的相关性，找到λ_t＝argmax_{j＝1，2，...，N}|＜Δ^(t-1)D_j＞|，其中D_j为D的第j列，D_j即称为矩阵D的第j个原子；

第三步：更新系数集：Λ_t＝Λ_t-1U{λ_t}，更新重建原子集合： ${\mathrm{\Ψ}}_t=\left[\begin{array}{cc}{\mathrm{\Ψ}}_{t-1}& D_{{\mathrm{\λ}}_t}\end{array}\right],$ 其中Ψ₀为空；

第四步：解最小二乘方程，得到需要估计的新的信号θ_t＝argmin_θ||Y-Ψ_tθ||₂；

第五步：更新残差：Δ^(t)＝Y-Ψ_tθ；

第六步：t＝t+1，重复第二步，直到t＝T。

最终迭代得到的θ_T即为所要求得的稀疏系数矩阵θ，然后根据公式(8)得到生物组织的光吸收分布A，所求得的光吸收分布反应了光声重建得到的组织内部的图像。

在具体求解时，可以将M×N大小的Y转为MN×1的列向量，将M×N×N×N四维矩阵K转为MN×NN二维矩阵。利用维数转换后的Y、K以及正交基Φ，结合OMP算法求解θ。最后利用公式(8)得到光吸收分布A。

本发明提供了一种基于压缩感知技术的光声图像重建方法，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。另外，本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (8)

1.一种基于压缩感知的光声图像重建方法，其特征在于，包括以下步骤： 

步骤一，在目标组织周围放置M个传感器； 

步骤二，用宽束短脉冲激光周期性照射组织W次，各传感器采集光声信号； 

步骤三，对各传感器采集的W次光声信号求平均，得到M组平均光声信号； 

步骤四，对每组平均光声信号做FFT，得到M组频谱数据； 

步骤五，设定M个随机向量，根据各随机向量对各组频谱数据进行N点随机抽样； 

步骤六，选取正交基矩阵，根据传感器脉冲响应及随机抽样得到的频率值计算测量矩阵； 

步骤七，利用压缩感知求解目标生物组织的光吸收分布，得到光声重建图像。 

2.根据权利要求1所述的一种基于压缩感知的光声图像重建方法，其特征在于，所述步骤一中放置传感器，可以任意放置在位于生物组织一侧的矩形网格中，或放置在包围生物组织的球体网格中。 

3.根据权利要求1所述的一种基于压缩感知的光声图像重建方法，其特征在于，所述步骤二激光发生器每发射一次激光，M个传感器同时采集一段光声信号，并重复W次，测量中，不需改变传感器位置。 

4.根据权利要求1所述的一种基于压缩感知的光声图像重建方法，其特征在于，所述步骤三分别对每个传感器采样得到的W组光声信号求平均，得到M组平均光声信号。 

5.根据权利要求1所述的一种基于压缩感知的光声图像重建方法，其特征在于，所述步骤四需要分别对每组平均光声信号进行快速傅里叶变换(FFT)，得到对应的频谱数据。 

6.根据权利要求1所述的一种基于压缩感知的光声图像重建方法，其特征在于，所述步骤五需要分别为各组频谱数据设定随机向量并对其进行N点随机抽样，将M组频谱抽样结果合并为一个M×N大小的频谱抽样矩阵Y。 

7.根据权利要求1所述的一种基于压缩感知的光声图像重建方法，其特征在于，所述步骤六中选择的正交基矩阵Φ可以是傅里叶正交基、小波正交基、曲波正交基中的一种，也可以选取其它正交基矩阵；所述步骤六中测量矩阵K可根据传感器和声源的位置、随机抽样得到的频率值及传感器脉冲响应计算得 到。 

8.根据权利要求1所述的一种基于压缩感知的光声图像重建方法，其特征在于，所述步骤七中利用压缩感知理论，使用压缩感知信号重建时的最优化求解方法，结合选择的正交基Φ、计算得到的测量矩阵K以及随机抽样得到的频谱抽样矩阵Y，求解得到生物组织的光吸收分布，即为光声图像的重建结果。 

Images