CN109671130B - 利用稀疏测量数据重建内窥式光声层析图像的方法及*** - Google Patents

Abstract

本发明公开了一种利用稀疏测量数据重建内窥式光声层析图像的方法及***，解决了利用稀疏的光声测量数据重建光吸收分布图像中出现严重伪影和失真的问题。首先构造可用于光声信号稀疏表示的自适应完备字典；其次利用随机测量矩阵测量到数据构建稀疏光声信号测量数据矩阵；然后根据测量矩阵和自适应完备字典计算感知矩阵；再者根据稀疏光声信号测量数据矩阵和感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵；最后利用系数矩阵和自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵，并采用图像重建公式，重建光吸收能量分布图，提高了成像精度。

Description

利用稀疏测量数据重建内窥式光声层析图像的方法及***

技术领域

本发明涉及医学成像技术领域，一种利用稀疏光声信号测量数据重建生物腔体横截面的光声层析图像的方法及***。

背景技术

光声层析成像(photoacoustic tomography，PAT)是一种基于生物组织光声效应的非电离功能成像方法，其成像参数是组织的光吸收系数和散射系数，可以实现高分辨率和高对比度的软组织深层成像。PAT的原理是采用短脉冲激光照射生物组织，组织吸收光能量后受热膨胀产生瞬时压力，并向外辐射宽带(10kHz～100MHz)超声波，即光声信号。声压的幅值与脉冲激光的强度成正比，反映组织的光吸收特性。超声换能器接收来自不同方向、不同位置的光声信号，送入计算机后采用合适的算法即可反演得到组织内部的初始声压或者光吸收能量的空间分布图，直观显示组织的内部结构。在此基础上还可估算组织的光学特性参数(光吸收系数和散射系数)的空间分布，反映组织的功能成分。

在实际应用中，特别是内窥式PAT(如血管内光声成像)，由于腔道内封闭成像几何的特殊性，受成像导管的机械结构、空间位置及成像时间等的限制，超声探测器往往只能进行有限角度的扫描，采集到稀疏的光声信号数据，进而导致重建出的光吸收分布图像中出现严重的伪影和失真，降低图像质量。因此，为了提高成像精度，需要解决利用稀疏的光声测量数据重建高质量PAT图像的问题。

发明内容

为了提高成像精度，本发明提供了一种利用稀疏测量数据重建内窥式光声层析图像的方法及***，解决了利用稀疏的光声测量数据重建光吸收分布图像中出现严重伪影和失真的问题。

为实现上述目的，本发明提供了如下方案：

一种利用稀疏测量数据重建内窥式光声层析图像的方法，所述方法包括：

构造可用于光声信号稀疏表示的自适应完备字典；

构建测量矩阵；

根据所述测量矩阵和所述自适应完备字典，计算感知矩阵；

利用所述测量矩阵获取稀疏光声信号测量数据，并构建稀疏光声信号测量数据矩阵；

根据所述稀疏光声信号测量数据矩阵和所述感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵；

利用所述系数矩阵和所述自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵；

根据所述完备光声信号数据矩阵，采用图像重建公式，重建光吸收能量分布图。

可选的，所述构造可用于光声信号稀疏表示的自适应完备字典，具体包括：

构建初始完备光声信号测量数据矩阵；所述初始完备光声信号测量数据矩阵中的元素为在N个测量位置处获取的完备光声信号测量数据，且在各测量位置处采集的光声信号的长度为L；所述初始完备光声信号测量数据矩阵为N×L维的矩阵P；

采用离散余弦变换基作为初始完备字典D₀；所述初始完备字典D₀是N×N维的初始完备字典，N是测量位置总数；

以所述初始完备光声信号测量数据矩阵为训练样本，按照公式

计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'；

其中，D是更新后的N×N维的完备字典；矩阵Θ₀是矩阵P基于初始完备字典D₀的初始稀疏表示的系数矩阵，大小是N×L维；向量d_j和d_k分别是更新后的完备字典D的第j列和第k列元素，即第j个原子和第k个原子；||·||₂表示2-范数；向量θ₀ ^j和θ₀ ^k分别是矩阵Θ₀的第j行和第k行元素；向量θ_0i是矩阵Θ₀的第i列元素，其中i＝1,2,...,L；

