一种光吸收能量分布图和声速分布图的重建方法及系统

1.本发明涉及血管内光声层析成像技术领域，特别是涉及一种光吸收能量分布图和声速分布图的重建方法及系统。


背景技术：

2.光声层析成像(photoacoustic tomography，pat)是一种新型复合功能成像技术，血管内光声层析成像(intravascular photoacoustic tomography，ivpat)是内窥式pat的典型应用，也是血管内超声的补充成像手段，对于粥样硬化性心血管疾病的临床诊断和介入治疗具有重要的研究意义和广阔的应用前景。
3.高质量的图像重建是ivpat技术的关键，标准重建方法(如滤波反投影和时间反演等)通常基于理想成像条件的假设，例如：超声波在组织中的传播速度恒定、探测器具有理想特性和足够带宽且能够获得全视角的完备测量数据、光声信号的时间和空间采样率足够高等。但在实际应用中，由于血管的复杂结构和成像装置等因素的限制，超声探测器往往只能在有限视角内采集不完备的光声信号。此外，考虑到血管组织复杂的光学和声学特性，超声波在不同组织成分中的传播速度通常存在较大差异，而多数情况下，在进行成像前无法预知组织中声速的空间分布情况，因此在重建图像的过程中，恒定声速的假设会导致重建图像的质量下降。


技术实现要素：

