AAR-Net: a deep neural network for photoacoustic image reconstruction in heterogeneous acoustic media
-
摘要: 在光声成像中，由于组织的吸收和扩散等引起的超声波衰减、由声速变化引起的相位偏差以及与声衰减相关的信号波形展宽都会降低图像的空间分辨率，针对该问题，提出一种基于深度学习的声学特性非均匀组织图像重建方法。通过将深度梯度下降（deep gradient descent，DGD）网络与U-Net相结合构建声伪影去除网络（acoustic artifacts removal network，AAR-Net）。DGD模块利用梯度信息减少非均匀声学特性对重建图像质量的影响，实现信号域到图像域的转换。U-Net模块实现对DGD模块输出的低质量图像的优化，实现图像域到图像域的转换。仿真、仿体和在体试验结果表明，与传统的非学习图像重建方法和最新的基于图像后处理的深度学习方法相比，采用该方法重建的图像结构相似度和峰值信噪比分别可提高约20%和10%。AAR-Net无需任何有关成像对象声学特性的先验知识，即可重建高质量图像。Abstract: Photoacoustic imaging suffers from degraded image quality owing to distorted and attenuated ultrasound waves propagating in an acoustic attenuating medium, phase deviation caused by changes in sound speed, and signal broadening related to acoustic attenuation. To address this issue, a deep learning method is proposed to reconstruct photoacoustic images of acoustically heterogeneous medium. A deep neural network is constructed, named acoustic artifacts removal network (AAR-Net) by combining deep gradient descent (DGD) network with U-Net. The DGD module aims to achieve the conversion from signal domain to image domain, which uses the gradient information to reduce the impact of heterogeneous acoustic properties on reconstructed image quality. U-Net module aims to optimize low-quality images output by the DGD module and realize the image-to-image conversion. The simulation, phantom and in vivo studies show that the proposed method outperforms traditional non-learning methods and the state-of-the-art post-processing based deep learning method. The image similarity and peak signal-to-noise ratio obtained by this method are improved by about 20% and 10%, respectively. AAR-Net enables reconstruction of high-quality images without any prior knowledge of acoustic properties of imaging objects.
-
光声成像是一种新型功能成像方法，由于辐射源是短脉冲激光，因此可以产生类似于其他光学成像的高对比度图像。同时，基于光声效应产生的超声波（即光声波）在生物组织中传播时具有低散射特性，因此还具有超声成像高穿透深度的优势[1]。
光声图像重建是根据超声换能器采集的光声信号，通过求解声学逆问题重建光吸收能量密度（absorbed optical energy density，AOED）空间分布图或初始声压分布图。传统方法包括滤波反投影（filtered back projection，FBP）[2]、时间反演（time reversal，TR）[3]和延迟求和（delay and sum，DAS）[4]等。为了简化问题，上述方法通常假设成像目标是具有均匀声学特性的无损介质，超声波以恒定的速度在介质中传播。然而在实际应用中，人体组织的声学特性表现出空间不均匀性。不同成分组织之间的声速差异会影响光声信号到达探测器的时间，进而影响根据传播时间计算的声源到探测器之间的距离，此时使用恒定声速重建的图像中就会出现目标畸变和错位等现象。此外，当光声波在随机声学异质介质中传播时，特别是在声学特性严重失配的组织边界处会发生吸收、扩散、反射和散射等多种声学效应，导致部分声波偏离其原始传播路径，探测器采集到衰减的声压信号。而且这种衰减是频率相关的，高频信号分量比低频信号分量的衰减更严重，最终导致成像深度降低，信号的波形展宽和带宽降低也限制了成像的空间分辨率。
光声图像质量与波束成形过程密切相关，波束成形技术基于对超声波传播时间的测量，然而这种测量没有考虑超声波在介质中传播时可能被多次反射，且每次的反射波都可能被超声换能器检测到的情况，进而在重建图像的过程中将混响信号映射到错误的位置，导致图像中出现反射伪影。反射伪影通常出现在比实际光吸收体更深的位置，若出现在感兴趣区域内，则可能被误判成有用目标，降低图像识别的精度[5]。此外，反射伪影可能比组织产生的原始光声信号强度更高，导致成像深度下降。由于反射杂波与组织的光学特性有关，因此不能通过简单的信号平均进行抑制。目前抑制反射伪影的方法主要有5种：1）杂波去相关理论[6]，仅适用于易产生变形的组织，而且需要额外的机械装置以及对操作者进行专门培训。2）短滞后空间相干度（short-lag spatial coherence，SLSC）法[7]，根据光声信号和反射杂波的SLSC之间的差异区分非相干杂波和有用信号，进而将SLSC作为重建图像时不同信号成分的加权因子，实现图像的无杂波重建。该方法无法区分空间相干性类似的有用信号和反射杂波，而且SLSC的计算依赖于对与波束成形相关的几何路径的时间测量，在存在高回声反射或强反射的情况下无法保证成像精度。3）多波长激励[8]，仅适用于光声光谱成像，当在不同深度处产生类似的光谱响应时，或者反射伪影与另一个光吸收体（或另一个反射伪影）出现在相同位置时，无法获得满意的成像结果。4）光声−超声双模态成像方法，利用聚焦超声波[9]或平面超声波[10]模拟向反射体传播的部分光声波，获得只包含反射波的超声图像，进而得到去除反射伪影的光声图像。此类方法需要采集和处理大量超声信号，并对光声和超声测量数据进行精准匹配，过程非常耗时，而且超声和光声信号具有相同接收路径的假设并非在任何情况下都成立。5）基于模型的图像重建[11]，通过求解前向成像模型的逆问题实现图像重建，建立成像模型时，通过充分考虑复杂的成像条件和成像目标特性，可以避免产生伪影的因素。然而对光声信号的产生以及在介质中的传播过程进行准确的建模是非常困难的，而且此类方法需要采用优化算法迭代求解前向模型的逆问题，计算成本很高。
除了反射以外，随机散射也是导致光声波衰减的主要原因之一。为了补偿声散射，Deán-Ben等[12]在成像目标和探测器之间放置了强吸收体，并根据探测器采集的光声信号估计成像目标中声散射体的空间分布。Yin等[13]在光声扫描之前利用超声扫描测量散射介质的声学特性，进而重构每个探测器位置的格林函数。前述方法需要额外的强吸收体或超声扫描，提高了成像复杂度，不适用于一些特殊的成像场景（如内窥式成像）。Rui等[14]设计了一种相关滤波器，用于分离声压信号中的直接波和散射波，但由于计算复杂度较高，其在实际场景中应用的可行性仍有待验证。
近年来，随着高性能并行计算技术和大数据技术的飞速发展，深度学习在光声成像领域的应用日益广泛，主要包括图像重建和图像后处理两方面[15-16]。前者是利用深度神经网络（deep neural network，DNN）实现从“信号”到“图像”的映射，将光声信号测量数据直接输入网络，通过训练使网络学习目标图像的先验知识，该过程非常复杂，且耗时较长，但可以重构出大部分图像特征[17-19]。后者是首先采用常规方法（如FBP、TR和DAS等）根据光声信号重建出低质量的初始图像，再将其输入训练后的神经网络中进行优化，最终输出高质量图像，实现“图像”到“图像”的映射[20-22]。该过程简单快速，但处理效果受到初始重建质量的限制，一些在初始重建过程中丢失的信息很难通过“学习”恢复。
本研究提出一种用于声学异质介质光声图像重建的方法，搭建声伪影去除网络（acoustic artifacts removal network，AAR-Net），不需要任何有关声衰减介质的先验知识，也无需修改成像装置，即可实现从探测器采集的声压信号到高质量图像之间的“端到端”映射。通过仿真、仿体和在体研究证明，与假定超声波无损传输的图像重建方法相比，本研究方法可以有效提高成像质量。
1. 声学异质介质的光声图像重建
如图1所示，AAR-Net的输入是探测器在成像平面内采集的声压信号矩阵，输出高质量的AOED分布图。
图 1 本研究方法流程Fig. 1 Flow chart of the proposed method
下载:
全尺寸图片
1.1 成像模型
为了简单起见，本研究对成像场景进行了如下假设：光源位置固定，成像平面内的光照均匀，不考虑非固定照明和非均匀光通量的影响。在二维照明和圆形扫描几何的情况下，成像平面内的图像重建可以看作是一个二维过程。采用超声换能器阵列采集光声信号，其中每个阵元被理想化为点探测器，忽略其有限孔径和有限带宽的影响，不考虑有限角度扫描和稀疏采样的问题，探测器可以采集到成像平面内的完备信号。
采用蒙特卡罗（Monte Carlo，MC）方法模拟入射光在组织中的传输过程，得到成像平面中各处的光辐射通量[23]，进而得到AOED：
$$ A({\boldsymbol{r}}) = {\mu _a}({\boldsymbol{r}})\varPhi ({\boldsymbol{r}}) $$ (1) 式中，
$A({\boldsymbol{r}})$ 、${\mu _a}({\boldsymbol{r}})$ 和$\varPhi ({\boldsymbol{r}})$ 分别是成像平面中位置r处的AOED、光吸收系数和光辐射通量。光吸收体产生的初始声压与AOED成正比例关系
$$ {p_0}({\boldsymbol{r}}) = \varGamma A({\boldsymbol{r}}) $$ (2) 式中：
${p_0}({\boldsymbol{r}})$ 是位置${\boldsymbol{r}}$ 处的初始声压；$\varGamma$ 是热膨胀系数（本研究将其设定为1）。采用如下波动方程描述光声波在声学特性非均匀介质中的传播[24]：
$$ \left\{ \begin{gathered} {{\partial {\boldsymbol{u}}({\boldsymbol{r}},t)} \mathord{\left/ {\vphantom {{\partial {\boldsymbol{u}}({\boldsymbol{r}},t)} {\partial t}}} \right. } {\partial t}} = - \frac{1}{{{\rho _0}}}\nabla p({\boldsymbol{r}},t) \\ {{\partial \rho ({\boldsymbol{r}},t)} \mathord{\left/ {\vphantom {{\partial \rho ({\boldsymbol{r}},t)} {\partial t}}} \right. } {\partial t}} = - {\rho _0}\nabla \cdot {\boldsymbol{u}}({\boldsymbol{r}},t) \\ p({\boldsymbol{r}},t) = {[c({\boldsymbol{r}})]^2}\rho ({\boldsymbol{r}},t)\left[1 + \tau ({\boldsymbol{r}})\frac{\partial }{{\partial t}}{( - {\nabla ^2})^{\frac{a}{2} - 1}} + \eta ({\boldsymbol{r}}){( - {\nabla ^2})^{\frac{{a + 1}}{2} - 1}}\right] \\ \end{gathered} \right. $$ (3) 式中：
$p({\boldsymbol{r}},t)$ 和${\boldsymbol{u}}({\boldsymbol{r}},t)$ 分别是位置${\boldsymbol{r}}$ 处、时刻t的声压和质点振动速度；$c({\boldsymbol{r}})$ 是声速；$\,\rho ({\boldsymbol{r}},t)$ 是非均匀介质的声学密度；$ {\,\rho _0} $ 是均匀介质密度；$ \nabla $ 是哈密顿（Hamiltonian）算子；$ a $ 是幂律指数（对于生物组织其取值范围是1≤a≤1.5）；$\tau ({\boldsymbol{r}})$ 和$\eta ({\boldsymbol{r}})$ 分别是声吸收系数和色散系数[24]$$ \left\{ {\begin{array}{*{20}{l}} {\tau ({\boldsymbol{r}}) = - 2\mu {{[c({\boldsymbol{r}})]}^{a - 1}}} \\ {\eta ({\boldsymbol{r}}) = 2\mu {{[c({\boldsymbol{r}})]}^a}\tan (\text{π} a/2)} \end{array}} \right. $$ (4) 式中：
$\mu = ({10^{ - 7}}/2\text{π} ){\rm{c}}{{\rm{m}}^{ - 1}}{\rm{ra}}{{\rm{d}}^{ - 1}}{\rm{s}}$ 是幂律吸收因子[25-26]。式(3)的初始条件为$$ \left\{ {\begin{array}{*{20}{l}} {p({\boldsymbol{r}},t){|_{t = 0}} = {p_0}({\boldsymbol{r}})} \\ {u({\boldsymbol{r}},t){|_{t = 0}} = 0} \end{array}} \right. $$ (5) 采用k空间伪谱法[27]对式(3)离散化，得到
$$ \left\{ {\begin{array}{*{20}{l}} {{{\partial p({\boldsymbol{r}},t)} \mathord{\left/ {\vphantom {{\partial p({\boldsymbol{r}},t)} {\partial \xi }}} \right. } {\partial \xi }} = {\mathcal{F}^{ - 1}}\{ {\text{j}}{k_\xi }\kappa \mathcal{F}[p({\boldsymbol{r}},t)]\} } \\ {{u_\xi }({\boldsymbol{r}},t + \Delta t) = {u_\xi }({\boldsymbol{r}},t) - \dfrac{{\Delta t}}{{{\rho _0}}}\dfrac{{\partial p({\boldsymbol{r}},t)}}{{\partial \xi }}} \\ {{{\partial {u_\xi }({\boldsymbol{r}},t + \Delta t)} \mathord{\left/ {\vphantom {{\partial {u_\xi }({\boldsymbol{r}},t + \Delta t)} {\partial \xi }}} \right. } {\partial \xi }} = {\mathcal{F}^{ - 1}}\{ {\text{j}}{k_\xi }\kappa \mathcal{F}[{u_\xi }({\boldsymbol{r}},t + \Delta t)]\} } \\ {{\rho _\xi }({\boldsymbol{r}},t + \Delta t) = {\rho _\xi }({\boldsymbol{r}},t) - \Delta t{\rho _0}{{\partial {u_\xi }({\boldsymbol{r}},t + \Delta t)} \mathord{\left/ {\vphantom {{\partial {u_\xi }({\boldsymbol{r}},t + \Delta t)} {\partial \xi }}} \right. } {\partial \xi }}} \\ {p({\boldsymbol{r}},t + \Delta t) = {{[c({\boldsymbol{r}})]}^2}\left[\displaystyle\sum\limits_\xi {{\rho _\xi }} ({\boldsymbol{r}},t + \Delta t) + {{a_{{\rm{bs}}}}} + {{d_{{\rm{isp}}}}}\right]} \end{array}} \right. $$ (6) 其中
$$ \kappa {\text{ = }}{\rm{sinc}}\left[\frac{1}{2}c({\boldsymbol{r}}){k_\xi }\Delta t\right] $$ (7) abs是吸收项
$$ {{a_{{\rm{bs}}}}} = \tau ({\boldsymbol{r}}){\mathcal{F}^{ - 1}}\left\{ {\left({k_\xi }\kappa \right)^{a - 2}}\mathcal{F}\left[\sum\limits_\xi {{\rho _\xi }\left({\boldsymbol{r}},t + \Delta t\right)} \right]\right\} $$ (8) disp是色散项
$$ {{d_{{\rm{isp}}}}} = \eta ({\boldsymbol{r}}){\mathcal{F}^{ - 1}}\left\{ {\left({k_\xi }\kappa \right)^{a - 1}}\mathcal{F}\left[\sum\limits_\xi {{\rho _\xi }\left({\boldsymbol{r}},t + \Delta t\right)} \right]\right\} $$ (9) 式中：j是虚数单位；
$ \xi = (x,y) $ ；$ {k_\xi } $ 是在ξ方向上的空间波数分量；$ \mathcal{F}[ \cdot ] $ 和$ {\mathcal{F}^{ - 1}}[ \cdot ] $ 分别是二维傅里叶变换和其逆变换；$ {\rho _\xi }({\boldsymbol{r}},t) $ 和$ {u_\xi }({\boldsymbol{r}},t) $ 分别是介质声学密度和振动速度在$ \xi $ 方向上的分量；$\Delta t = {{C_{{\rm{FL}}}}} \cdot \Delta x/{c_0}$ 是时间步长；$ {c_0} $ 是均匀介质中的声速；$ \Delta x $ 是网格尺寸；CFL=0.3是仿真精确度和计算速度之间的折中系数。初始声学密度为$$ {\displaystyle \sum _{\xi }{\rho }_{\xi }({\boldsymbol{r}},t)}|{}_{t=0}=\frac{1}{{[c({\boldsymbol{r}})]}^{2}}{p}_{0}({\boldsymbol{r}},t)|{}_{t=0} $$ (10) 由AOED到产生光声信号的前向过程，用算子
$ H( \cdot ) $ 表示为$$ {\boldsymbol{p}} = H({\boldsymbol{A}}) $$ (11) 式中：
${\boldsymbol{A}}$ 是成像平面内的AOED矩阵；${\boldsymbol{p}}$ 是到达探测器表面的声压矩阵。图像重建就是下式的逆过程，用算子表示为$$ \hat {\boldsymbol{A}} = {H^*}({{\boldsymbol{p}}_{\text{m}}}) $$ (12) 式中：
$ {H^*}( \cdot ) $ 是$ H( \cdot ) $ 的伴随算子，$\hat {\boldsymbol{A}}$ 是重建的AOED分布图，${{\boldsymbol{p}}_{\text{m}}}$ 是探测器测量的声压矩阵。1.2 AAR-Net的结构
如图2所示，AAR-Net由深度梯度下降（deep gradient descent，DGD）模块[28]和U-Net模块构成。DGD模块用于根据光声信号测量数据重建初始图像，实现从信号域到图像域的转换。U-Net模块用于对DGD模块输出的初始图像进行优化，得到高质量图像，实现从图像域到图像域的转换。
图 2 AAR-Net体系结构Fig. 2 Architecture of AAR-Net下载:
全尺寸图片
DGD模块包含9个卷积层，分为编码层和解码层，编码层的作用是通过卷积在提取特征的同时增加层间通道数量，实现对前一个图层特征的结构细化；解码层的作用是通过减少层间滤波核的数量合并图层的相似特征。卷积操作采用的是滤波核尺寸为3×3×T、步长为1的same卷积，T是当前层输入的特征图数量，2个输入图像的滤波核初始数量即特征通道数是16，激活函数采用线性整流函数ReLU。通过远跳连接将初始输入与输出图像进行像素值求和，得到当前层的输出图像。
DGD结构单元的操作表示为
$$ {{\boldsymbol{A}}_1} = {N_\theta }({{\boldsymbol{A}}_0},{\boldsymbol{G}}) $$ (13) 式中：
${{\boldsymbol{A}}_1}$ 是DGD模块输出的AOED分布图；$ {N_\theta }( \cdot ) $ 表示学习参数为$ \theta $ 的DGD网络。DGD网络的2个输入分别是初始AOED分布图和光声信号测量值与理论值之间的拟合梯度${\boldsymbol{G}}$ $$ {{\boldsymbol{A}}_0} = {H^*}({{\boldsymbol{P}}_{\rm{m}}}) $$ (14) $$ {\boldsymbol{G}} = \nabla d({{\boldsymbol{P}}_{\rm{m}}},H({{\boldsymbol{A}}_0})) $$ (15) 式中：
${{\boldsymbol{P}}_{\text{m}}}$ 是光声信号测量值矩阵；$ \nabla d( \cdot ) $ 是数据拟合梯度运算；$H({{\boldsymbol{A}}_0})$ 是前向成像算子输出的光声信号理论值。U-Net模块共包含26个卷积层，其中卷积操作采用的是滤波核尺寸为3×3×T、步长为1的same卷积，最后一次卷积操作采用滤波核尺寸为1×1×T的1×1卷积。滤波核初始数量即特征通道数是64。激活函数采用ReLU，池化操作采用滤波核尺寸为2×2×T、步长为2的最大池化，上采样操作采用的滤波核尺寸为2×2×T，步长为2。
1.3 训练网络
训练AAR-Net时，采用误差反向传播算法计算网络中每一层的梯度，利用自适应矩估计（adaptive moment estimation，Adam）梯度下降算法[29]沿各层参数梯度方向对网络参数进行优化。
AAR-Net的总损失函数定义为
$$ \mathcal{L} = {\mathcal{L}_{{\rm{DGD}}}} + {\mathcal{L}_{\rm{U}}} $$ (16) 式中：
$ {\mathcal{L}_{{\rm{DGD}}}} $ 和$ {\mathcal{L}_{\rm{U}}} $ 分别是DGD模块和U-Net模块的损失函数$$ \left\{ {\begin{array}{*{20}{l}} {{\mathcal{L}_{{\rm{DGD}}}}{\text{ = }}\dfrac{1}{M}\displaystyle\sum\limits_{i = 1}^M {\left\| {{{\boldsymbol{A}}_{1,i}} - {{\boldsymbol{A}}_{{\rm{true}},i}}} \right\|_2^2} } \\ {{\mathcal{L}_{\rm{U}}} = \dfrac{1}{M}\displaystyle\sum\limits_{i = 1}^M {\left\| {{{\boldsymbol{A}}_i} - {{\boldsymbol{A}}_{{\rm{true}},i}}} \right\|_2^2} } \end{array}} \right. $$ (17) 式中：
$ M $ 是小批量训练集中的样本个数；${{\boldsymbol{A}}_{{\rm{true}},i}}$ 是训练集中第$ i $ 个样本的期望输出图像；${{\boldsymbol{A}}_i}$ 是U-Net模块输出的图像；${{\boldsymbol{A}}_{1,i}}$ 是第i个样本输入DGD模块后输出的图像：$$ {{\boldsymbol{A}}_{1,i}} = {N_\theta }({{\boldsymbol{A}}_{0,i}},\nabla d({{\boldsymbol{P}}_{\rm{m}}},H({{\boldsymbol{A}}_{0,i}}))) $$ (18) 式中：
${{\boldsymbol{A}}_{0,i}}$ 是第i个样本的输入图像；$ \nabla d( \cdot ) $ 是光声信号测量值和理论值的拟合梯度$$ \nabla d({{\boldsymbol{P}}_{\rm{m}}},H({{\boldsymbol{A}}_{0,i}})) = {H^ * }(H({{\boldsymbol{A}}_{0,i}}) - {{\boldsymbol{P}}_{\rm{m}}}) $$ (19) 在AAR-Net中，DGD模块实现“信号”到“图像”映射，U-Net模块实现“图像”到“图像”的映射，二者任务特征、输入数据和模型结构复杂性之间的差异导致其学习速度不平衡，U-Net模块的收敛速度比DGD模块快得多。因此在训练开始时就直接优化整个框架的效率并不高。为了加快整体网络的收敛速度，本研究先对DGD进行了预训练，待其达到稳定的性能之后，再通过最小化损失函数式(16)训练整个联合学习框架。
网络训练中设置学习率为0.0002，批处理大小为16，epoch最大值是200。损失函数曲线如图3所示，可见随着epoch的增加，损失函数值逐渐减小，当epoch>200时，变化趋于平缓。
图 3 网络的损失函数曲线Fig. 3 Loss of network training下载:
全尺寸图片
1.4 数据集
本研究分别建立了仿真、仿体和在体数据集，用于训练、验证和测试AAR-Net。
1.4.1 仿真数据集
建立包含不同组织成分的计算机数值仿体，并设置不同组织成分的特性参数，如表1所示，不同组织成分的声速和密度服从以表1中所列值为均值、方差分别为5和0.05的高斯分布。采用MATLAB MCXLAB工具包仿真得到成像平面中的光辐射通量，其中入射光波长是1.7 μm，脉冲宽度是20 ns，入射光子总数是106，进而根据预设的光吸收系数得到成像平面内AOED的空间分布。采用MATLAB K-wave工具包[30]求解式(6)得到由声学异质组织产生并最终到达探测器表面的声压信号，其中加入信噪比为30 dB的随机噪声。将格式为灰度图像的仿真光声信号分布图和AOED分布图分别作为样本的输入和期望输出，构造仿真数据集，部分样本示例如图4所示，成像平面大小是67 mm × 67 mm，图像尺寸是256像素×256像素。
表 1 数值仿体中不同组织成分的特性参数Table 1 Characteristic parameters of different tissue components in numerical phantoms组织名称 折射率 反射率 吸收系数/cm‒1 散射系数/cm‒1 各向异性因子 声速/(m·s‒1) 密度/(kg·L‒1) 径向厚度/mm 外壁 1.42 0.08 0.70 6 0.85 1600 1.026 0.00~0.50 内壁 1.42 0.05 0.40 6 0.85 1580 1.073 0.10~0.40 纤维帽 1.47 0.24 0.78 25 0.78 1610 1.058 0.20~0.50 钙化 1.47 0.34 0.83 560 0.78 1540 1.668 0.30~0.50 脂质池 1.47 0.42 0.90 520 0.82 1650 0.947 0.10~0.30 混合钙化 1.47 0.23 0.96 550 0.80 1650 0.947 0.01~0.30 内腔 1.35 0.12 0 610 0.99 1540 1.065 1.20~1.70 
图 4 仿真数据集中的样本图像示例Fig. 4 Examples of samples in simulation datasets下载:
全尺寸图片
原始数据集中共包含2000对样本，将所有样本随机打乱，按照6∶2∶2的比例划分为训练集、验证集和测试集，利用验证集选择最佳网络模型。为了避免产生过拟合，采用数据增强的方法对训练集进行扩充，包括对样本进行随机旋转（‒180°~180°）、随机平移（不超过图像宽度的10%）和翻转，将样本数量扩充至4800对。
1.4.2 仿体数据集
采用琼脂、明胶、蓖麻油和乳化剂等材料制作了14个直径约为20 mm、高度约为80 mm的柱状仿体，其中包含不同形状、不同位置、不同材质（如印度墨水、铅芯和铁丝等）的嵌入物，旨在模拟具有不同光吸收特性的吸收体，部分仿体照片如图5所示。
图 5 仿体的实物照片Fig. 5 Photos of phantoms下载:
全尺寸图片
用于采集仿体光声信号的成像系统采用Nd:YAG激光器输出的全波波长（680 ~900 nm）可调谐激光作为光源，脉冲重复率是10 Hz，脉冲宽度是7 ns，最大脉冲能量为120 mJ。超声探测装置由256个中心频率为5 MHz、带宽为60%的聚焦超声换能器阵元排列成覆盖角度为270°的环形阵列。对14个仿体进行扫描共获得4379个切片的扫描数据，随机打乱后，按照6∶3∶1的比例划分为训练集、验证集和测试集。采用数据增强技术对训练集进行扩增，将样本总数扩充至7610对。
1.4.3 在体数据集
在体数据集由对活体小鼠进行全身扫描采集的光声信号构成。成像系统采用Nd:YAG激光器输出的715、730、760、800、850 nm可调谐激光作为光源，在730 nm处的最大入射脉冲能量约为70 mJ，脉冲重复率和宽度分别是10 Hz和8 ns。超声探测装置由128个中心频率为5 MHz、带宽为60%的聚焦超声换能器阵元排列成曲率半径为40 mm、覆盖角度为270°的环形阵列。对1只麻醉后的雄性裸鼠（6~8周龄，仰卧位和俯卧位）进行扫描，依次得到头部、胸部和腹部共计1970个切片的扫描数据。按照6∶3∶1的比例划分为训练集、验证集和测试集，增强后的训练集中共包括3900个切片。
2. 方法的性能评价
分别采用计算机仿真、仿体和在体数据验证所提方法的可行性，并评价其性能。在仿真试验中，成像目标的特性参数是已知的，因此该试验旨在定量评价AAR-Net根据光声信号测量数据重建AOED分布图的精度。仿体和在体试验用于验证AAR-Net在真实场景中应用的可行性。
选择TR[3]、SLSC[7]、U-Net后处理法[21]（以下简称U-Net）和基于DGD网络的图像重建方法[18]（以下简称DGD）作为基线方法，与AAR-Net进行比较，评估AAR-Net的性能。其中DGD网络包含4个结构单元，U-Net输入的初始图像采用TR法得到。采用相同的数据集训练AAR-Net、DGD网络和U-Net，DGD与AAR-Net的训练参数如下：Batch size为64，epoch为200；U-Net的训练参数为：Batch size为1，epoch为200。U-Net、DGD和AAR-Net的训练时间分别是2.4、6.5、8.9 h。
在仿真试验中，以前向仿真得到的AOED分布图作为标准图像，采用峰值信噪比（peak signal to noise ratio，PSNR）和结构相似度（structural similarity，SSIM）作为量化指标，客观评价重建图像的质量。在仿体和在体试验中无法获知AOED的真实值，因此采用对比度（contrast ratio，CR）和对比噪声比（contrast-to-noise ratio，CNR）作为量化指标评价图像质量。
本研究试验中的计算机配置为AMD Ryzen 7 4800H CPU、8GB RAM、Radeon Graphics和Windows 10 64位操作系统。神经网络的搭建、训练和测试的运行环境为Ubuntu 16.04，所用GPU是NVIDIA公司生产的GeForce 3090Ti，显存为16 GB。编程语言是Python 3.7，深度学习框架是TensorFlow 2.6.0和Keras 2.0。
3. 试验结果与分析
3.1 仿真试验结果
分别采用不同方法对仿真光声信号进行图像重建，部分结果如图6所示。由图6可知，采用TR法重建的图像中存在明显的模糊和畸变。SLSC法重建的图像中存在部分组织被当作杂波清除掉的问题，这是由于致密组织产生的非相干信号的空间相干性类似于或者低于杂波，采用SLSC法无法对其进行有效区分。DGD网络重建的图像中，管壁区域存在少量模糊，对复杂成分组织区域的重建效果不理想。经过U-Net处理的图像中，在管腔边界处仍然存在模糊，导致边界不清晰。AAR-Net的重建图像质量最优，组织边界最清晰，且对于成分复杂的组织区域也表现出相对较好的重建效果，更接近标准图像。
图 6 根据仿真声压信号重建图像的结果Fig. 6 Results of image reconstruction from simulated acoustic pressure signals下载:
全尺寸图片
3.2 仿体试验结果
对仿体的图像重建结果如图7所示，与照片对比，采用AAR-Net、DGD、U-Net、SLSC和TR方法重建的图像都能够显示嵌入目标的位置。但在DGD、U-Net、SLSC和TR重建图像中，嵌入物周围存在伪影，而AAR-Net重建的图像具有最佳的视觉效果和最高的对比度。由定量评价结果可知，AAR-Net的重建误差更小且结构准确度更高。
图 7 仿体图像重建结果Fig. 7 Results of image reconstruction for phantoms下载:
全尺寸图片
3.3 在体试验结果
图8是包含复杂结构的活体小鼠4个胸腹切片的图像重建结果。由图8可知，与仿真和仿体研究类似，采用SLSC和TR法重建图像的整体亮度以及横截面边缘的对比度都明显低于深度学习方法重建的图像。在AAR-Net重建图像中，内部结构区分度更高，轮廓更清晰，箭头标示的组织细节的对比度均高于其他方法的重建图像，证明AAR-Net重建复杂精细目标的性能优于非学习方法以及DGD和U-Net方法。
图 8 在体小鼠全身扫描图像重建结果Fig. 8 Results of image reconstruction from in vivo whole-body scanning data of mice下载:
全尺寸图片
4. 讨论
4.1 DGD单元数对重建图像质量的影响
为了分析AAR-Net中的DGD结构单元数对网络性能的影响，在其他条件和参数不变的情况下，改变DGD单元数，对仿真光声信号进行图像重建，结果如图9所示，训练网络的损失函数曲线如图10所示。由图可知，增加DGD单元可以改善重建图像的质量评价指标，但是这种改善并不显著，图像视觉效果的差异也不明显。这是由于AAR-Net采用了先初始重建再优化的策略，DGD模块输出的图像还需输入U-Net模块进行优化，因而改变DGD单元数并不会对最终的图像质量造成很大影响。但是增加DGD单元会显著延长图像重建时间，如表2所列。为了取得重建质量和时间成本的折中，本研究最终采用1个DGD单元构建AAR-Net。
图 9 采用由不同数目的DGD单元构成的AAR-Net根据仿真信号重建图像的结果Fig. 9 Results of image reconstruction from simulated signals using AAR-Net consisting of different number of DGD units下载:
全尺寸图片
图 10 训练由不同数量的DGD单元组成的AAR-Net的损失曲线Fig. 10 Loss of training AAR-Net consisting of different number of DGD units下载:
全尺寸图片
表 2 分别采用包含不同数量的DGD单元的AAR-Net重建一幅（256像素×256像素）像素图像所需的时间Table 2 The time required to reconstruct a （256×256） pixel image using AAR-Net with different number of DGD unitss DGD单元数 重建时间 1 10.6 2 13.4 3 16.2 4 19.0 4.2 所提方法的优势和不足
在光声成像中，由组织不均匀声学特性所导致的图像质量下降问题是光声波在介质中传播时的反射、散射、吸收和扩散等多种声学效应共同作用的结果。与吸收和扩散相比，反射和散射是引起声衰减的主要因素。目前的解决方法多依赖于改善成像系统，例如：杂波去相关法需要用超声探头触诊待成像组织使其产生变形[6]；双模态成像法需要在光声成像前进行超声扫描，定量测定组织的声学特性[9-10]；多波长激励法需要分别利用不同波长的短脉冲激光照射样品得到光声光谱图像序列，并分析图像中各吸收体的光谱响应等[8]。与单纯进行光声扫描相比，额外的扫描过程明显延长了图像采集时间，提高了成像过程的复杂度。本研究方法基于图像增强的策略，可以在不改变成像装置的前提下补偿由于反射或散射等因素所致的光声信号失真，抑制与组织声学异质性相关的图像伪影，重建高质量的AOED分布图。本研究方法降低了对扫描和检测设备的要求，使构建更紧凑的成像系统成为可能。此外，与基于模型的图像重建方法相比，本研究方法无需任何与成像目标特性相关的先验知识，不需要选择正则化工具，即可完成从“信号”到“图像”的映射。该方法可以推广用于光声成像的幂律衰减补偿，不受探测几何形状的限制。
除了由于组织的复杂声学特性导致的反射和散射伪影之外，未优化的光学器件导致的模糊伪影、探测器的有限带宽伪影以及有限的探测视角导致的有限视角伪影都会影响图像的视觉效果，降低成像精度。本研究方法在前向成像算子中充分考虑了组织的不均匀声学特性，但并未包含其他非理想成像条件，如非均匀光覆盖、探测器的有限角度扫描、有限带宽和脉冲响应等。为了提高成像质量，后续工作中需要全面考虑实际成像场景中的非理想情况，避免产生伪影的因素。
5. 结束语
本研究提出的光声图像重建方法解决了由于介质的不均匀声学特性所致的图像质量下降问题。所搭建的神经网络由DGD模块和U-Net模块组成，其中前者实现了从原始光声信号到AOED分布图的映射，将前向成像算子和其伴随算子嵌入到网络中，但算子的计算过程与网络训练分开，降低了网络的复杂性，减少了所需的测量数据；后者用于优化初始重建图像，最后输出高质量的图像。仿真、仿体以及在体试验结果表明，采用所提方法重建图像的视觉效果和评价指标均优于TR、SLSC、U-Net后处理法和基于DGD的图像重建方法。讨论结果表明，虽然增加AAR-Net中的DGD结构单元数，可使重建图像质量略有改善，但同时也会大大提高时间成本，需要在图像质量和时间成本之间做出权衡。在未来的工作中，考虑搭建体系结构更为复杂的DNN，将光声信号的去噪、图像重建和图像优化整合起来，根据去噪后的光声信号重建初始图像，再利用图像增强技术对其进行后处理，实现从信号到高质量图像的映射。
-
图 1 本研究方法流程
Fig. 1 Flow chart of the proposed method
下载:
全尺寸图片
图 2 AAR-Net体系结构
Fig. 2 Architecture of AAR-Net
下载:
全尺寸图片
图 3 网络的损失函数曲线
Fig. 3 Loss of network training
下载:
全尺寸图片
图 4 仿真数据集中的样本图像示例
Fig. 4 Examples of samples in simulation datasets
下载:
全尺寸图片
图 5 仿体的实物照片
Fig. 5 Photos of phantoms
下载:
全尺寸图片
图 6 根据仿真声压信号重建图像的结果
Fig. 6 Results of image reconstruction from simulated acoustic pressure signals
下载:
全尺寸图片
图 7 仿体图像重建结果
Fig. 7 Results of image reconstruction for phantoms
下载:
全尺寸图片
图 8 在体小鼠全身扫描图像重建结果
Fig. 8 Results of image reconstruction from in vivo whole-body scanning data of mice
下载:
全尺寸图片
图 9 采用由不同数目的DGD单元构成的AAR-Net根据仿真信号重建图像的结果
Fig. 9 Results of image reconstruction from simulated signals using AAR-Net consisting of different number of DGD units
下载:
全尺寸图片
图 10 训练由不同数量的DGD单元组成的AAR-Net的损失曲线
Fig. 10 Loss of training AAR-Net consisting of different number of DGD units
下载:
全尺寸图片
表 1 数值仿体中不同组织成分的特性参数
Table 1 Characteristic parameters of different tissue components in numerical phantoms
组织名称 折射率 反射率 吸收系数/cm‒1 散射系数/cm‒1 各向异性因子 声速/(m·s‒1) 密度/(kg·L‒1) 径向厚度/mm 外壁 1.42 0.08 0.70 6 0.85 1600 1.026 0.00~0.50 内壁 1.42 0.05 0.40 6 0.85 1580 1.073 0.10~0.40 纤维帽 1.47 0.24 0.78 25 0.78 1610 1.058 0.20~0.50 钙化 1.47 0.34 0.83 560 0.78 1540 1.668 0.30~0.50 脂质池 1.47 0.42 0.90 520 0.82 1650 0.947 0.10~0.30 混合钙化 1.47 0.23 0.96 550 0.80 1650 0.947 0.01~0.30 内腔 1.35 0.12 0 610 0.99 1540 1.065 1.20~1.70 表 2 分别采用包含不同数量的DGD单元的AAR-Net重建一幅（256像素×256像素）像素图像所需的时间
Table 2 The time required to reconstruct a （256×256） pixel image using AAR-Net with different number of DGD units
s DGD单元数 重建时间 1 10.6 2 13.4 3 16.2 4 19.0 -
[1] 李娇, 李帅, 陈冀景, 等. 非接触光声成像研究进展及其在生物医学上的应用[J]. 中国激光, 2021, 48(19): 284–302. LI Jiao, LI Shuai, CHEN Jijing, et al. Progress and biomedical application of non-contact photoacoustic imaging[J]. Chinese journal of lasers, 2021, 48(19): 284–302. [2] 王雅慧, 王浩全, 任时磊. 滤波反投影算法在光声层析成像中的应用[J]. 国外电子测量技术, 2020, 39(5): 1–4. doi: 10.19652/j.cnki.femt.2002001 WANG Yahui, WANG Haoquan, REN Shilei. Application of filter back projection algorithm in photoacoustic tomography[J]. Foreign electronic measurement technology, 2020, 39(5): 1–4. doi: 10.19652/j.cnki.femt.2002001 [3] 韩朵朵, 孙正, 苑园. 血管内光声图像的时间反演重建方法[J]. 中国图象图形学报, 2016, 21(4): 442–450. doi: 10.11834/jig.20160405 HAN Duoduo, SUN Zheng, YUAN Yuan. Time-reversal based reconstruction of intravascular photoacoustic images[J]. Journal of image and graphics, 2016, 21(4): 442–450. doi: 10.11834/jig.20160405 [4] POUDEL J, LOU Yang, ANASTASIO M A. A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography[J]. Physics in medicine & biology, 2019, 64(14): 14TR01. [5] 段爽, 孙正. 声学特性不均匀组织的光声层析图像重建研究进展[J]. 生物医学工程学杂志, 2019, 36(3): 486–492. DUAN Shuang, SUN Zheng. Review on photoacoustic tomographic image reconstruction for acoustically heterogeneous tissues[J]. Journal of biomedical engineering, 2019, 36(3): 486–492. [6] PETROSYAN T, THEODOROU M, BAMBER J, et al. Rapid scanning wide-field clutter elimination in epi-optoacoustic imaging using comb LOVIT[J]. Photoacoustics, 2018, 10: 20–30. doi: 10.1016/j.pacs.2018.02.001 [7] LEDIJU B M A, KUO N, SONG D Y, et al. Short-lag spatial coherence beamforming of photoacoustic images for enhanced visualization of prostate brachytherapy seeds[J]. Biomedical optics express, 2013, 4(10): 1964–1977. doi: 10.1364/BOE.4.001964 [8] NGUYEN H N Y, HUSSAIN A, STEENBERGEN W. Reflection artifact identification in photoacoustic imaging using multi-wavelength excitation[J]. Biomedical optics express, 2018, 9(10): 4613. doi: 10.1364/BOE.9.004613 [9] KUNIYIL A S M, STEENBERGEN W. Photoacoustic-guided focused ultrasound (PAFUSion) for identifying reflection artifacts in photoacoustic imaging[J]. Photoacoustics, 2015, 3(4): 123–131. doi: 10.1016/j.pacs.2015.09.001 [10] SCHWAB H M, BECKMANN M F, SCHMITZ G. Photoacoustic clutter reduction by inversion of a linear scatter model using plane wave ultrasound measurements[J]. Biomedical optics express, 2016, 7(4): 1468. doi: 10.1364/BOE.7.001468 [11] HIRSCH L, GONZÁLEZ M G, REY V L. On the robustness of model-based algorithms for photoacoustic tomography: comparison between time and frequency domains[J]. Review of scientific instruments, 2021, 92(11): 114901. doi: 10.1063/5.0065966 [12] DEÁN-BEN X L, ÖZBEK A, LÓPEZ-SCHIER H, et al. Acoustic scattering mediated single detector optoacoustic tomography[J]. Physical review letters, 2019, 123(17): 174301. doi: 10.1103/PhysRevLett.123.174301 [13] YIN Jie, TAO Chao, CAI Peng, et al. Photoacoustic tomography based on the Green’s function retrieval with ultrasound interferometry for sample partially behind an acoustically scattering layer[J]. Applied physics letters, 2015, 106(23): 234101. doi: 10.1063/1.4922386 [14] RUI Wei, TAO Chao, LIU Xiaojun. Photoacoustic imaging in scattering media by combining a correlation matrix filter with a time reversal operator[J]. Optics express, 2017, 25(19): 22840. doi: 10.1364/OE.25.022840 [15] RAJENDRAN P, SHARMA A, PRAMANIK M. Photoacoustic imaging aided with deep learning: a review[J]. Biomedical engineering letters, 2022, 12(2): 155–173. doi: 10.1007/s13534-021-00210-y [16] MANWAR R, ZAFAR M, XU Qiuyun. Signal and image processing in biomedical photoacoustic imaging: a review[J]. Optics, 2020, 2(1): 1–24. doi: 10.3390/opt2010001 [17] FARNIA P, MOHAMMADI M, NAJAFZADEH E, et al. High-quality photoacoustic image reconstruction based on deep convolutional neural network: towards intra-operative photoacoustic imaging[J]. Biomedical physics & engineering express, 2020, 6(4): 045019. [18] SUN Zheng, WANG Xinyu, YAN Xiangyang. An iterative gradient convolutional neural network and its application in endoscopic photoacoustic image formation from incomplete acoustic measurement[J]. Neural computing and applications, 2021, 33(14): 8555–8574. doi: 10.1007/s00521-020-05607-x [19] JOSEPH F K, ARORA A, KANCHARLA P, et al. Generative adversarial network-based photoacoustic image reconstruction from bandlimited and limited-view data[C]//Proc SPIE 11642, Photons Plus Ultrasound: Imaging and Sensing. [S.l.: s.n.], 2021: 208−213. [20] VU T, LI Mucong, HUMAYUN H, et al. A generative adversarial network for artifact removal in photoacoustic computed tomography with a linear-array transducer[J]. Experimental biology and medicine, 2020, 245(7): 597–605. doi: 10.1177/1535370220914285 [21] SUN Zheng, YAN Xiangyang. A deep learning method for limited-view intravascular photoacoustic image reconstruction[J]. Journal of medical imaging and health informatics, 2020, 10(11): 2707–2713. doi: 10.1166/jmihi.2020.3204 [22] GODEFROY G, ARNAL B, BOSSY E. Compensating for visibility artefacts in photoacoustic imaging with a deep learning approach providing prediction uncertainties[J]. Photoacoustics, 2021, 21: 100218. doi: 10.1016/j.pacs.2020.100218 [23] 宋宜昌, 邢达. 用蒙特卡罗方法模拟光在多层组织中的吸收特性[J]. 激光生物学报, 2002, 11(1): 10–14. doi: 10.3969/j.issn.1007-7146.2002.01.003 [24] TREEBY B E, COX B T. Modeling power law absorption and dispersion for acoustic propagation using the fractional Laplacian[J]. The journal of the acoustical society of America, 2010, 127(5): 2741–2748. doi: 10.1121/1.3377056 [25] LA RIVIÈRE P J, ZHANG Jin, ANASTASIO M A. Image reconstruction in optoacoustic tomography for dispersive acoustic media[J]. Optics letters, 2006, 31(6): 781. doi: 10.1364/OL.31.000781 [26] MOHAMMADI L, BEHNAM H, TAVAKKOLI J, et al. Skull’s photoacoustic attenuation and dispersion modeling with deterministic ray-tracing: towards real-time aberration correction[J]. Sensors, 2019, 19(2): 345. doi: 10.3390/s19020345 [27] 齐梦雨, 严壮志, 赵丽丽. 基于有限元和k空间伪谱法的光声成像仿真平台[J]. 中国医疗器械杂志, 2018, 42(6): 413–416,427. QI Mengyu, YAN Zhuangzhi, ZHAO Lili. Simulation platform of photoacoustic imaging based on finite element and k-space pseudospectral method[J]. Chinese journal of medical instrumentation, 2018, 42(6): 413–416,427. [28] ADLER J, ÖKTEM O. Solving ill-posed inverse problems using iterative deep neural networks[J]. Inverse problems, 2017, 33(12): 124007. doi: 10.1088/1361-6420/aa9581 [29] KINGMA D P, BA J. Adam: a method for stochastic optimization[EB/OL]. (2014−12−22)[2022−01−01]. https://arxiv.org/abs/1412.6980.pdf. [30] TREEBY B E, COX B T. K-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields[J]. Journal of biomedical optics, 2010, 15(2): 021314. doi: 10.1117/1.3360308