计算矩阵P基于更新后的完备字典D的稀疏表示的系数矩阵Θ；系数矩阵Θ的大小是N×L维；

判断矩阵P更新后的完备字典D以及系数矩阵Θ是否满足

若是则停止迭代，并将当前迭代对应的更新后的完备字典D确定为自适应完备字典；

若否则将更新后的初始稀疏表示的系数矩阵Θ₀'代替初始稀疏表示的系数矩阵矩阵Θ₀，返回计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'步骤。

可选的，所述计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'，具体包括：

以所述初始完备光声信号测量数据矩阵为训练样本，按照公式

采用奇异值分解算法，计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'。

可选的，所述构建测量矩阵，具体包括：

生成N×N维的哈达玛矩阵；

在所述哈达玛矩阵中随机地选取M行向量，构建测量矩阵；所述构建测量矩阵为一个M×N维的矩阵，其中M<<N，M是稀疏测量位置数。

可选的，所述根据所述稀疏光声信号测量数据矩阵和所述感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵，具体包括：

步骤1：根据稀疏光声信号测量数据矩阵Y和感知矩阵A，得到完备光声信号矩阵P_Y的稀疏表示的系数矩阵Θ_Y；

步骤2：初始化参数；r₀＝y_i，size＝4，stage＝1，m＝1，

其中各参数的含义为：向量y_i是稀疏光声信号测量数据矩阵Y的第i列元素；r₀是初次余量，r_m是第m次迭代的余量；size是搜索步长；stage是阶段；m是迭代次数；U是相关系数矩阵；J是相关系数矩阵U按照从大到小排列，前size个值对应感知矩阵A的列序号集合；J₀是每次迭代找到的列序号集合；Λ_m是第m次迭代从感知矩阵A中选出的列向量的列序号集合；A_m是按索引Λ_m选出的感知矩阵A的列向量集合；

是空集；

步骤3：计算当前迭代余量r_m；

步骤4：判断当前迭代余量r_m是否满足||r_m||₂≤ε₁；其中ε₁是控制迭代次数的阈值，若是，则停止迭代，并根据上次迭代确定的每个完备光声信号稀疏表示的系数构建系数矩阵

若否，则转向步骤5；

步骤5：计算相关系数矩阵U；所述相关系数矩阵U的计算公式为U＝{u_k|u_k＝<r_m-1,α_k>|,k＝1,2,3,…,N}；其中，向量u_k是相关系数矩阵U的第k列，向量α_k是感知矩阵A的第k列，<·>是向量的内积运算；

步骤6：将相关系数矩阵U中的元素按照从大到小进行排列，选择前size个元素对应感知矩阵A的列序号来构成集合J；

步骤7：正则化；在集合J中寻找满足约束条件的子集J₀，并在所有满足约束条件的子集J₀中选择具有最大能量的子集；所述约束条件为|u_i|≤2|u_j|，for alli,j∈J₀；其中，u_i和u_j分别是相关系数矩阵U中的第i列和第j列；

步骤8：根据Λ_m＝Λ_m-1∪J₀、A_m＝A_m-1∪a_j，j∈J₀，计算A_m；

步骤9：根据以下公式计算当前迭代次数对应的每个完备光声信号稀疏表示的系数；计算公式为：

向量

是完备光声信号稀疏表示的系数矩阵

的第i列元素；向量y_i是矩阵Y的第i列元素；向量θ_yi是矩阵Θ_Y的第i列元素

步骤10：计算下一次迭代余量；下一迭代余量的计算公式为

r_m+1为下一迭代余量；

步骤11：判断下一次迭代余量是否满足||r_m+1-r_m||≤ε₂，其中ε₂是控制阶段转换的阈值，若是则令stage＝stage+1，size＝size·stage，转向步骤6；若否，则令m＝m+1转向步骤4。

可选的，所述利用所述系数矩阵和所述自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵，具体包括：

采用公式

恢复出完备光声信号数据矩阵；其中，矩阵

是恢复出的完备光声信号数据矩阵，大小是N×L维；

是完备光声信号稀疏表示的系数矩阵；

是自适应完备字典。

可选的，所述根据所述完备光声信号数据矩阵，采用图像重建公式，重建光吸收能量分布图，具体包括：

将恢复出的完备光声信号数据矩阵