4.本发明的目的是提供一种光吸收能量分布图和声速分布图的重建方法及系统，能够根据超声探测器在血管腔内采集的有限角度光声信号，同时重建血管横截面上的光吸收能量分布图和声速分布图，并提高重建图像的质量。
5.为实现上述目的，本发明提供了如下方案：
6.一种光吸收能量分布图和声速分布图的重建方法，包括：
7.获取待成像平面的测量成像数据；所述测量成像数据包括有限角度光声信号测量数据和常数声速值；
8.将所述测量成像数据输入到图像重建模型中，得到待成像平面处的光吸收能量重建分布图和声速重建分布图；所述图像重建模型是根据历史成像数据，利用梯度下降算法对卷积神经网络进行训练得到的。
9.可选的，在所述获取待成像平面的成像数据之前，还包括：
10.获取训练集；所述训练集包括待测血管多个平面的历史成像数据、光吸收能量历史分布图和声速历史分布图；
11.以所述历史成像数据为输入，以所述光吸收能量历史分布图和所述声速历史分布图为期望输出，利用梯度下降算法对卷积神经网络进行训练，得到图像重建模型。
12.可选的，所述以所述历史成像数据为输入，以所述光吸收能量历史分布图和所述声速历史分布图为期望输出，利用梯度下降算法对卷积神经网络进行训练，得到图像重建
模型，具体包括：
13.将所述训练集划分为多个训练子集；
14.确定训练前的卷积神经网络为第i
×
epoch
×
i阶卷积神经网络；
15.令网络层序号i＝1；
16.令训练集遍历次数epoch＝0；
17.确定任一训练子集为当前训练子集；
18.令当前训练子集的索引i＝1；
19.将所述当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
(i 1)阶卷积神经网络；
20.更新当前训练子集，令i的数值增加1，并返回步骤“将所述当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
(i 1)阶卷积神经网络”直至历遍训练子集，得到第i
×
(epoch 1)
×
1阶卷积神经网络；
21.令epoch的数值增加1并返回步骤“令当前训练子集的索引i＝1”直至epoch大于训练集遍历次数阈值，得到第(i 1)
×0×
1阶卷积神经网络；
22.令i的数值增加1并返回步骤“令训练集遍历次数epoch＝0”直至i大于网络层序号阈值，得到i
×
epoch
×
i阶卷积神经网络作为图像重建模型。
23.可选的，所述将所述当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
(i 1)阶卷积神经网络，具体包括：
24.将所述当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
i次迭代时的光吸收能量历史预测分布图和第i
×
epoch
×
i次迭代时的声速历史预测分布图；
25.根据光吸收能量历史分布图和第i
×
epoch
×
i次迭代时的光吸收能量历史预测分布图，确定第i
×
epoch
×
i阶卷积神经网络执行光吸收能量分布时的损失为第一损失；
26.根据声速历史分布图和第i
×
epoch
×
i次迭代时的声速历史预测分布图，确定第i
×
epoch
×
i阶卷积神经网络执行声速分布时的损失为第二损失；
27.根据所述第一损失确定第i
×
epoch
×
i阶卷积神经网络输出层的第一误差向量；
28.根据所述第二损失确定第i
×
epoch
×
i阶卷积神经网络输出层的第二误差向量；
29.根据所述第一误差向量和第一递推公式确定第i
×
epoch
×
i阶卷积神经网络每一层的第一误差向量；
30.根据所述第二误差向量和第二递推公式确定第i
×
epoch
×
i阶卷积神经网络每一层的第二误差向量；
31.根据第一权重、第二权重、多个第一误差向量和多个第二误差向量，确定第i
×
epoch
×
i阶卷积神经网络每一层参数的梯度；所述第一权重为卷积神经网络执行光吸收能量分布时的权重；所述第二权重为卷积神经网络执行声速分布时的权重；
32.基于第i
×
epoch
×
i阶卷积神经网络每一层参数的梯度，利用adam算法调整第i
×
epoch
×
i阶卷积神经网络中每一层参数的尺度，得到第i
×
epoch
×
(i 1)阶卷积神经网络。
33.可选的，在所述利用adam算法调整第i
×
epoch
×
i阶卷积神经网络中每一层参数的尺度之后，还包括：
34.利用梯度正则化算法调整第一权重和第二权重。
35.可选的，
36.所述第一损失为：
37.所述第二损失为：
38.所述卷积神经网络输出层的第一误差向量为：
39.所述卷积神经网络输出层的第二误差向量为：
40.所述第一递推公式为：δ
1,c
＝(w
1,c 1
)
t
δ
1,c 1
⊙
σ'(z
1,c
)；
41.所述第二递推公式为：δ
2,c
＝(w
2,c 1
)
t
δ
2,c 1
⊙
σ'(z
2,c
)；
42.卷积神经网络每一层参数的梯度：
[0043][0044]
其中，c1(w,b)和c2(w,b)分别为第一损失和第二损失；w为网络的权重参数集；b为网络的偏置参数集；k为第i组训练子集中样本的序号，k＝1,2,3,...,m；m为当前训练子集的样本数量；a
true,k
和c
true,k
分别为第k个样本的光吸收能量历史分布图和声速历史分布图；a
i,k
和c
i,k
分别为第k个样本经过第i次迭代后的光吸收能量历史预测分布图和声速历史预测分布图；θ
1,i
是执行光吸收能量分布的第i个网络结构单元的参数，是针对执行光吸收能量分布的结构相同且具有不同学习参数θ
1,i
的卷积神经网络，是第i个网络结构单元中光吸收能量分布的梯度，pk为历史成像数据中的第k个样本的有限角度光声信号测量数据，a
i-1,k
为第k个样本经过第i-1次迭代后的光吸收能量历史预测分布图；θ
2,i
是执行声速分布的第i个网络结构单元的参数，是针对执行声速分布时结构相同且具有不同学习参数θ
2,i
的卷积神经网络，是第i个网络结构单元中声速分布的梯度，c
i-1,k
为第k个样本经过第i-1次迭代后的声速历史预测分布图，h(c
i-1,k
)是与c
i-1,k
有关的算子；δ
1,21
和δ
2,21
分别是卷积神经网络输出层的第一误差向量和第二误差向量；z
1,21
和z
2,21
分别是执行光吸收能量分布和执行声速分布时第21层网络产生的线性激活响应；第21层网络为卷积神经网络输出层；
“⊙”
是矩阵的hadamard积，表示矩阵对应元素点对点的乘法；σ'(
·
)为relu函数的一阶导数；δ
1,c
和δ
2,c
分别是执行光吸收能量分布和执行声速分布时第c层网络的第一误差向量和第二误差向量；c为特征图层序号，在第一递推公式、第二
递推公式中c＝20,19,...,3,2，在求卷积神经网络每一层参数的梯度的公式中c＝2,3,...,21；w
1,c 1
和w
2,c 1
分别为执行光吸收能量分布和执行声速分布时第c 1层网络的权重参数矩阵；上标“t”表示求矩阵的转置；δ
1,c 1
和δ
2,c 1
分别是执行光吸收能量分布和执行声速分布时第c 1层网络的第一误差向量和第二误差向量；z
1,c
和z
2,c
分别为前向传播中执行光吸收能量分布和执行声速分布时第c层网络产生的线性激活响应；c(w,b)为总损失；λ1和λ2分别为第一权重和第二权重；w
1,c
和w
2,c
分别为执行光吸收能量分布和执行声速分布时第c层网络的权重参数矩阵；b
1,c
和b
2,c
分别为执行光吸收能量分布和执行声速分布时第c层网络的偏置参数向量；a
1,c-1
和a
2,c-1
分别为执行光吸收能量分布和执行声速分布时第c-1层网络的输出；wc和bc分别是第c层网络的权重参数矩阵和偏置参数向量。
[0045]
一种光吸收能量分布图和声速分布图的重建系统，包括：
[0046]
测量成像数据获取模块，用于获取待成像平面的测量成像数据；所述测量成像数据包括有限角度光声信号测量数据和常数声速值；
[0047]
分布图重建模块，用于将所述测量成像数据输入到图像重建模型中，得到待成像平面处的光吸收能量重建分布图和声速重建分布图；所述图像重建模型是根据历史成像数据，利用梯度下降算法对卷积神经网络进行训练得到的。
[0048]
可选的，所述系统，还包括：
[0049]
训练集获取模块，用于获取训练集；所述训练集包括待测血管多个平面的历史成像数据、光吸收能量历史分布图和声速历史分布图；
[0050]
图像重建模型确定模块，用于以所述历史成像数据为输入，以所述光吸收能量历史分布图和所述声速历史分布图为期望输出，利用梯度下降算法对卷积神经网络进行训练，得到图像重建模型。
[0051]
可选的，所述图像重建模型确定模块，具体包括：
[0052]
训练子划分单元，用于将所述训练集划分为多个训练子集；
[0053]
初始卷积神经网络确定单元，用于确定训练前的卷积神经网络为第i
×
epoch
×
i阶卷积神经网络；
[0054]
网络层序号赋值单元，用于令网络层序号i＝1；
[0055]
训练集遍历次数赋值单元，用于令训练集遍历次数epoch＝0；
[0056]
当前训练子集确定单元，用于确定任一训练子集为当前训练子集；
[0057]
索引赋值单元，用于令当前训练子集的索引i＝1；
[0058]
训练单元，用于将所述当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
(i 1)阶卷积神经网络；
[0059]
当前训练子集更新单元，用于更新当前训练子集，令i的数值增加1，并调用所述训练单元直至历遍训练子集，得到第i
×
(epoch 1)
×
1阶卷积神经网络；
[0060]
训练集遍历次数更新单元，用于令epoch的数值增加1并调用所述索引赋值单元直至epoch大于训练集遍历次数阈值，得到第(i 1)
×0×
1阶卷积神经网络；
[0061]
网络层序号更新单元，用于令i的数值增加1并调用所述训练集遍历次数赋值单元直至i大于网络层序号阈值，得到i
×
epoch
×
i阶卷积神经网络作为图像重建模型。
[0062]
可选的，所述训练单元，具体包括：
[0063]
历史预测分布图确定子单元，用于将所述当前训练子集输入到第i
×
epoch
×
i阶
卷积神经网络中，得到第i
×
epoch
×
i次迭代时的光吸收能量历史预测分布图和第i
×
epoch
×
i次迭代时的声速历史预测分布图；
[0064]
第一损失确定子单元，用于根据光吸收能量历史分布图和第i
×
epoch
×
i次迭代时的光吸收能量历史预测分布图，确定第i
×
epoch
×
i阶卷积神经网络执行光吸收能量分布时的损失为第一损失；
[0065]
第二损失确定子单元，用于根据声速历史分布图和第i
×
epoch
×
i次迭代时的声速历史预测分布图，确定第i
×
epoch
×
i阶卷积神经网络执行声速分布时的损失为第二损失；
[0066]
第一误差向量确定子单元，用于根据所述第一损失确定第i
×
epoch
×
i阶卷积神经网络输出层的第一误差向量；
[0067]
第二误差向量确定子单元，用于根据所述第二损失确定第i
×
epoch
×
i阶卷积神经网络输出层的第二误差向量；
[0068]
第一误差向量组确定子单元，用于根据所述第一误差向量和第一递推公式确定第i
×
epoch
×
i阶卷积神经网络每一层的第一误差向量；
[0069]
第二误差向量组确定子单元，用于根据所述第二误差向量和第二递推公式确定第i
×
epoch
×
i阶卷积神经网络每一层的第二误差向量；
[0070]
梯度确定子单元，用于根据第一权重、第二权重、多个第一误差向量和多个第二误差向量，确定第i
×
epoch
×
i阶卷积神经网络每一层参数的梯度；所述第一权重为卷积神经网络执行光吸收能量分布时的权重；所述第二权重为卷积神经网络执行声速分布时的权重；
[0071]
参数尺度调整子单元，用于基于第i
×
epoch
×
i阶卷积神经网络每一层参数的梯度，利用adam算法调整第i
×
epoch
×
i阶卷积神经网络中每一层参数的尺度，得到第i
×
epoch
×
(i 1)阶卷积神经网络。
[0072]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0073]
本发明提供的一种光吸收能量分布图和声速分布图的重建方法及系统，利用梯度下降算法对卷积神经网络进行训练得到图像重建模型，能够根据超声探测器在血管腔内采集的有限角度光声信号，同时重建血管横截面上的光吸收能量分布图和声速分布图，并提高重建图像的质量。
附图说明
[0074]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0075]
图1为本发明实施例中光吸收能量分布图和声速分布图的重建方法的流程图。
[0076]
图2为本发明实施例中ivpat成像原理示意图。图2(a)为本发明实施例中超声探测器在血管腔内的测量位置示意图，图2(b)为本发明实施例中成像平面示意图；
[0077]
图3为本发明实施例中前向传播jir-net的结构示意图；
[0078]
图4为本发明实施例中搭建jir-net的流程图；
[0079]
图5为本发明实施例中训练jir-net的流程图。
具体实施方式
[0080]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0081]
本发明的目的是提供一种光吸收能量分布图和声速分布图的重建方法及系统，能够根据超声探测器在血管腔内采集的有限角度光声信号，同时重建血管横截面上的光吸收能量分布图和声速分布图，并提高重建图像的质量。
[0082]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0083]
如图1，本发明提供了一种光吸收能量分布图和声速分布图的重建方法，包括：
[0084]
获取待成像平面的测量成像数据；测量成像数据包括有限角度光声信号测量数据和常数声速值；
[0085]
将测量成像数据输入到图像重建模型中，得到待成像平面处的光吸收能量重建分布图和声速重建分布图；图像重建模型是根据历史成像数据，利用梯度下降算法对卷积神经网络进行训练得到的。
[0086]
在获取待成像平面的成像数据之前，还包括：
[0087]
获取训练集；训练集包括待测血管多个平面的历史成像数据、光吸收能量历史分布图和声速历史分布图；
[0088]
以历史成像数据为输入，以光吸收能量历史分布图和声速历史分布图为期望输出，利用梯度下降算法对卷积神经网络进行训练，得到图像重建模型。
[0089]
具体的，以历史成像数据为输入，以光吸收能量历史分布图和声速历史分布图为期望输出，利用梯度下降算法对卷积神经网络进行训练，得到图像重建模型，具体包括：
[0090]
将训练集划分为多个训练子集；
[0091]
确定训练前的卷积神经网络为第i
×
epoch
×
i阶卷积神经网络；
[0092]
令网络层序号i＝1；
[0093]
令训练集遍历次数epoch＝0；
[0094]
确定任一训练子集为当前训练子集；
[0095]
令当前训练子集的索引i＝1；
[0096]
将当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
(i 1)阶卷积神经网络；
[0097]
更新当前训练子集，令i的数值增加1，并返回步骤“将当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
(i 1)阶卷积神经网络”直至历遍训练子集，得到第i
×
(epoch 1)
×
1阶卷积神经网络；
[0098]
令epoch的数值增加1并返回步骤“令当前训练子集的索引i＝1”直至epoch大于训练集遍历次数阈值，得到第(i 1)
×0×
1阶卷积神经网络；
[0099]
令i的数值增加1并返回步骤“令训练集遍历次数epoch＝0”直至i大于网络层序号
阈值，得到i
×
epoch
×
i阶卷积神经网络作为图像重建模型。
[0100]
进一步的，将当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
(i 1)阶卷积神经网络，具体包括：
[0101]
将当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
i次迭代时的光吸收能量历史预测分布图和第i
×
epoch
×
i次迭代时的声速历史预测分布图；
[0102]
根据光吸收能量历史分布图和第i
×
epoch
×
i次迭代时的光吸收能量历史预测分布图，确定第i
×
epoch
×
i阶卷积神经网络执行光吸收能量分布时的损失为第一损失；
[0103]
根据声速历史分布图和第i
×
epoch
×
i次迭代时的声速历史预测分布图，确定第i
×
epoch
×
i阶卷积神经网络执行声速分布时的损失为第二损失；
[0104]
根据第一损失确定第i
×
epoch
×
i阶卷积神经网络输出层的第一误差向量；
[0105]
根据第二损失确定第i
×
epoch
×
i阶卷积神经网络输出层的第二误差向量；
[0106]
根据第一误差向量和第一递推公式确定第i
×
epoch
×
i阶卷积神经网络每一层的第一误差向量；
[0107]
根据第二误差向量和第二递推公式确定第i
×
epoch
×
i阶卷积神经网络每一层的第二误差向量；
[0108]
根据第一权重、第二权重、多个第一误差向量和多个第二误差向量，确定第i
×
epoch
×
i阶卷积神经网络每一层参数的梯度；第一权重为卷积神经网络执行光吸收能量分布时的权重；第二权重为卷积神经网络执行声速分布时的权重；
[0109]
基于第i
×
epoch
×
i阶卷积神经网络每一层参数的梯度，利用adam算法调整第i
×
epoch
×
i阶卷积神经网络中每一层参数的尺度，得到第i
×
epoch
×
(i 1)阶卷积神经网络。
[0110]
本发明提供的的光吸收能量分布图和声速分布图的重建方法，在利用adam算法调整第i
×
epoch
×
i阶卷积神经网络中每一层参数的尺度之后，还包括：
[0111]
利用梯度正则化算法调整第一权重和第二权重。
[0112]
其中，
[0113]
第一损失为：
[0114]
第二损失为：
[0115]
卷积神经网络输出层的第一误差向量为：
[0116]
卷积神经网络输出层的第二误差向量为：
[0117]
第一递推公式为：δ
1,c
＝(w
1,c 1
)
t
δ
1,c 1
⊙
σ'(z
1,c
)；
[0118]
第二递推公式为：δ
2,c
＝(w
2,c 1
)
t
δ
2,c 1
⊙
σ'(z
2,c
)；
[0119]
卷积神经网络每一层参数的梯度：
[0120][0121]
其中，c1(w,b)和c2(w,b)分别为第一损失和第二损失；w为网络的权重参数集；b为网络的偏置参数集；k为第i组训练子集中样本的序号，k＝1,2,3,...,m；m为当前训练子集的样本数量；a
true,k
和c
true,k
分别为第k个样本的光吸收能量历史分布图和声速历史分布图；a
i,k
和c
i,k
分别为第k个样本经过第i次迭代后的光吸收能量历史预测分布图和声速历史预测分布图；θ
1,i
是执行光吸收能量分布的第i个网络结构单元的参数，是针对执行光吸收能量分布的结构相同且具有不同学习参数θ
1,i
的卷积神经网络，是第i个网络结构单元中光吸收能量分布的梯度，pk为历史成像数据中的第k个样本的有限角度光声信号测量数据，a
i-1,k
为第k个样本经过第i-1次迭代后的光吸收能量历史预测分布图；θ
2,i
是执行声速分布的第i个网络结构单元的参数，是针对执行声速分布时结构相同且具有不同学习参数θ
2,i
的卷积神经网络，是第i个网络结构单元中声速分布的梯度，c
i-1,k
为第k个样本经过第i-1次迭代后的声速历史预测分布图，h(c
i-1,k
)是与c
i-1,k
有关的算子；δ
1,21
和δ
2,21
分别是卷积神经网络输出层的第一误差向量和第二误差向量；z
1,21
和z
2,21
分别是执行光吸收能量分布和执行声速分布时第21层网络产生的线性激活响应；第21层网络为卷积神经网络输出层；
“⊙”
是矩阵的hadamard积，表示矩阵对应元素点对点的乘法；σ'(
·
)为relu函数的一阶导数；δ
1,c
和δ
2,c
分别是执行光吸收能量分布和执行声速分布时第c层网络的第一误差向量和第二误差向量；c为特征图层序号，在第一递推公式、第二递推公式中c＝20,19,...,3,2，在求卷积神经网络每一层参数的梯度的公式中c＝2,3,...,21；w
1,c 1
和w
2,c 1
分别为执行光吸收能量分布和执行声速分布时第c 1层网络的权重参数矩阵；上标“t”表示求矩阵的转置；δ
1,c 1
和δ
2,c 1
分别是执行光吸收能量分布和执行声速分布时第c 1层网络的第一误差向量和第二误差向量；z
1,c
和z
2,c
分别为前向传播中执行光吸收能量分布和执行声速分布时第c层网络产生的线性激活响应；c(w,b)为总损失；λ1和λ2分别为第一权重和第二权重；w
1,c
和w
2,c
分别为执行光吸收能量分布和执行声速分布时第c层网络的权重参数矩阵；b
1,c
和b
2,c
分别为执行光吸收能量分布和执行声速分布时第c层网络的偏置参数向量；a
1,c-1
和a
2,c-1
分别为执行光吸收能量分布和执行声速分布时第c-1层网络的输出；wc和bc分别是第c层网络的权重参数矩阵和偏置参数向量。
[0122]
在采用基于模型的方法重建光声图像的过程中，采用传统优化方法迭代求解非线性最小二乘问题时往往非常耗时，且需要合理选择适当的正则化方法及其参数，重建精度很大程度上依赖于有关成像目标的先验假设模型。针对此问题，本发明提出一种卷积神经网络(convolutional neural network，cnn)模型，即联合迭代重建网络(joint iterative reconstruction network，jir-net)，将梯度下降信息拟合到cnn中，并利用反向传播的梯度下降求解非线性最小二乘问题，实现根据超声探测器在血管腔内采集的有限角度光声信
号，同时重建血管横截面上的光吸收能量分布图和声速分布图。
[0123]
本发明所称问题是以下述技术方案实现的：
[0124]
首先，采用计算机仿真的方法，建立适用于有限角度ivpat图像重建的深度学习数据集；然后，搭建jir-net模型，并预训练每一层网络，评估网络性能，进而优化各层网络的超参数；最后，将有限角度光声信号测量数据和预先设定的常数声速值共同输入训练后的jir-net中，输出光吸收能量和声速分布的估计值，并用灰度图像表示。
[0125]
如图1所示，本发明方法的具体实施步骤如下。
[0126]
一、构造用于有限角度ivpat图像重建的深度学习数据集
[0127]
步骤1-1建立血管横截面计算机仿真模型
[0128]
如图2所示，建立血管横截面计算机仿真模型，其中成像导管位于血管横截面中心，周围由内向外依次是血管腔、粥样硬化斑块和血管壁(三层结构，包括内膜、中膜和外膜)。
[0129]
步骤1-2建立ivpat的前向仿真模型
[0130]
采用蒙特卡罗模拟的方法，得到在短脉冲激光照射下的血管横截面上的光辐射通量，进而根据预设的组织光吸收系数得到光吸收能量密度：
[0131]
a(r)＝μa(r)φ(r)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0132]
式中，r是成像平面中的位置，a(r)是位置r处的光吸收能量密度，μa(r)是位置r处的光吸收系数，φ(r)是位置r处的光辐射通量。
[0133]
组织中的光吸收体基于光声效应产生的初始声压为：
[0134]
p0(r)＝γa(r)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0135]
式中，p0(r)是位置r处的初始声压，γ是组织的热膨胀系数。光声信号在均匀组织中的传播用如下波动方程描述：
[0136][0137]
式中，p(r,t)是位置r在时刻t的声压，i(t)是入射激光脉冲能量函数，c(r)是位置r处的声速，βe和c
p
分别是组织的等压膨胀系数和一定压力下的比热容。上述过程可用算子表示为：
[0138]
p＝h(c)a
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0139]
其中，p是到达探测器表面的声压向量，a是光吸收能量密度向量，h(c)是与声速c有关的算子。
[0140]
步骤1-3构建数据集
[0141]
采用上述方法分别仿真出全角度(360
°
)稀疏采样(采样率依次是0.1、0.3、0.5、0.7和0.9)光声信号以及有限角度(扫描角度依次是90
°
、180
°
和270
°
)稀疏采样(采样率依次是0.1、0.3、0.5、0.7和0.9)和密集采样(采样率等于1)光声信号数据，将其和预设的常数声速值共同作为样本输入，将前向仿真得到的光吸收能量分布和预设的声速分布作为样本的目标输出，得到ivpat图像重建所需的深度学习数据集。
[0142]
步骤1-4划分数据集
[0143]
将数据集中的所有样本随机打乱，按照6:2:2的比例将数据集划分为训练集、验证集和测试集，利用验证集选择最佳网络模型。把训练集中的样本平均分成n组，每组是一个
小批量训练集，共包含m个样本。
[0144]
二、设计前向传播的jir-net模型
[0145]
本发明方法搭建的jir-net结构如图3所示，搭建流程如图4所示，具体步骤如下：
[0146]
步骤2-1确定网络架构和参数。如附图3所示，jir-net中卷积操作的滤波核尺寸是3
×3×
s，s是当前层输入的特征图数量，采用步长为1的same卷积，网络单元的四个输入的滤波核初始数量(即特征通道数channels)均为16。激活函数采用线性整流函数relu。通过远跳连接将网络单元的输入与输出图像求和，得到当前网络单元的输出图像。
[0147]
步骤2-2初始化网络参数：设置最大网络层数i
max
＝4，网络的层序号i＝1，特征图层序号c＝1。
[0148]
步骤2-3初始化网络参数，令卷积次数conv＝0，特征通道数channels＝16。
[0149]
步骤2-4计算第i层网络结构单元中光吸收能量密度的梯度：
[0150][0151]
其中，(简写为)是第i层网络单元中光吸收能量密度的梯度，pk是训练集中第k个样本的光声信号测量值，a
i-1,k
是训练集中第k个样本经过i-1次迭代后的光吸收能量密度估计值，c
i-1,k
是训练集中第k个样本经过i-1次迭代后的声速估计值，h(c
i-1,k
)是与声速c
i-1,k
有关的算子，上标“t”表示求矩阵的转置。
[0152]
步骤2-5计算第i层网络结构单元中声速的梯度：
[0153]
伴随光声波动方程定义如下：
[0154][0155]
式中，q(r,t)是时刻t时位置r处的光声波的伴随波，t是采集光声信号的时间长度，s(r,t)是伴随源：
[0156][0157]
式中，r
l
是探测器的第l个测量位置，l是探测器测量位置的总数，是时刻t时探测器在位置r
l
处测得的光声信号测量值，p(r
l
,t)是时刻t时位置r
l
处的光声信号理论值，由式(3)求得，δ(r-r
l
)是狄拉克函数。式(6)的初始条件是：
[0158][0159]
根据式(3)和式(6)得到声速c(r)的梯度：
[0160][0161]
根据式(9)计算成像平面中各处的声速梯度，得到第i层网络单元中声速分布的梯度简写为
[0162]
步骤2-6将第1组小批量训练集中的m个样本输入jir-net，经c-1次卷积后产生特征图层。
[0163]
步骤2-7a
i-1,k
和经过conv次卷积之后的特征图层分别是map
c1
和map
c2
，c
i-1,k
和经过conv次卷积之后的特征图层分别是map
c3
和map
c4
。分别对map
c1
、map
c2
、map
c3
和map
c4
通过多个特征通道并行地执行卷积操作，产生多组线性激活响应，然后计算线性激活响应的relu函数值，得到卷积操作的特征图层。
[0164]
步骤2-8令c
←
c 1，conv
←
conv 1，channels
←
channels
×
2。
[0165]
步骤2-9执行步骤2-7，然后执行步骤2-10。
[0166]
步骤2-10令c
←
c 1，conv
←
conv 1，map
c,a
＝map
c1
map
c2
，map
c,c
＝map
c3
map
c4
。其中，map
c,a
和map
c,c
分别是a
i-1,k
和c
i-1,k
进行梯度拟合后的特征图层。
[0167]
步骤2-11将map
c,a
和map
c,c
进行拼接，得到特征图层mapc，且channels
←
channels
×
2。
[0168]
步骤2-12特征通道的数量不变，即令channels
←
channels。
[0169]
步骤2-13对mapc通过多个特征通道并行地执行分组卷积操作(group＝2)，产生多组线性激活响应。然后计算线性激活响应的relu函数值，得到分组卷积操作的特征图层。
[0170]
步骤2-14令c
←
c 1，conv
←
conv 1，channels
←
channels/2。
[0171]
步骤2-15执行步骤2-13，然后执行步骤2-16。
[0172]
步骤2-16令c
←
c 1，conv
←
conv 1，沿通道方向拆分mapc为map
c,a
和map
c,c
，且map
c,a
和map
c,c
的特征通道数量均为channels/2。
[0173]
步骤2-17令channels＝1，分别对map
c,a
和map
c,c
通过多个特征通道并行地执行卷积操作，产生多组线性激活响应。然后计算线性激活响应的relu函数值，得到卷积操作的特征图层。
[0174]
步骤2-18令c
←
c 1，conv
←
conv 1，map
c,a
＝map
c,a
a
i-1,k
，map
c,c
＝map
c,c
c
i-1,k
，之后计算线性激活响应的relu函数值，得到卷积操作的特征图层。
[0175]
步骤2-19令i
←
i 1。如果i》4，执行步骤2-20；否则返回步骤2-3。
[0176]
步骤2-20输出m个样本的预测图像。此时，c＝21。
[0177]
三、优化网络参数
[0178]
jir-net的训练流程如图5所示，具体步骤如下：
[0179]
步骤3-1初始化网络的层序号i，令i＝1。
[0180]
步骤3-2初始化遍历训练集中所有小批量训练集的次数epoch＝0。
[0181]
步骤3-3初始化小批量训练集的索引i＝1。
[0182]
步骤3-4设计损失函数
[0183]
将估算光吸收能量分布和声速分布的任务分别命名为任务1和任务2。
[0184]
任务1的损失函数为：
[0185][0186]
其中：
[0187]
pk＝h(c
0,k
)a
0,k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0188]
式中，w是网络的权重参数集，b是网络的偏置参数集，k是第i组小批量训练集中样本的序号(k＝1,2,3,
…
,m)，a
true,k
是第k个样本的光吸收能量密度理论值(由前向仿真得到)，a
i,k
是第k个样本经过第i次迭代后的光吸收能量密度估计值，θ
1,i
是任务1第i个网络结构单元的参数，是针对任务1、结构相同且具有不同学习参数θ
1,i
的卷积神经网络，c
0,k
是第1次迭代时输入的声速初始值(即预设的常数声速)，a
0,k
是第1次迭代时输入的光吸收能量密度初始值。
[0189]
任务2的损失函数为：
[0190][0191]
其中：
[0192][0193]
式中，c
true,k
是第k个样本的声速理论值，c
i,k
是第k个样本经过第i次迭代后的声速估计值，θ
2,i
是任务2第i个网络结构单元的参数，是针对任务2结构相同且具有不同学习参数θ
2,i
的卷积神经网络。
[0194]
网络的总损失函数为：
[0195]
c(w,b)＝λ1c1(w,b) λ2c1(w,b)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0196]
其中，λ1和λ2是权重参数。
[0197]
步骤3-5将第i组小批量训练集中的m个样本输入搭建好的jir-net中进行前向传播。
[0198]
步骤3-6根据递推关系，计算网络输出层各参数的梯度。
[0199]
任务1输出层各参数的梯度为：
[0200][0201]
任务2输出层各参数的梯度为：
[0202][0203]
其中，z
1,21
和z
2,21
分别是任务1和任务2第21层网络产生的线性激活响应，w
1,21
和w
2,21
分别是任务1和任务2第21层网络的权重参数矩阵，b
1,21
和b
2,21
分别是任务1和任务2第21层网络的偏置参数向量，a
1,20
和a
2,20
分别是任务1和任务2第20层网络的输出，
“⊙”
是矩阵的hadamard积，表示矩阵对应元素点对点的乘法，σ'(
·
)是relu函数的一阶导数。
[0204]
步骤3-7分别计算任务1和任务2网络各层的误差向量。
[0205]
任务1网络输出层(即第21层)的误差向量为：
[0206][0207]
任务2网络输出层(即第21层)的误差向量为：
[0208][0209]
其中，δ
1,21
和δ
2,21
分别是任务1和任务2网络输出层的误差向量。
[0210]
根据相邻层之间的误差递推关系分别求得任务1和任务2网络其余层的误差向量：
[0211]
δ
1,c
＝(w
1,c 1
)
t
δ
1,c 1
⊙
σ'(z
1,c
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0212]
δ
2,c
＝(w
2,c 1
)
t
δ
2,c 1
⊙
σ'(z
2,c
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(21)
[0213]
式中，c＝20,19,...,3,2，δ
1,c
和δ
2,c
分别是任务1和任务2第c层网络的误差向量，w
1,c 1
和w
2,c 1
分别是任务1和任务2第c 1层网络的权重参数矩阵，z
1,c
和z
2,c
分别是前向传播中任务1和任务2第c层网络产生的线性激活响应。
[0214]
步骤3-8计算网络各层参数的梯度。
[0215]
任务1各层网络参数的梯度为：
[0216][0217]
任务2各层网络参数的梯度为：
[0218][0219]
其中，c＝2,3,
…
,21，c1(w,b)和c2(w,b)分别是任务1和任务2网络的损失函数，w
1,c
和w
2,c
分别是任务1和任务2第c层网络的权重参数矩阵，b
1,c
和b
2,c
分别是任务1和任务2第c层网络的偏置参数向量，a
1,c-1
和a
2,c-1
分别是前向传播中任务1和任务2第c-1层网络的输出。
[0220]
总网络各层参数的梯度为：
[0221][0222]
其中，c＝2,3,
…
,21，wc和bc分别是第c层网络的权重参数矩阵和偏置参数向量。
[0223]
步骤3-9利用adam算法调整网络各层参数的尺度，利用梯度正则化(gradient normalization,gradnorm)算法调整参数λ1和λ2，使其达到最优。
[0224]
步骤3-10令i
←
i 1。若i《n 1，则转向步骤3-4；否则令epoch
←
epoch 1，并执行步骤3-11。
[0225]
步骤3-11判断是否停止遍历：若epoch《a，a是epoch的最大值，则转向步骤3-3；否
则停止遍历，执行步骤3-12。
[0226]
步骤3-12令i
←
i 1，若i《5，则转向步骤3-2，否则执行步骤3-13。
[0227]
步骤3-13确定网络的各个参数，得到最终的网络模型。
[0228]
四、利用优化后的jir-net重建光吸收能量分布图和声速分布图
[0229]
将超声探测器实际采集的有限角度光声信号测量数据和预先设置的常数声速值共同输入训练好的jir-net中，输出即为高质量的光吸收能量分布图和声速分布图。
[0230]
此外，本发明还提供了一种光吸收能量分布图和声速分布图的重建系统，包括：
[0231]
测量成像数据获取模块，用于获取待成像平面的测量成像数据；测量成像数据包括有限角度光声信号测量数据和常数声速值；
[0232]
分布图重建模块，用于将测量成像数据输入到图像重建模型中，得到待成像平面处的光吸收能量重建分布图和声速重建分布图；图像重建模型是根据历史成像数据，利用梯度下降算法对卷积神经网络进行训练得到的。
[0233]
训练集获取模块，用于获取训练集；训练集包括待测血管多个平面的历史成像数据、光吸收能量历史分布图和声速历史分布图；
[0234]
图像重建模型确定模块，用于以历史成像数据为输入，以光吸收能量历史分布图和声速历史分布图为期望输出，利用梯度下降算法对卷积神经网络进行训练，得到图像重建模型。
[0235]
具体的，图像重建模型确定模块，具体包括：
[0236]
训练子划分单元，用于将训练集划分为多个训练子集；
[0237]
初始卷积神经网络确定单元，用于确定训练前的卷积神经网络为第i
×
epoch
×
i阶卷积神经网络；
[0238]
网络层序号赋值单元，用于令网络层序号i＝1；
[0239]
训练集遍历次数赋值单元，用于令训练集遍历次数epoch＝0；
[0240]
当前训练子集确定单元，用于确定任一训练子集为当前训练子集；
[0241]
索引赋值单元，用于令当前训练子集的索引i＝1；
[0242]
训练单元，用于将当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
(i 1)阶卷积神经网络；
[0243]
当前训练子集更新单元，用于更新当前训练子集，令i的数值增加1，并调用训练单元直至历遍训练子集，得到第i
×
(epoch 1)
×
1阶卷积神经网络；
[0244]
训练集遍历次数更新单元，用于令epoch的数值增加1并调用索引赋值单元直至epoch大于训练集遍历次数阈值，得到第(i 1)
×0×
1阶卷积神经网络；
[0245]
网络层序号更新单元，用于令i的数值增加1并调用训练集遍历次数赋值单元直至i大于网络层序号阈值，得到i
×
epoch
×
i阶卷积神经网络作为图像重建模型。
[0246]
进一步的，训练单元，具体包括：
[0247]
历史预测分布图确定子单元，用于将当前训练子集输入到第i
×
epoch
×
i阶卷积神经网络中，得到第i
×
epoch
×
i次迭代时的光吸收能量历史预测分布图和第i
×
epoch
×
i次迭代时的声速历史预测分布图；
[0248]
第一损失确定子单元，用于根据光吸收能量历史分布图和第i
×
epoch
×
i次迭代时的光吸收能量历史预测分布图，确定第i
×
epoch
×
i阶卷积神经网络执行光吸收能量分
布时的损失为第一损失；
[0249]
第二损失确定子单元，用于根据声速历史分布图和第i
×
epoch
×
i次迭代时的声速历史预测分布图，确定第i
×
epoch
×
i阶卷积神经网络执行声速分布时的损失为第二损失；
[0250]
第一误差向量确定子单元，用于根据第一损失确定第i
×
epoch
×
i阶卷积神经网络输出层的第一误差向量；
[0251]
第二误差向量确定子单元，用于根据第二损失确定第i
×
epoch
×
i阶卷积神经网络输出层的第二误差向量；
[0252]
第一误差向量组确定子单元，用于根据第一误差向量和第一递推公式确定第i
×
epoch
×
i阶卷积神经网络每一层的第一误差向量；
[0253]
第二误差向量组确定子单元，用于根据第二误差向量和第二递推公式确定第i
×
epoch
×
i阶卷积神经网络每一层的第二误差向量；
[0254]
梯度确定子单元，用于根据第一权重、第二权重、多个第一误差向量和多个第二误差向量，确定第i
×
epoch
×
i阶卷积神经网络每一层参数的梯度；第一权重为卷积神经网络执行光吸收能量分布时的权重；第二权重为卷积神经网络执行声速分布时的权重；
[0255]
参数尺度调整子单元，用于基于第i
×
epoch
×
i阶卷积神经网络每一层参数的梯度，利用adam算法调整第i
×
epoch
×
i阶卷积神经网络中每一层参数的尺度，得到第i
×
epoch
×
(i 1)阶卷积神经网络。
[0256]
其中，参数尺度调整子单元还用于在所述利用adam算法调整第i
×
epoch
×
i阶卷积神经网络中每一层参数的尺度之后，利用梯度正则化算法调整第一权重和第二权重。
[0257]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0258]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。