的每个元素

代入以下公式重建光吸收能量分布图；所述公式为

其中，k＝1,2,...,N，i＝1,2,...,L，

是

的第k行、第i列元素，Φ(r)是位置r处的光吸收能量，C_p是组织的比热容，c是超声波在组织中的传播速度，β是组织的体积膨胀温度系数，r₀是探测器与图像中心点之间的距离矢量，φ₀是探测器与x轴之间的夹角。

一种利用稀疏测量数据重建内窥式光声层析图像的***，包括：

自适应完备字典构造模块，用于构造可用于光声信号稀疏表示的自适应完备字典；

测量矩阵构建模块，用于构建测量矩阵；

感知矩阵计算模块，用于根据所述测量矩阵和所述自适应完备字典，计算感知矩阵；

稀疏光声信号测量数据矩阵构建模块，用于利用所述测量矩阵获取稀疏光声信号测量数据，并构建稀疏光声信号测量数据矩阵；

系数矩阵计算模块，用于根据所述稀疏光声信号测量数据矩阵和所述感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵；

完备光声信号数据计算模块，用于利用所述系数矩阵和所述自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵；

光吸收能量分布图重建模块，用于根据所述完备光声信号数据矩阵，采用图像重建公式，重建光吸收能量分布图。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明提供了一种利用稀疏测量数据重建内窥式光声层析图像的方法及***，解决了利用稀疏的光声测量数据重建光吸收分布图像中出现严重伪影和失真的问题，提高了成像精度。该方法首先构造可用于光声信号稀疏表示的自适应完备字典；其次利用随机测量矩阵测量到数据构建稀疏光声信号测量数据矩阵；然后根据测量矩阵和自适应完备字典计算感知矩阵；再者根据稀疏光声信号测量数据矩阵和感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵；最后利用系数矩阵和自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵，并采用图像重建公式，重建光吸收能量分布图。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例利用稀疏测量数据重建内窥式光声层析图像方法的流程示意图；

图2为本发明实施例内窥式PAT成像及图像重建示意图；其中，图2(a)是成像示意图，图2(b)是超声探测器接收组织产生的光声信号和图像重建示意图；

图3为本发明实施例利用稀疏测量数据重建内窥式光声层析图像***的结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明提供了一种利用稀疏测量数据重建内窥式光声层析图像的方法及***，以提高成像精度。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

文中各符号为：N为测量位置数；L为在各测量位置处采集的光声信号的长度；D₀为N×N维的初始完备字典；P为在N个测量位置处采集的完备光声信号数据矩阵，大小是N×L维；Θ₀为矩阵P基于D₀的初始稀疏表示的N×L维的系数矩阵；θ_0i为矩阵Θ₀的第i列元素，其中i＝1,2,...,L；θ₀ ^j为矩阵Θ₀的第j行元素，其中j＝1,2,...,N；θ₀ ^k为Θ₀的第k行元素，其中k＝1,2,...,N；T₀为稀疏表示系数中非零元素的个数；||·||₂、2-范数；||·||₀为0-范数；D为N×N维的更新后的完备字典；d_j为D的第j列元素，即第j个原子；d_k为D的第k列元素，即第k个原子；E_0k为去掉d_k后在P中造成的误差矩阵；E_0k'为E_0k去零后的结果矩阵；θ₀ ^k'为θ₀ ^k去零后的结果向量；E_0k'^T为E_0k'的转置矩阵；Q为由E_0k'^TE_0k'的单位特征向量组成的N×N维的正交矩阵；S为由E_0k'的奇异值组成的N×L维对角矩阵；S(1,1)为矩阵S的第一行、第一列位置上的元素，即E_0k'的最大奇异值；V为由E_0k'^TE_0k'的单位特征向量组成的L×L维正交矩阵；V^T为V的转置矩阵；Θ₀'为更新后的初始稀疏表示；Θ为矩阵P基于D的稀疏表示的N×L维系数矩阵；ε为允许的最大重构误差；

为自适应完备字典；Φ为M×N维的随机测量矩阵；A为M×N维的感知矩阵；M为稀疏测量位置数，且M<<N；α_k为A的第k列元素，即第k个原子；P_Y为N×L维完备光声信号数据矩阵；Y为M×L维的稀疏(即不完备)光声信号测量数据矩阵；y_i为Y的第i列元素；Θ_Y为完备光声信号稀疏表示的N×L维系数矩阵；θ_yi为Θ_Y的第i列元素；m为迭代次数；

为空集；r₀为信号的初始余量(残差)；r_m-1为第m-1次迭代的余量；r_m为第m次迭代的余量；size为搜索步长；stage为阶段；<·>为向量的内积运算；U为相关系数矩阵；u_k为U的第k列元素；J为U中size个最大值对应A的索引值(即列序号)集合；J₀为每次迭代找到的索引值(即列序号)集合；u_i为U的第i列元素，其中i∈J₀；u_j为U的第j列元素，其中j∈J₀；Λ_m-1为第m-1次迭代从A中选出的列向量的索引(即列序号)集合；Λ_m为第m次迭代从A中选出的列向量的索引(即列序号)集合；A_m-1为按索引Λ_m-1选出的A的列向量集合；A_m为按索引Λ_m选出的A的列向量集合；a_j为A的第j列元素，其中j∈J₀；ε₁为控制迭代次数的阈值；ε₂为控制阶段转换的阈值；

为完备光声信号稀疏表示的N×L维系数矩阵；

为

的第i列元素；

为恢复出的N×L维完备光声信号数据矩阵；

为矩阵

的第k行、第i列元素，其中k＝1,2,...,N，i＝1,2,...,L；r为成像点与图像中心点之间的距离矢量；Φ(r)为位置r处的光吸收能量；C_p为组织的比热容；c为超声波在组织中的传播速度；β为组织的体积膨胀温度系数；r₀为探测器与图像中心点之间的距离矢量；r_i为成像点与探测器之间的距离矢量；φ₀为探测器与x轴之间的夹角；γ为超声探测器的有效接收角度；x、y分别为x-y平面直角坐标系的横轴和纵轴。

实施例1

图1为本发明实施例利用稀疏测量数据重建内窥式光声层析图像方法的流程示意图，如图1所示，本发明实施例提供的方法具体包括以下步骤：

步骤101：构造可用于光声信号稀疏表示的自适应完备字典。

步骤102：构建测量矩阵。

步骤103：根据所述测量矩阵和所述自适应完备字典，计算感知矩阵。

步骤104：利用所述测量矩阵获取稀疏光声信号测量数据，并构建稀疏光声信号测量数据矩阵。

步骤105：根据所述稀疏光声信号测量数据矩阵和所述感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵。

步骤106：利用所述系数矩阵和所述自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵。

步骤107：根据所述完备光声信号数据矩阵，采用图像重建公式，重建光吸收能量分布图。

其中，步骤101具体包括以下步骤：

步骤1：建立完备光声信号测量数据集；

首先，在N个测量位置处获取完备的光声信号测量数据集，在各测量位置处采集的光声信号的长度为L，则完备光声信号测量数据集表示为N×L维的矩阵P。

步骤2：初始化完备字典；

采用离散余弦变换基作为初始完备字典D₀，其矩阵形式为：

其中，矩阵D₀是N×N维的初始完备字典，N是测量位置的总数。

步骤3：完备光声信号测量数据集的稀疏表示；

计算矩阵P基于初始完备字典D₀的初始稀疏表示的系数矩阵Θ₀；其计算公式为

其中，矩阵Θ₀是P的基于初始字典D₀的初始稀疏表示的系数矩阵，大小是N×L维；向量θ_0i是矩阵Θ₀的第i列元素，其中i＝1,2,...,L；T₀为稀疏表示系数中非零元素的个数；||·||₂表示2-范数；||·||₀表示0-范数。

步骤4：更新完备字典；

将矩阵P作为训练样本，采用以下公式计算更新后的完备字典D。其计算公式为：

其中，矩阵D是更新后的N×N维的完备字典；向量d_j和向量d_k分别是矩阵D的第j列和第k列元素，即第j个原子和第k个原子；||·||₂表示2-范数；向量θ₀ ^j和θ₀ ^k分别是矩阵Θ₀的第j行和第k行元素。

本发明实施例选择了奇异值分解法来进行具体求解，其目的作用是采用奇异值分解法得到更新后的d_k，进而得到更新后的完备字典D。具体为：对误差矩阵进行奇异值分解。

去掉θ₀ ^k中所有的0元素，仅保留非零值，所得结果向量为θ₀ ^k'。接着，得到误差矩阵：

式中，E_0k是去掉d_k后在P中造成的误差矩阵，大小是N×L维。仅保留E_0k中d_k与θ₀ ^k的非零值乘积项所在的列得到E_0k'，大小是N×L维。最后，对E_0k'进行奇异值分解：E_0k'＝QSV^T(5)；

其中，Q是由E_0k'^TE_0k'的单位特征向量组成的N×N维的正交矩阵，其中第一列元素是向量d_k；S是由E_0k'的奇异值组成的N×L维对角矩阵；V是由E_0k'^TE_0k'的单位特征向量组成的L×L维正交矩阵；V^T是V的转置矩阵。用V的第一列和S(1,1)的乘积更新θ₀ ^k'，其中S(1,1)是矩阵S的第一行、第一列位置上的元素，即E_0k'的最大奇异值。在逐列更新完成后，得到更新后的字典D和更新后的初始稀疏表示的系数矩阵Θ₀'；

步骤5：根据公式(6)判断是否收敛；公式(6)为

其中，ε是允许的最大重构误差。

若矩阵P、矩阵D和矩阵Θ₀'满足公式(6)的收敛条件，则停止迭代，将自适应完备字典

否则转向步骤3，用更新后的字典D代替初始完备字典D₀，更新稀疏表示。

步骤102具体包括：构造M×N维的测量矩阵Φ。

具体方法是：首先，生成N×N维的哈达玛矩阵；然后，在哈达玛矩阵中随机地选取M行向量，构成一个M×N维的矩阵，其中M<<N，M是稀疏测量位置数。

步骤103具体包括：

根据测量矩阵Φ和自适应完备字典

得到M×N维的感知矩阵A。

其中，感知矩阵A为

A满足限定等距性(Restricted IsometryProperty,RIP)条件。

步骤104具体包括：

超声换能器根据测量矩阵Φ在M个测量位置处测得稀疏的(即不完备的)光声信号测量数据，并构建稀疏光声信号测量数据矩阵Y。

稀疏光声信号测量数据矩阵Y为Y＝Φ·P_Y(8)；

其中，矩阵Y是M×L维的稀疏(即不完备)光声信号测量数据矩阵；矩阵P_Y是N×L维的完备光声信号矩阵。

步骤105具体包括以下几个步骤：

步骤1：根据稀疏光声信号测量数据矩阵Y和感知矩阵A，得到完备光声信号矩阵P_Y的稀疏表示的系数矩阵Θ_Y。

其中，

s.t.Y＝A·Θ_Y(9)；式中||·||₀是矩阵的0-范数；Θ_Y是完备光声信号稀疏表示的系数矩阵，大小是N×L维。

步骤2：初始化参数；r₀＝y_i，size＝4，stage＝1，m＝1，

其中各参数的含义为：向量y_i是稀疏光声信号测量数据矩阵Y的第i列元素；r₀初次余量(残差)，r_m是第m次迭代的余量(残差)；size是搜索步长；stage是阶段；m是迭代次数；U是相关系数矩阵；J是相关系数矩阵U按照从大到小排列，前size个值对应感知矩阵A的索引值(列序号)集合；J₀是每次迭代找到的索引值(列序号)集合；Λ_m是第m次迭代从感知矩阵A中选出的列向量的索引(列序号)集合；A_m是按索引Λ_m选出的感知矩阵A的列向量集合；

是空集。

步骤3：计算当前迭代余量r_m。其计算公式为：

步骤4：判断当前迭代余量r_m是否满足||r_m||₂≤ε₁；其中ε₁是控制迭代次数的阈值，若是，则停止迭代，并根据上次迭代确定的每个完备光声信号稀疏表示的系数构建系数矩阵

若否，则转向步骤5。

步骤5：计算相关系数矩阵U；所述相关系数矩阵U的计算公式为U＝{u_k|u_k＝<r_m-1,α_k>|,k＝1,2,3,…,N}(10)；其中，向量u_k是相关系数矩阵U的第k列，向量α_k是感知矩阵A的第k列，<·>是向量的内积运算。

步骤6：将相关系数矩阵U中的元素按照从大到小进行排列，选择前size个元素对应感知矩阵A的列序号来构成集合J。

步骤7：正则化；在集合J中寻找满足约束条件的子集J₀，并在所有满足约束条件的子集J₀中选择具有最大能量的子集；所述约束条件为|u_i|≤2|u_j|，for alli,j∈J₀(11)；其中，u_i和u_j分别是相关系数矩阵U中的第i列和第j列。

能量计算公式为

步骤8：令Λ_m＝Λ_m-1∪J₀，A_m＝A_m-1∪a_j(j∈J₀)，计算A_m。

步骤9：根据以下公式计算每个完备光声信号稀疏表示的系数；其计算公式为：

其中，向量

是完备光声信号稀疏表示的系数矩阵

的第i列元素；向量y_i是矩阵Y的第i列元素；向量θ_yi是矩阵Θ_Y的第i列元素。

步骤10：计算下一次迭代的余量；下一迭代余量的计算公式为

步骤11：判断下一次迭代余量是否满足||r_m+1-r_m||≤ε₂，其中ε₂是控制阶段转换的阈值，若是则更新stage和size，即stage＝stage+1，size＝size·stage，转向步骤6；若否则，令m＝m+1转向步骤4。

步骤106具体包括利用完备光声信号稀疏表示的系数矩阵

和自适应完备字典

通过稀疏反变换恢复出完备光声信号数据矩阵。

其计算公式为

其中，矩阵

是恢复出的完备光声信号数据矩阵，大小是N×L维。

步骤107具体包括：将矩阵

的每个元素

代入下式重建光吸收能量分布图。

其中，k＝1,2,...,N，i＝1,2,...,L，

是

的第k行、第i列元素，Φ(r)是位置r处的光吸收能量，C_p是组织的比热容，c是超声波在组织中的传播速度，β是组织的体积膨胀温度系数，r₀是探测器与图像中心点之间的距离矢量，φ₀是探测器与x轴之间的夹角，如图2所示。

实施例2

为了实现上述目的，本发明还提供了一种利用稀疏测量数据重建内窥式光声层析图像的***。

图3为本发明实施例利用稀疏测量数据重建内窥式光声层析图像***的结构示意图，如图3所示，本发明实施例提供的***包括以下几个模块。

自适应完备字典构造模块100，用于构造可用于光声信号稀疏表示的自适应完备字典。

测量矩阵构建模块200，用于构建测量矩阵。

感知矩阵计算模块300，用于根据所述测量矩阵和所述自适应完备字典，计算感知矩阵。

稀疏光声信号测量数据矩阵构建模块400，用于利用所述测量矩阵获取稀疏光声信号测量数据，并构建稀疏光声信号测量数据矩阵。

系数矩阵计算模块500，用于根据所述稀疏光声信号测量数据矩阵和所述感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵。

完备光声信号数据计算模块600，用于利用所述系数矩阵和所述自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵。

光吸收能量分布图重建模块700，用于根据所述完备光声信号数据矩阵，采用图像重建公式，重建光吸收能量分布图。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的***而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (7)

1.一种利用稀疏测量数据重建内窥式光声层析图像的方法，其特征在于，所述方法包括：

构造可用于光声信号稀疏表示的自适应完备字典；

构建测量矩阵；

根据所述测量矩阵和所述自适应完备字典，计算感知矩阵；

利用所述测量矩阵获取稀疏光声信号测量数据，并构建稀疏光声信号测量数据矩阵；

根据所述稀疏光声信号测量数据矩阵和所述感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵；

利用所述系数矩阵和所述自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵；

根据所述完备光声信号数据矩阵，采用图像重建公式，重建光吸收能量分布图；

所述构造可用于光声信号稀疏表示的自适应完备字典，具体包括：

构建初始完备光声信号测量数据矩阵；所述初始完备光声信号测量数据矩阵中的元素为在N个测量位置处获取的完备光声信号测量数据，且在各测量位置处采集的光声信号的长度为L；所述初始完备光声信号测量数据矩阵为N×L维的矩阵P；

采用离散余弦变换基作为初始完备字典D₀；所述初始完备字典D₀是N×N维的初始完备字典，N是测量位置总数；

以所述初始完备光声信号测量数据矩阵为训练样本，按照公式

计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'；

其中，D是更新后的N×N维的完备字典；矩阵Θ₀是矩阵P基于初始完备字典D₀的初始稀疏表示的系数矩阵，大小是N×L维；向量d_j和d_k分别是更新后的完备字典D的第j列和第k列元素，即第j个原子和第k个原子；||·||₂表示2-范数；向量θ₀ ^j和θ₀ ^k分别是矩阵Θ₀的第j行和第k行元素；向量θ_0i是矩阵Θ₀的第i列元素，其中i＝1,2,...,L；

计算矩阵P基于更新后的完备字典D的稀疏表示的系数矩阵Θ；系数矩阵Θ的大小是N×L维；

判断矩阵P更新后的完备字典D以及系数矩阵Θ是否满足

ε为允许的最大重构误差；

若是则停止迭代，并将当前迭代对应的更新后的完备字典D确定为自适应完备字典；

若否则将更新后的初始稀疏表示的系数矩阵Θ₀'代替初始稀疏表示的系数矩阵矩阵Θ₀，返回计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'步骤。

2.根据权利要求1所述的利用稀疏测量数据重建内窥式光声层析图像的方法，其特征在于，所述计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'，具体包括：

以所述初始完备光声信号测量数据矩阵为训练样本，按照公式

采用奇异值分解算法，计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'。

3.根据权利要求1所述的利用稀疏测量数据重建内窥式光声层析图像的方法，其特征在于，所述构建测量矩阵，具体包括：

生成N×N维的哈达玛矩阵；

在所述哈达玛矩阵中随机地选取M行向量，构建测量矩阵；所述构建测量矩阵为一个M×N维的矩阵，其中M<<N，M是稀疏测量位置数。

4.根据权利要求1所述的利用稀疏测量数据重建内窥式光声层析图像的方法，其特征在于，所述根据所述稀疏光声信号测量数据矩阵和所述感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵，具体包括：

步骤1：根据稀疏光声信号测量数据矩阵Y和感知矩阵A，得到完备光声信号矩阵P_Y的稀疏表示的系数矩阵Θ_Y；

步骤2：初始化参数；r₀＝y_i，size＝4，stage＝1，m＝1，

其中各参数的含义为：向量y_i是稀疏光声信号测量数据矩阵Y的第i列元素；r₀是初次余量，r_m是第m次迭代的余量；size是搜索步长；stage是阶段；m是迭代次数；U是相关系数矩阵；J是相关系数矩阵U按照从大到小排列，前size个值对应感知矩阵A的列序号集合；J₀是每次迭代找到的列序号集合；Λ_m是第m次迭代从感知矩阵A中选出的列向量的列序号集合；A_m是按索引Λ_m选出的感知矩阵A的列向量集合；

是空集；

步骤3：计算当前迭代余量r_m；

步骤4：判断当前迭代余量r_m是否满足||r_m||₂≤ε₁；其中ε₁是控制迭代次数的阈值，若是，则停止迭代，并根据上次迭代确定的每个完备光声信号稀疏表示的系数构建系数矩阵

若否，则转向步骤5；

步骤5：计算相关系数矩阵U；所述相关系数矩阵U的计算公式为U＝{u_k|u_k＝<r_m-1,α_k>|,k＝1,2,3,…,N}；其中，向量u_k是相关系数矩阵U的第k列，向量α_k是感知矩阵A的第k列，<·>是向量的内积运算；

步骤6：将相关系数矩阵U中的元素按照从大到小进行排列，选择前size个元素对应感知矩阵A的列序号来构成集合J；

步骤7：正则化；在集合J中寻找满足约束条件的子集J₀，并在所有满足约束条件的子集J₀中选择具有最大能量的子集；所述约束条件为|u_i|≤2|u_j|，for alli,j∈J₀；其中，u_i和u_j分别是相关系数矩阵U中的第i列和第j列；

步骤8：根据Λ_m＝Λ_m-1∪J₀、A_m＝A_m-1∪a_j，j∈J₀，计算A_m；a_j为感知矩阵A的第j列元素；

步骤9：根据以下公式计算当前迭代次数对应的每个完备光声信号稀疏表示的系数；计算公式为：

向量

是完备光声信号稀疏表示的系数矩阵

的第i列元素；向量y_i是矩阵Y的第i列元素；向量θ_yi是矩阵Θ_Y的第i列元素

步骤10：计算下一次迭代余量；下一迭代余量的计算公式为

r_m+1为下一迭代余量；

步骤11：判断下一次迭代余量是否满足||r_m+1-r_m||≤ε₂，其中ε₂是控制阶段转换的阈值，若是则令stage＝stage+1，size＝size·stage，转向步骤6；若否，则令m＝m+1转向步骤4。

5.根据权利要求1所述的利用稀疏测量数据重建内窥式光声层析图像的方法，其特征在于，所述利用所述系数矩阵和所述自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵，具体包括：

采用公式

恢复出完备光声信号数据矩阵；其中，矩阵

是恢复出的完备光声信号数据矩阵，大小是N×L维；

是完备光声信号稀疏表示的系数矩阵；

是自适应完备字典。

6.根据权利要求1所述的利用稀疏测量数据重建内窥式光声层析图像的方法，其特征在于，所述根据所述完备光声信号数据矩阵，采用图像重建公式，重建光吸收能量分布图，具体包括：

将恢复出的完备光声信号数据矩阵

的每个元素

代入以下公式重建光吸收能量分布图；所述公式为

其中，k＝1,2,...,N，i＝1,2,...,L，

是

的第k行、第i列元素，Φ(r)是位置r处的光吸收能量，C_p是组织的比热容，c是超声波在组织中的传播速度，β是组织的体积膨胀温度系数，r₀是探测器与图像中心点之间的距离矢量，φ₀是探测器与x轴之间的夹角。

7.一种利用稀疏测量数据重建内窥式光声层析图像的***，其特征在于，所述***包括：

自适应完备字典构造模块，用于构造可用于光声信号稀疏表示的自适应完备字典；

测量矩阵构建模块，用于构建测量矩阵；

感知矩阵计算模块，用于根据所述测量矩阵和所述自适应完备字典，计算感知矩阵；

稀疏光声信号测量数据矩阵构建模块，用于利用所述测量矩阵获取稀疏光声信号测量数据，并构建稀疏光声信号测量数据矩阵；

系数矩阵计算模块，用于根据所述稀疏光声信号测量数据矩阵和所述感知矩阵，优化计算得到完备光声信号稀疏表示的系数矩阵；

完备光声信号数据计算模块，用于利用所述系数矩阵和所述自适应完备字典，通过稀疏反变换恢复出完备光声信号数据矩阵；

光吸收能量分布图重建模块，用于根据所述完备光声信号数据矩阵，采用图像重建公式，重建光吸收能量分布图；

所述构造可用于光声信号稀疏表示的自适应完备字典，具体包括：

构建初始完备光声信号测量数据矩阵；所述初始完备光声信号测量数据矩阵中的元素为在N个测量位置处获取的完备光声信号测量数据，且在各测量位置处采集的光声信号的长度为L；所述初始完备光声信号测量数据矩阵为N×L维的矩阵P；

采用离散余弦变换基作为初始完备字典D₀；所述初始完备字典D₀是N×N维的初始完备字典，N是测量位置总数；

以所述初始完备光声信号测量数据矩阵为训练样本，按照公式

计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'；

其中，D是更新后的N×N维的完备字典；矩阵Θ₀是矩阵P基于初始完备字典D₀的初始稀疏表示的系数矩阵，大小是N×L维；向量d_j和d_k分别是更新后的完备字典D的第j列和第k列元素，即第j个原子和第k个原子；||·||₂表示2-范数；向量θ₀ ^j和θ₀ ^k分别是矩阵Θ₀的第j行和第k行元素；向量θ_0i是矩阵Θ₀的第i列元素，其中i＝1,2,...,L；

计算矩阵P基于更新后的完备字典D的稀疏表示的系数矩阵Θ；系数矩阵Θ的大小是N×L维；

判断矩阵P更新后的完备字典D以及系数矩阵Θ是否满足

ε为允许的最大重构误差；

若是则停止迭代，并将当前迭代对应的更新后的完备字典D确定为自适应完备字典；

若否则将更新后的初始稀疏表示的系数矩阵Θ₀'代替初始稀疏表示的系数矩阵矩阵Θ₀，返回计算更新后的完备字典D和更新后的初始稀疏表示的系数矩阵Θ₀'步骤。

Images