Contents Introduction 引言
Real-time integrated photoacoustic and ultrasound (PAUS) imaging is a promising approach to bring the molecular sensitivity of optical contrast mechanisms into clinical ultrasound (US) systems. Laser pulses transmitted into tissue induce light absorption from endogenous chromophores or exogenous contrast agents, which in turn launch acoustic waves according to the photoacoustic (PA) effect that can be used for imaging. We have recently developed a customized system for simultaneous PA and US imaging using interleaved techniques at fast scan rates [1]. One of the potential clinical applications is real-time, quantitative monitoring of blood oxygenation in the microvasculature [2], [3]. It requires not only multiple measurements at different optical wavelengths, but also high image quality to preserve microvascular topology.
实时集成的光声和超声（PAUS）成像是一种有前景的方法，可以将光学对比机制的分子敏感性引入临床超声系统（US）。激光脉冲传输到组织中会引起内源性色素或外源性对比剂的光吸收，进而根据光声（PA）效应启动声波，可用于成像。我们最近开发了一种用于同时进行 PA 和US成像的定制系统，使用交错技术以快速扫描率 [1] 。其中一个潜在的临床应用是实时、定量监测微血管中的血氧饱和度 [2] ， [3] 。这不仅需要在不同的光学波长进行多次测量，还需要高质量的图像以保持微血管拓扑结构。
Similar to US beamforming, PA reconstruction widely uses a traditional delay-and-sum (DAS) algorithm [4] for simplicity. However, the limited view and relatively narrow bandwidth of clinical US arrays greatly degrade image quality due to the ultra-broad bandwidth nature of PA signals [5]. This ill-posed problem causes structure loss, low contrast, and diverse artifacts making image interpretation difficult.
类似于US波束形成，PA 重建广泛使用传统的延迟和求和（DAS）算法 [4] ，因为它简单。然而，临床US阵列的有限视野和相对狭窄的带宽极大地降低了图像质量，这是由于 PA 信号的超宽带宽特性 [5] 。这个问题导致结构丢失、对比度低和多样的伪影，使得图像解释困难。
To address these challenges, many groups have adopted reconstruction techniques from US or radar imaging to PA imaging. In particular, adaptive approaches such as a Minimum Variance (MV) method was developed to reduce off-axis signal and sidelobe artifacts by assigning apodization weights based on statistics [6], [7]. Reconstruction using Delay Multiply and Sum (DMAS) methods can achieve higher image contrast by enhancing signal coherence nonlinearly [8], [9]. Also, iterative techniques for the inverse reconstruction problem have been consistently developed for PA imaging using signal sparsity and low-rankness [10], [11]. All these methods may improve image quality by adopting sophisticated models based on system physics, data statistics, and underlying object properties. However, the main drawbacks are high computational complexity and the requirement of carefully selecting a handful of parameters.
为了应对这些挑战，许多团队采用了从雷达成像到光声成像的重建技术。特别是，开发了最小方差（MV）方法作为一种自适应方法，通过基于统计数据分配权重来减少轴外信号和旁瓣伪影。使用延迟乘积和求和（DMAS）方法的重建可以通过非线性增强信号一致性来实现更高的图像对比度。此外，针对光声成像的逆重建问题，一直在持续开发基于信号稀疏性和低秩性的迭代技术。所有这些方法都可能通过采用基于系统物理、数据统计和底层对象属性的复杂模型来提高图像质量。然而，主要缺点是计算复杂度高和需要仔细选择少数参数。
Currently, a new category of reconstruction methods has been inspired by the field of machine learning (ML). Due to the great success of deep convolutional neural networks (CNN) in computer vision, reconstruction using supervised learning is an emerging research area in medical imaging [12]. The network extracts best features via learning weights/filters to mitigate ill-posedness in inverse problems. The most popular ML framework uses learning in the image-domain, where training inputs are corrupted images processed by a standard reconstruction method under ill-posed conditions [13]–[17]. Using this approach, the network avoids trying to capture detailed reconstruction operations and concentrates on filtering artifacts and noise. However, since the details of actual detector data are lost after image reconstruction, ML applied to reconstructed images often cannot recover weak signals and fine structures can be lost.
目前，机器学习（ML）领域的启发已催生了一种新的重建方法类别。由于深度卷积神经网络（CNN）在计算机视觉中的巨大成功，使用监督学习进行重建已成为医学成像中的一个新兴研究领域 [12] 。网络通过学习权重/滤波器来提取最佳特征，以减轻逆问题中的不适定性。最流行的 ML 框架使用图像域中的学习，其中训练输入是由标准重建方法处理的在不适定条件下的损坏图像 [13] – [17] 。采用这种方法，网络避免尝试捕捉详细的重建操作，专注于过滤伪影和噪声。然而，由于实际探测器数据在图像重建后丢失，应用于重建图像的 ML 通常无法恢复微弱信号，细微结构可能会丢失。
Some frameworks are based on iterative schemes to train regularizations in Compressed Sensing (CS), but their extensive computations restrict real-time clinical applications [18], [19]. Zhu et al proposed a framework starting with acquired data without prior knowledge of physics, but a fully connected layer requires a large number of weighting parameters for large data sets [20]. Allman et al employed PA raw data, but the application was only limited to the classification of point-like targets from artifacts [21].
一些框架基于迭代方案来训练压缩感知（CS）中的正则化，但它们的大量计算限制了实时临床应用。朱等人提出了一个框架，从获取的数据开始，不需要物理学的先验知识，但是一个完全连接的层需要大数据集的大量权重参数。Allman等人使用了PA原始数据，但应用仅限于从伪影中分类点状目标。
Here we explore practical PA image reconstruction based on a deep-learning technique suitable for real-time PAUS imaging. We first examine the link between model-based methods and basic neural network layers to help design and interpret the learning structure. As discussed below, this study led us to modify 2-D raw data (with time and detector dimensions) into a 3-D array (with two spatial dimensions and a channel dimension), where a channel packet corresponds to the propagation delay profile for one spatial point, as an input to the neural network. The delay operation simplifies the learning process and the extension to channel dimension retains more information and increases learning accuracy.
在这里，我们探讨了基于深度学习技术的实用光声成像重建方法，适用于实时光声超声成像。我们首先检验了基于模型的方法与基本神经网络层之间的联系，以帮助设计和解释学习结构。如下所述，这项研究使我们将二维原始数据（具有时间和探测器维度）修改为三维数组（具有两个空间维度和一个通道维度），其中一个通道包对应于一个空间点的传播延迟轮廓，作为神经网络的输入。延迟操作简化了学习过程，而通道维度的扩展保留了更多信息并提高了学习精度。
Our subsequent architecture is based on U-net [22], where dyadic scale decomposition can access data in multi-resolution support. The structure can extract comprehensive features from the transformed 3-D array, replacing hand-crafted functions and generalizing standard filtering techniques. For training, we restrict the scope of absorber types to microvessels and create synthetic datasets using simulation. Operators transforming ground-truth to radio frequency (RF) array data are based on our current fast-swept PAUS system [1], i.e. take into account the spectral bandwidth and geometry of a real imaging probe. However, it is not limited to only one imaging system and can be applied to any PA system with known geometry and characteristics.
我们随后的架构基于U-net，其中二进制尺度分解可以访问多分辨率支持的数据。该结构能从变换后的三维数组中提取全面的特征，替代手工制作的功能，并推广标准的过滤技术。在训练中，我们将吸收体类型的范围限制在微血管，并使用模拟创建合成数据集。将地面真实数据转换为射频（RF）数组数据的操作基于我们当前的快速扫描PAUS系统，即考虑了真实成像探头的光谱带宽和几何形状。然而，它不仅限于一个成像系统，还可以应用于任何具有已知几何形状和特性的PA系统。
To demonstrate the performance of the CNN-based method, we first compare it to standard methods using synthetic data. Then, we performed phantom experiments using the fast-swept PAUS system, and finally imaged a human finger in vivo.
为了展示基于 CNN 方法的性能，我们首先使用合成数据与标准方法进行比较。然后，我们使用快速扫描的 PAUS 系统进行了幻影实验，并最终对人体手指进行了体内成像。
Signal and System Model 信号与系统模型
A. PA Forward Operation A. PA 前向操作
The spatiotemporal pressure change
\begin{equation*} \left({\nabla ^{2} - \frac {1}{v_{s}^{2}} {\frac {\partial ^{2}}{\partial t^{2}}}}\right) p(t, \mathbf {r}) = -\frac {\beta }{C_{p}}\frac {\partial H(t, \mathbf {r})}{\partial t},\tag{1}\end{equation*}
where \begin{equation*} p(t, \mathbf {r}') = \frac {\Gamma }{4\pi v_{s}^{2}} \frac {\partial }{\partial t} \left[{ \int \frac {d \mathbf {r}}{\vert \mathbf {r} - \mathbf {r}'\vert } s(\mathbf {r}) \delta \left({t-\frac {\vert \mathbf {r} - \mathbf {r}'\vert }{v_{s}}}\right) }\right]\tag{2}\end{equation*}
\begin{equation*} {y(t, \mathbf {r}'_{j}) = h(t)*\big [p_{d}(t, \mathbf {r}'_{j})+ n(t, \mathbf {r}'_{j})\big] }\tag{3}\end{equation*}
短激光脉冲激发下的时空压力变化
\begin{equation*} \left({\nabla ^{2} - \frac {1}{v_{s}^{2}} {\frac {\partial ^{2}}{\partial t^{2}}}}\right) p(t, \mathbf {r}) = -\frac {\beta }{C_{p}}\frac {\partial H(t, \mathbf {r})}{\partial t},\tag{1}\end{equation*}
\begin{equation*} p(t, \mathbf {r}') = \frac {\Gamma }{4\pi v_{s}^{2}} \frac {\partial }{\partial t} \left[{ \int \frac {d \mathbf {r}}{\vert \mathbf {r} - \mathbf {r}'\vert } s(\mathbf {r}) \delta \left({t-\frac {\vert \mathbf {r} - \mathbf {r}'\vert }{v_{s}}}\right) }\right]\tag{2}\end{equation*}
\begin{equation*} {y(t, \mathbf {r}'_{j}) = h(t)*\big [p_{d}(t, \mathbf {r}'_{j})+ n(t, \mathbf {r}'_{j})\big] }\tag{3}\end{equation*}
B. Limitations of Handheld Linear Array System
手持线性阵列系统的局限性
Reconstruction methods have been developed based on the measurement geometry. Filtered back-projection (FBP) is derived from the inverse of the PA forward operation in the spatiotemporal domain or k-space (frequency) domain [25]. Exact FBP formulas were demonstrated for a 3-D absorber distribution where the detection geometry is spherical, planar or cylindrical [26]. Imaging for a 2D spatial plane (slice) was adopted in standard tomographic scans assuming the detector (transducer) is focused in the plane [27], [28]. A circular transducer for a 2D source distribution produces an accurate reconstruction if the number of detector elements provides enough spatial sampling density, as shown in Fig. 1 (B) [15], [29]–[31], and their bandwidth is not limited. In contrast, a linear sensor geometry (
基于测量几何已经开发了重建方法。滤波反投影（FBP）是从时空域或 k 空间（频率）域的 PA 正向操作的逆推导出的 [25] 。对于检测几何是球形、平面或圆柱形的 3-D 吸收体分布，已经演示了精确的 FBP 公式 [26] 。在标准断层扫描中采用了 2D 空间平面（切片）的成像，假设探测器（换能器）聚焦在平面上 [27] ， [28] 。一个圆形换能器对于 2D 源分布产生准确的重建，如果探测器元素的数量提供足够的空间采样密度，如 Fig. 1 (B) [15] ， [29] – [31] 所示，且其带宽不受限制。相比之下，线性传感器几何（
Simulation results using standard filtered back-projection reconstruction. (A) presents two example object shapes. (B-D) shows reconstructions when the measurement conditions are (B) circular array with full bandwidth, (C) linear array with full bandwidth, and (D) linear array with limited bandwidth (11-19 MH). Arrows indicate structural losses and artifacts. All images are visualized on a log-scale colormap (40 dB range).
使用标准滤波反投影重建的模拟结果。(A)展示了两个示例物体形状。(B-D)展示了测量条件下的重建结果：(B)圆形阵列全带宽，(C)线性阵列全带宽，以及(D)线性阵列有限带宽(11-19 MH)。箭头指示了结构损失和伪影。所有图像均在对数尺度色图上可视化（40 dB范围）。
This ill-posed problem for a finite bandwidth linear array is more understandable in the frequency domain. For simplicity, assume that the directivity function is constant and the noise power is zero. Then, the 2-D Fourier transform of Eq. 3 for this geometry can be represented as [32]
\begin{equation*} Y(f,k_{x}) = \alpha \frac {f sgn(f)}{\sqrt {\left({\frac {f}{c}}\right)^{2}-k_{z}^{2}}} S(k_{z}, k_{x})\tag{4}\end{equation*}
where 对于有限带宽线性阵列的这一病态问题，在频域中更易于理解。为简化问题，假设指向性函数是常数，噪声功率为零。那么，对于这种几何形状， Eq. 3 的二维傅里叶变换可以表示为 [32]
\begin{equation*} Y(f,k_{x}) = \alpha \frac {f sgn(f)}{\sqrt {\left({\frac {f}{c}}\right)^{2}-k_{z}^{2}}} S(k_{z}, k_{x})\tag{4}\end{equation*}
Reconstruction simulations in k-space using one circular object and two simple linear objects rotated by 90 degrees. All images visualize absolute pixel values on a log-scale colormap (40 dB range). The maximum value in each image is represented as pure white. (A) Ground-truth images. (B) K-domain-GT obtained by 2-D Fourier transforming the ground truth. (C) K-domain-data obtained by 2-D Fourier transforming raw data. Raw x-t data obtained with the forward model. The dotted lines indicate
使用一个圆形物体和两个简单的线性物体（旋转 90 度）在 k 空间中进行重建模拟。所有图像在对数比例颜色图（40 dB 范围）上显示绝对像素值。每幅图像中的最大值显示为纯白色。（A）基准真实图像。（B）通过二维傅里叶变换基准真实得到的 K 域-GT。（C）通过二维傅里叶变换原始数据得到的 K 域数据。使用前向模型得到的原始 x-t 数据。虚线表示
A special case of a vertical line source (bottom line in Fig. 2) exaggerates the problem. Its k-space spectrum is almost totally filtered, and only low amplitude spectral sidelobes invisible in the ground truth image survive. Thus, only top and bottom source points are visualized in the reconstruction.
垂直线源的一个特殊案例（ Fig. 2 中的底线）夸大了问题。其 k 空间谱几乎完全被滤除，只有在基准真实图像中不可见的低幅度谱侧瓣存活。因此，只有顶部和底部的源点在重建中可视化。
In this paper, the target objects of interest are microvessels, where the shape can be represented as a sum of straight and curvy lines. Since the signal components of the typical vascular structure are widely distributed in the k-domain, and a sufficient fraction are maintained even after limited view/bandwidth induced filtration, there is the possibility to reconstruct the entire object shape. However, note that it will be very challenging to recover vertical portions, as shown in Fig. 1 and Fig. 2.
在本文中，感兴趣的目标对象是微血管，其形状可以表示为直线和曲线的总和。由于典型血管结构的信号成分在 k 域中广泛分布，并且即使在受限视图/带宽诱导的过滤后仍保持足够的比例，因此有可能重建整个对象形状。然而，请注意，恢复垂直部分将是非常具有挑战性的，如 Fig. 1 和 Fig. 2 所示。
Conventional Methods in PA Image Reconstruction
PA 图像重建中的传统方法
A. Delay and Sum (DAS)
A. 延迟与求和 (DAS)
Most commercial US systems use DAS beamforming [4] for real-time image reconstruction. This procedure applies a delay due to the propagation distances between an observation point in the image and each transducer element prior to summation of signals across all array elements. Likewise, since the PA signal is based on one-way acoustic propagation, the DAS method can be applied as
\begin{equation*} \tilde {s} (\mathbf {r}) = \sum _{j} w(\mathbf {r}, \mathbf {r}'_{j}) y \left({\frac {\vert \mathbf {r} - \mathbf {r}'_{j}\vert }{v_{s}}, \mathbf {r}'_{j}}\right)\tag{5}\end{equation*}
where \begin{equation*} \tilde {s} (\mathbf {r}) = \sum _{j} w(\mathbf {r}, \mathbf {r}'_{j}) f (\mathbf {r}, j)\tag{6}\end{equation*}
大多数商业US系统使用 DAS 波束形成 [4] 进行实时图像重建。这一过程涉及到在信号跨所有阵列元素求和之前，由于图像中观察点与每个换能器元素之间的传播距离而产生的延迟。同样，由于 PA 信号基于单向声波传播，DAS 方法可以作为
\begin{equation*} \tilde {s} (\mathbf {r}) = \sum _{j} w(\mathbf {r}, \mathbf {r}'_{j}) y \left({\frac {\vert \mathbf {r} - \mathbf {r}'_{j}\vert }{v_{s}}, \mathbf {r}'_{j}}\right)\tag{5}\end{equation*}
\begin{equation*} \tilde {s} (\mathbf {r}) = \sum _{j} w(\mathbf {r}, \mathbf {r}'_{j}) f (\mathbf {r}, j)\tag{6}\end{equation*}
(A) Measurement geometry. A 2-D image plane with respect to a linear array transducer is defined as the z-x plane. (B) 2-D measurement data. The curved lines indicate propagation delay profiles of particular image points at different depths. (C) 3-D transformed data. Channel packets correspond to the delay profiles indicated by straight lines.
（A）测量几何。相对于线性阵列换能器定义为 z-x 平面的 2-D 图像平面。（B）2-D 测量数据。曲线表示不同深度的特定图像点的传播延迟轮廓。（C）3-D 转换数据。通道包对应于由直线指示的延迟轮廓。
B. Minimum Variance (MV) B. 最小方差（MV）
MV is also based on Eq. 6 but the weighting is adaptive [6]. For a position
\begin{equation*} \min _{\bar { \mathbf {w}}} \bar { \mathbf {w}}^{T} E[\mathbf {f}_{i} \mathbf {f}_{i}^{T}] \bar { \mathbf {w}} \;\; s.t. \;\; \bar { \mathbf {w}}^{T} \mathbf {1} = 1,\tag{7}\end{equation*}
where \begin{equation*} \bar { \mathbf {w}} = \frac { \mathbf {R}_{i}^{-1} \mathbf {1}}{\mathbf {1}^{H} \mathbf {R}_{i}^{-1} \mathbf {1}}\tag{8}\end{equation*}
\begin{equation*} \hat { \mathbf {R}}_{i} = \frac {1}{(2N+1)(J-L)} \sum _{n=-N}^{N}\sum _{l=0}^{J-L} \bar { \mathbf {f}}_{i-n}^{(l)}\bar { \mathbf {f}}_{i-n}^{(l) T}\tag{9}\end{equation*}
MV 同样基于 Eq. 6 ，但权重是自适应的 [6] 。对于一个位置
\begin{equation*} \min _{\bar { \mathbf {w}}} \bar { \mathbf {w}}^{T} E[\mathbf {f}_{i} \mathbf {f}_{i}^{T}] \bar { \mathbf {w}} \;\; s.t. \;\; \bar { \mathbf {w}}^{T} \mathbf {1} = 1,\tag{7}\end{equation*}
\begin{equation*} \bar { \mathbf {w}} = \frac { \mathbf {R}_{i}^{-1} \mathbf {1}}{\mathbf {1}^{H} \mathbf {R}_{i}^{-1} \mathbf {1}}\tag{8}\end{equation*}
\begin{equation*} \hat { \mathbf {R}}_{i} = \frac {1}{(2N+1)(J-L)} \sum _{n=-N}^{N}\sum _{l=0}^{J-L} \bar { \mathbf {f}}_{i-n}^{(l)}\bar { \mathbf {f}}_{i-n}^{(l) T}\tag{9}\end{equation*}
C. Delay Multiply and Sum (DMAS)
延迟乘积求和法（DMAS）
The multi-channel array
\begin{equation*} \tilde {s}(\mathbf {r}_{i}) = \sum _{j_{1},j_{2}} sgn(\bar {f}_{j_{1},j_{2}}(\mathbf {r}_{i}))\sqrt {\vert \bar {f}_{j_{1},j_{2}}(\mathbf {r}_{i})\vert }\tag{10}\end{equation*}
where 多通道阵列
\begin{equation*} \tilde {s}(\mathbf {r}_{i}) = \sum _{j_{1},j_{2}} sgn(\bar {f}_{j_{1},j_{2}}(\mathbf {r}_{i}))\sqrt {\vert \bar {f}_{j_{1},j_{2}}(\mathbf {r}_{i})\vert }\tag{10}\end{equation*}
Schematic diagram illustrating reconstruction methods. 2-D measurement data are transformed into a 3-D array as shown in Fig. 3 followed by reconstruction. (A) DAS / MV methods. An image pixel is determined by weighting and summing channel samples at a corresponding pixel position. Weights vary with position. Unlike DAS, MV adaptively assigns weights depending on data statistics. (B) DMAS method. Channel samples are coupled and multiplied before summation. This additional nonlinear operation is required to prevent a dimensional problem. (C) Iterative method. This is based on the
示意图说明重建方法。2D 测量数据如 Fig. 3 所示转换为 3D 数组，随后进行重建。(A) DAS/MV 方法。图像像素的确定是通过对应像素位置的通道样本进行加权和求和来实现的。权重随位置变化。与 DAS 不同，MV 根据数据统计自适应分配权重。(B) DMAS 方法。通道样本在求和前进行耦合和乘法运算。这种额外的非线性操作是为了防止维度问题。(C) 迭代方法。这是基于压缩感知(CS)中的
D. Iterative Method With Compressed Sensing (CS)
D. 迭代方法与压缩感知（CS）
Iterative methods attempt to solve the inverse problem by adding regularization to overcome the ill-posed condition [10], [11]. The matrix form of Eq. 3 can be expressed as
\begin{equation*} \bar { \mathbf {s}} = arg\min _{ \mathbf {s}} {\vert \vert \mathbf {y}- \mathbf {H} \mathbf {s} \vert \vert _{2}^{2} + \lambda \vert \vert \mathbf {W} ^{T} \mathbf {s}\vert \vert _{1} }\tag{11}\end{equation*}
where \begin{equation*} { \mathbf {s}}_{k+1} = \mathbf {W}\Theta _{\lambda \tau } \big ((\mathbf {W}^{T}-\tau (\mathbf {H} \mathbf {W})^{T} \mathbf {H}) \mathbf {s}_{k} + \tau (\mathbf {H} \mathbf {W})^{T} \mathbf {y}\big)\tag{12}\end{equation*}
迭代方法试图通过添加正则化来解决逆问题，以克服病态条件 [10] ， [11] 。矩阵形式的 Eq. 3 可以表示为
\begin{equation*} \bar { \mathbf {s}} = arg\min _{ \mathbf {s}} {\vert \vert \mathbf {y}- \mathbf {H} \mathbf {s} \vert \vert _{2}^{2} + \lambda \vert \vert \mathbf {W} ^{T} \mathbf {s}\vert \vert _{1} }\tag{11}\end{equation*}
\begin{equation*} { \mathbf {s}}_{k+1} = \mathbf {W}\Theta _{\lambda \tau } \big ((\mathbf {W}^{T}-\tau (\mathbf {H} \mathbf {W})^{T} \mathbf {H}) \mathbf {s}_{k} + \tau (\mathbf {H} \mathbf {W})^{T} \mathbf {y}\big)\tag{12}\end{equation*}
The main disadvantage of MV, DMAS and CS is the high computational complexity for real-time imaging even though modifications have been proposed to reduce the burden. The selection of statistical operators or feature bases is a crucial step in model-based schemes. If the selection does not agree with the inherent properties of PA data, imaging is inaccurate. As an alternative, deep-learning approaches can build optimal feature maps through training and provide a practical reconstruction framework due to fast computation. In addition, they are well-known for strong generalization over parameters such as SNR [36].
MV、DMAS 和 CS 的主要缺点是即使进行了修改以减轻负担，实时成像的计算复杂性仍然很高。在基于模型的方案中，选择统计算子或特征基是关键步骤。如果选择与 PA 数据的固有属性不符，成像将不准确。作为一种替代方案，深度学习方法可以通过训练构建最优特征图，并由于计算速度快，提供一个实用的重建框架。此外，它们以强大的参数泛化能力而闻名，例如 SNR [36] 。
Deep-Learning Reconstruction
深度学习重建
Recently, researchers have begun to find connections between conventional model-based approaches and deep convolutional neural networks (CNN) for inverse problems [13], [33]. We have also explored these links to build an efficient CNN for PA imaging using limited view/bandwidth arrays so that the network takes full advantage of signal characteristics such as row-rank, coherence and sparsity during learning.
最近，研究人员开始发现传统模型基方法与深度卷积神经网络（CNN）在逆问题上的联系。我们也探索了这些联系，以构建一个高效的CNN，用于使用有限视角/带宽阵列的光声成像，使得网络在学习过程中能充分利用信号特性，如行秩、一致性和稀疏性。
A. Preprocessing A. 预处理
It is not clear that a CNN can reconstruct absorption structures directly from PA data. Given the data dimension, learning would be extremely complex since the network architecture must encode the underlying PA forward operation. A popular approach to reduce this burden is preprocessing raw data with simple DAS reconstruction and using the resultant rough images as training input [14], [15]. The CNN can then focus on learning the characteristics of artifacts in input images. However, since these images can lose detailed information on object structure, the CNN output would not be perfect [20].
尚不清楚 CNN 是否能直接从 PA 数据重建吸收结构。考虑到数据维度，学习将极其复杂，因为网络架构必须编码潜在的 PA 前向操作。减轻这一负担的一种流行方法是使用简单的 DAS 重建预处理原始数据，并使用结果粗略图像作为训练输入 [14] ， [15] 。然后 CNN 可以专注于学习输入图像中的伪影特征。然而，由于这些图像可能丢失关于对象结构的详细信息，CNN 的输出将不会完美 [20] 。
Our strategy is to use the transformed 3-D data
我们的策略是使用图 Fig. 3 中展示的转换后的三维数据
This preprocessing is associated with row rank and sparsity in CS. The spatial domains combining axial and lateral dimensions
这种预处理与 CS 中的行秩和稀疏性相关。结合轴向和横向尺寸的空间域
Note that MV and DMAS methods access the coherence using channel-sample correlation to indirectly enhance low-rank (high-coherence) signals while suppressing high-rank (less-coherence) artifacts and noise. The success of these methods suggests that preprocessing would help the CNN to find optimal filter bases (weights) for regression using a realistic number of examples in the training set.
请注意，MV 和 DMAS 方法使用通道-样本相关性间接访问一致性，以增强低秩（高一致性）信号，同时抑制高秩（低一致性）伪影和噪声。这些方法的成功表明，预处理将帮助 CNN 找到用于回归的最佳滤波基（权重），使用训练集中的现实数量的示例。
B. Structure B. 结构
1) Operation and Notation:
1) 操作和符号：
2) Encoder-Decoder: 2) 编码器-解码器：
A standard encoder structure in CNN can be represented as
\begin{equation*} v({ \mathbf {C}_{j}}) = \mathbf {\Phi }^{T} \eta \left({v\left({\sum _{i} \mathbf {F}_{i} \circledast \Psi _{i,j} + \alpha _{j}\mathbf {1}_{n_{1},n_{2}}}\right) }\right)\tag{13}\end{equation*}
where \begin{equation*} \mathbf {Z}_{j} = \eta \left({\sum _{i} v_{n_{1},n_{2}}^{-1}\bigg (\big (\bar {\mathbf {\Phi }} v(\bar { \mathbf {C}}_{i}) \big)\bigg) \circledast \Psi _{i,j} + \beta _{j} \mathbf {1}_{n_{1},n_{2}} }\right)\tag{14}\end{equation*}
CNN 中的标准编码器结构可以表示为
\begin{equation*} v({ \mathbf {C}_{j}}) = \mathbf {\Phi }^{T} \eta \left({v\left({\sum _{i} \mathbf {F}_{i} \circledast \Psi _{i,j} + \alpha _{j}\mathbf {1}_{n_{1},n_{2}}}\right) }\right)\tag{13}\end{equation*}
\begin{equation*} \mathbf {Z}_{j} = \eta \left({\sum _{i} v_{n_{1},n_{2}}^{-1}\bigg (\big (\bar {\mathbf {\Phi }} v(\bar { \mathbf {C}}_{i}) \big)\bigg) \circledast \Psi _{i,j} + \beta _{j} \mathbf {1}_{n_{1},n_{2}} }\right)\tag{14}\end{equation*}
The encoder-decoder convolution layer is similar to the standard reconstruction methods shown in Fig. 4. The common structure is weighting (filtering) channel samples at a location
编码器-解码器卷积层与 Fig. 4 中显示的标准重建方法类似。常见的结构是在位置
Although a CNN layer shares identical weights over space due to the convolution operation, the framework of multi-channel weights per layer and multi-layers increases the expressive power. The iterative method consists of matrix multiplications with no compacting support. While a CNN uses fixed filter size (usually
尽管 CNN 层由于卷积操作在空间上共享相同的权重，但每层的多通道权重和多层的框架增加了表达能力。迭代方法包括矩阵乘法，没有压缩支持。虽然 CNN 使用固定的滤波器大小（通常
Currently, deep learning approaches have been investigated to understand the mathematical framework needed to solve inverse problems. Yin et al proposed the low-rank Hankel matrix approach using a combination of nonlocal basis and local basis for sparse representation of signals [37]. The framelet method has been successfully applied to image processing tasks since matrix decomposition reflects both local and nonlocal behavior of the signal. Ye et al discovered that the encoder-decoder framework of CNN generalizes the framelet representation [38]. In particular, the neural network can decompose the Hankel matrix of 3-D input data and shrink its rank to achieve a rank-deficient ground-truth (See Appendix II).
目前，深度学习方法已被研究以理解解决逆问题所需的数学框架。Yin 等人提出了使用非局部基和局部基的组合的低秩 Hankel 矩阵方法，用于信号的稀疏表示 [37] 。自从矩阵分解反映了信号的局部和非局部行为以来，framelet 方法已成功应用于图像处理任务。Ye 等人发现 CNN 的编码器-解码器框架概括了 framelet 表示 [38] 。特别是，神经网络可以分解 3-D 输入数据的 Hankel 矩阵并缩小其等级，以实现一个等级不足的基准（见 Appendix II ）。
3) Implementation Details:
3) 实施细节：
Our network is based on U-net. Fig. 5 presents the architecture. The left and right sides of the U-shape network correspond to successive encoders and decoders, respectively. Following a
我们的网络基于 U-net。 Fig. 5 展示了架构。U 形网络的左右两侧分别对应连续的编码器和解码器。在一个
Deep-learning architecture for PA image reconstruction. Raw data are converted into a 3-D array by a lookup table (LUT). The array is used as a multi-channel input to the network (first box). Each box represents a multi-channel feature map. The number of channels is denoted on the top or bottom of the box. Feature map sizes decrease and increase via max-pooling and upsampling, respectively. All convolutional layers consist of
PA 图像重建的深度学习架构。原始数据通过查找表（LUT）转换为 3-D 数组。该数组用作网络的多通道输入（第一个框）。每个框代表一个多通道特征图。通道数在框的顶部或底部标注。特征图大小通过最大池化和上采样分别减小和增加。所有卷积层由
C. Network Training C. 网络训练
Training was performed by minimizing a loss function as
\begin{equation*} \hat {\Upsilon } = arg\min _{\Upsilon } {\sum _{k=1}^{K} \vert \vert \mathbf {s}^{(k)}-\Upsilon (\mathbf {F}^{(k)})\vert \vert _{F}^{2} }\tag{15}\end{equation*}
where 通过最小化损失函数来进行训练，其中
表 I 训练网络的参数
D. Training Data D. 训练数据
The supervised learning framework requires data at a large scale. However, it is mostly impractical to obtain clinical raw data accompanied by real ground-truth vascular maps. Therefore, we trained the network by creating synthetic data mimicking typical microvessel networks, as shown in Fig. 5.
监督学习框架需要大规模的数据。然而，获取伴随真实地面真实血管图的临床原始数据通常是不切实际的。因此，我们通过创建模仿典型微血管网络的合成数据来训练网络，如 Fig. 5 所示。
The simulation transforming ground-truth to RF array data is based completely on our PAUS system. The impulse response function
\begin{equation*} {d(\theta _{j})} = \frac {\sin \left({\frac {\pi l}{\lambda } \sin \theta _{j}}\right)}{\frac {\pi l}{\lambda } \sin \theta _{j}}, \;\; \theta _{j} = \tan ^{-1} \left({\frac {x-x_{j}'}{z}}\right),\tag{16}\end{equation*}
where 将地面真实数据转换为 RF 阵列数据的模拟完全基于我们的 PAUS 系统。冲激响应函数
\begin{equation*} {d(\theta _{j})} = \frac {\sin \left({\frac {\pi l}{\lambda } \sin \theta _{j}}\right)}{\frac {\pi l}{\lambda } \sin \theta _{j}}, \;\; \theta _{j} = \tan ^{-1} \left({\frac {x-x_{j}'}{z}}\right),\tag{16}\end{equation*}
Reference vascular images were obtained from the fundus oculi drive [40]. The database contains retina color images captured by camera that can be used for vessel extraction. We used only binary images (manually extracted images), where white pixels denote the segmentation of blood vessels. These images were randomly partitioned, re-sized, rotated and combined with each other to augment the training numbers. Next, every binary image was modified to a gray-scale image where the dynamic range of the vessel signal intensity is 20 dB. Lastly, every image was amplified with different values to obtain measurements involving a wide range of SNR. Table I summarizes all parameters.
参考血管图像是从眼底驱动 [40] 获得的。该数据库包含由相机捕获的视网膜彩色图像，可用于血管提取。我们仅使用二值图像（手动提取的图像），其中白色像素表示血管的分割。这些图像被随机分割、调整大小、旋转并相互结合以增加训练数量。接下来，每个二值图像被修改为灰度图像，其中血管信号强度的动态范围为 20 dB。最后，每个图像都被放大了不同的值，以获得涉及广泛 SNR 范围的测量。 Table I 总结了所有参数。
Simulation Results 模拟结果。
We tested the reconstruction methods using both simulation and experimental data. We compared CNN-based approaches with other common reconstruction methods including DAS, DMAS and/or CS. Here, we call a network using DAS results (without Hilbert transform) a single-channel input ‘UNET’ and a network using 3-D transformed arrays a multi-channel input ‘upgUNET’. As shown in Fig. 5, UNET learns in the image-domain (2D projection) while upgUNET learns in the data-domain (3D array). DAS employs a rectangular apodization function where the activated aperture size is determined by f-number (=0.5). For image display, it uses the Hilbert transform following summation. The iterative method exploits total-variation (TV) regularization and wavelet transforms for sparsity dictionaries. We used the general criteria that the iteration stops when the norm of the gradient is less than a threshold value. We checked that more iterations barely change image quality or metric values. We empirically selected all parameter values minimizing the error between ground-truth and reconstruction results. Table II summarizes all parameters for the reconstructions.
我们使用模拟和实验数据测试了重建方法。我们将基于 CNN 的方法与其他常见的重建方法（包括 DAS、DMAS 和/或 CS）进行了比较。这里，我们称使用 DAS 结果（不使用希尔伯特变换）的网络为单通道输入的“UNET”，使用 3-D 变换数组的网络为多通道输入的“upgUNET”。如 Fig. 5 所示，UNET 在图像域（2D 投影）中学习，而 upgUNET 在数据域（3D 数组）中学习。DAS 采用矩形权重函数，激活的孔径大小由 f-number（=0.5）决定。在图像显示中，它在求和后使用希尔伯特变换。迭代方法利用总变差（TV）正则化和小波变换进行稀疏字典。我们使用的一般标准是，当梯度的范数小于阈值时，迭代停止。我们检查了更多的迭代几乎不会改变图像质量或度量值。我们根据最小化地面真实数据和重建结果之间的误差，经验性地选择了所有参数值。 Table II 总结了重建的所有参数。
表 II 重建方法的参数
The same procedure was used to generate testing and training datasets. However, two sets were generated from independent objects for independent verification. Fig. 6 shows imaging results using the selected reconstruction methods from two particular examples where object shapes and data SNR are totally different. As expected, standard DAS (f-number =0.5) followed by Hilbert transformation provides low-contrast, poor-resolution images. While MV, DMAS and CS improve contrast in general, they often suppress weak signals. CNN-based methods restore most of the vasculature with stronger contrast and higher resolution. For upgUNET processing, fine vessels are more clearly visible, as shown in the circled areas of Figs. 6 (M) and (N). Lost structures are mostly vertically-extended vessels because their signal power is extremely low.
生成测试和训练数据集使用了相同的程序。然而，两组数据是从独立的对象中生成，用于独立验证。 Fig. 6 展示了使用选定的重建方法从两个特定示例中得到的成像结果，其中对象形状和数据信噪比完全不同。如预期，标准的 DAS（f 数=0.5）接着希尔伯特变换提供了低对比度、分辨率差的图像。虽然 MV、DMAS 和 CS 通常能改善对比度，但它们经常抑制弱信号。基于 CNN 的方法恢复了大部分血管结构，具有更强的对比度和更高的分辨率。对于 upgUNET 处理，如 Figs. 6 (M) and (N) 中圈出的区域所示，细小血管更清晰可见。丢失的结构主要是垂直延伸的血管，因为它们的信号功率极低。
Reconstruction results using synthetic data. Two particular objects are tested and all images are displayed using a log-scale colormap. (A, B) Ground-truth images. (C, D) Delay-and-sum results. Hilbert transform is applied as post-processing. (E, F) Minimum variance results. (G, H) Delay-multiply-and-sum results. (I, J) Iterative CS method results. Wavelet dictionaries and total-variation regularization are used for compressed sensing. (K, L) Deep-learning results. An input is a 2-D array using DAS. (M, N) Deep-learning results. An input is a 3-D multi-channel array.
使用合成数据的重建结果。测试了两个特定的对象，所有图像都使用对数尺度颜色图显示。（A, B）真实图像。（C, D）延迟求和结果。希尔伯特变换作为后处理应用。（E, F）最小方差结果。（G, H）延迟乘积求和结果。（I, J）迭代 CS 方法结果。压缩感知使用了小波字典和总变分正则化。（K, L）深度学习结果。输入是一个使用 DAS 的 2-D 数组。（M, N）深度学习结果。输入是一个 3-D 多通道数组。
In addition to these qualitative comparisons, we quantified performance differences employing the peak-signal-to-noise ratio (PSNR) and structure similarity (SSIM) metrics. The PSNR is defined via the MSE as
\begin{equation*} \eta _{PSNR} = 10\log _{10} \left({\frac {n_{1} n_{2} I_{max}^{2}}{ \vert \vert \mathbf {s}-\bar { \mathbf {s}}\vert \vert _{F}^{2}} }\right)\tag{17}\end{equation*}
where \begin{equation*} \eta _{SSIM} = \frac {(2~\mu _{ \mathbf {s}}\mu _{\bar { \mathbf {s}}}+c_{1})(2\sigma _{ \mathbf {s},\bar { \mathbf {s}}}+c_{2}) }{(\mu _{ \mathbf {s}}^{2} + \mu _{\bar { \mathbf {s}}}^{2} + c_{1})(\sigma _{ \mathbf {s}}^{2} + \sigma _{\bar { \mathbf {s}}}^{2} + c_{2}) }\tag{18}\end{equation*}
除了这些定性比较外，我们还使用峰值信噪比（PSNR）和结构相似性（SSIM）指标来量化性能差异。PSNR 通过 MSE 定义为
\begin{equation*} \eta _{PSNR} = 10\log _{10} \left({\frac {n_{1} n_{2} I_{max}^{2}}{ \vert \vert \mathbf {s}-\bar { \mathbf {s}}\vert \vert _{F}^{2}} }\right)\tag{17}\end{equation*}
\begin{equation*} \eta _{SSIM} = \frac {(2~\mu _{ \mathbf {s}}\mu _{\bar { \mathbf {s}}}+c_{1})(2\sigma _{ \mathbf {s},\bar { \mathbf {s}}}+c_{2}) }{(\mu _{ \mathbf {s}}^{2} + \mu _{\bar { \mathbf {s}}}^{2} + c_{1})(\sigma _{ \mathbf {s}}^{2} + \sigma _{\bar { \mathbf {s}}}^{2} + c_{2}) }\tag{18}\end{equation*}
Table III presents average PSNR and SSIM values computed from 1,000 samples where every vascular shape is different. The numbers in parentheses are the standard deviations. Deep-learning approaches offer significant gain over standard methods, strongly suggesting that the network provides quantitatively better image quality. Note that the upgUNET produces the best values, in agreement with visual inspection.
Table III 展示了从 1,000 个样本中计算得出的平均 PSNR 和 SSIM 值，每个血管形状都不同。括号中的数字是标准差。深度学习方法相比标准方法有显著提升，强烈表明网络提供了更好的图像质量。注意，upgUNET 产生了最佳值，这与视觉检查一致。
表 III 不同方法的定量比较
Experimental Results 实验结果
A. PAUS System A. PAUS 系统
Our customized system for spectroscopic PA imaging is illustrated in Fig. 7. The scanner (Vantage, Verasonics, WA, USA) is programmed to record RF data at different wavelengths and fiber positions. A compact diode-pumped laser (TiSon GSA, Laser-export, Russia) transmits a light pulse at any arbitrary wavelength ranging from 700 to 900 nm at a 1 kHz rate. Transmitted pulses are sequentially delivered to 20 fiber terminals mounted around the top and bottom surfaces of a linear array transducer (LA 15/128-1633, Vermon S.A. France). The transducer center frequency is 15 MHz and the 3 dB bandwidth is around 8 MHz. US firings are interspersed with laser firings such that a full PA/US image frame at a fixed optical wavelength is recorded every 20 ms, producing a 50 Hz display rate for integrated US and PA images. System details are described in [1].
我们为光谱 PA 成像定制的系统如 Fig. 7 所示。扫描仪（Vantage，Verasonics，WA，美国）被编程以在不同的波长和光纤位置记录 RF 数据。一种紧凑的二极管泵浦激光器（TiSon GSA，Laser-export，俄罗斯）以 1 kHz 的速率在 700 到 900 nm 的任意波长传输光脉冲。传输的脉冲依次传送到围绕线性阵列换能器（LA 15/128-1633，Vermon S.A. 法国）顶部和底部表面安装的 20 个光纤端子。换能器中心频率为 15 MHz，3 dB 带宽约为 8 MHz。US次发射与激光发射交错，以便每 20 ms 记录一个固定光学波长的完整 PA/US图像帧，为集成US和 PA 图像产生 50 Hz 的显示率。系统细节在 [1] 中描述。
Our customized PAUS system. An US scanner triggers a compact diode-pumped laser such that it emits pulses (around 1 mJ energy) at a 1 kHz rate with switching wavelength ranging from 700 nm to 900 nm. The laser is delivered to integrated fibers arranged on the two sides of a linear array transducer. A motor controlled by the scanner allows laser pulses to couple with different fibers sequentially. The scanner records PA signals that originate from light propagation into tissue.
我们定制的PAUS系统。一个US扫描器触发一个紧凑的二极管泵浦激光器，使其以1 kHz的频率发射脉冲（能量约为1 mJ），波长切换范围从700 nm到900 nm。激光被输送到排列在线性阵列换能器两侧的集成光纤中。扫描器控制的电机允许激光脉冲依次与不同的光纤耦合。扫描器记录由光在组织中传播产生的PA信号。
Here, for the purpose of reconstruction tests, we acquired data using one wavelength at 795 nanometers. The sampling rate for acoustic array data
在这里，为了重建测试的目的，我们使用 795 纳米的一个波长获取数据。声学阵列数据
B. Phantom Study B. 幻影研究
We constructed a phantom containing a metal wire acting as an optical absorption target. As shown in Fig. 8, the wire shape approximates the letter ‘W’. It was suspended from a cubical container such that it appears as the ‘W’ shape in the z-x imaging plane. The container was filled with an intralipid solution (Fresenius Kabi, Deerfield, USA) acting as a scattering medium. The concentration of the intralipid is around 2% and the effective attenuation coefficient is around 0.1 mm−1. Channel data were recorded with our customized system. Reconstruction results are shown in Fig. 8, which compares deep-learning methods with standard methods. For DAS, we tested a small f-number (0.1) in addition to the default value (0.5). The lower number corresponds to larger aperture size, which can access some information from diagonal lines at the expense of SNR in the entire field. In agreement with simulation results, DMAS and CS improve contrast but suppress weak structures.
我们构建了一个包含金属线的幻影，该金属线作为光吸收目标。如 Fig. 8 所示，线形近似于字母‘W’。它被悬挂在一个立方体容器中，使其在 z-x 成像平面上呈现为‘W’形状。容器内充满了作为散射介质的 Intralipid 溶液（Fresenius Kabi，Deerfield，美国）。Intralipid 的浓度约为 2%，有效衰减系数约为 0.1 mm −1 。通道数据由我们的定制系统记录。重建结果显示在 Fig. 8 中，比较了深度学习方法和标准方法。对于 DAS，我们测试了一个小的 f 数（0.1）以及默认值（0.5）。较低的数字对应于更大的孔径大小，可以访问对角线上的一些信息，但以整个场的信噪比为代价。与仿真结果一致，DMAS 和 CS 提高了对比度但抑制了弱结构。
Reconstruction results. All images are displayed using a log-scale colormap (35 dB range). A ‘W’ shape wire is scanned by the PAUS system. (A) Delay-and-sum result. The f-number is 0.5. (B) Delay-and-sum result. The f-number is 0.1. (C) Delay-multiply-and-sum results. (D) Iterative CS method results. (E) UNET deep-learning result. An input is a 2-D array using DAS. (F) upgUNET deep-learning result. An input is a 3-D multi-channel tensor.
重建结果。所有图像均使用对数刻度颜色映射（35 dB 范围）显示。一个“W”形的线被 PAUS 系统扫描。（A）延迟求和结果。F 数为 0.5。（B）延迟求和结果。F 数为 0.1。（C）延迟乘求和结果。（D）迭代 CS 方法结果。（E）UNET 深度学习结果。输入是一个使用 DAS 的二维数组。（F）upgUNET 深度学习结果。输入是一个三维多通道张量。
In the deep-learning imaging results, object shapes are more distinct with higher resolution. In particular, the preferred upgUNET method restores most wire structure. One flaw in the deep-learning approaches is that the networks sometimes produce artifacts near objects, as shown in Figs. 8 (E) and (F). We believe they arise from low-level reverberations not modeled in synthetic training data.
在深度学习成像结果中，对象形状更加清晰，分辨率更高。特别是，首选的 upgUNET 方法恢复了大部分线结构。深度学习方法的一个缺陷是，网络有时会在对象附近产生伪影，如 Figs. 8 (E) and (F) 所示。我们认为这些伪影是由于合成训练数据中未建模的低级反射引起的。
C. In-Vivo Test C. 体内测试
Lastly, we scanned and imaged a human finger to study the feasibility of our suggested method for in-vivo vascular imaging in real-time. These studies were approved by the Institutional Review Board (IRB) of the University of Washington (Study# 00009196) and used both optical and acoustic energies well within ANSI (optical) and FDA (US) guidelines. The PAUS system recorded US and PA measurements using an interleaved pulse sequence for 50 Hz frame-rate imaging.
最后，我们扫描并成像了一个人类手指，以研究我们建议的方法在实时体内血管成像中的可行性。这些研究已获得华盛顿大学机构审查委员会（IRB）的批准（研究编号 00009196），并在 ANSI（光学）和 FDA（US）指南范围内使用了光学和声学能量。PAUS 系统记录了 US 和 PA 测量，使用交错脉冲序列进行 50 Hz 帧率成像。
Figs. 9 (A) and (B) show two longitudinal cross-sections of the finger as US B-mode images. We applied DAS, UNET and upgUNET to PA reconstructions for visual comparison (See Fig. 9 (C-H)). In this limited test, we found little difference between UNET and upgUNET reconstructions. We presume that upgUNET presents more realistic microvascular structures than B-mode images. However, it is obvious that deep learning images provide markedly superior contrast and resolution than equivalent images reconstructed with DAS.
Figs. 9 (A) and (B) 显示了手指的两个纵向横截面作为 US B 模式图像。我们对 PA 重建应用了 DAS、UNET 和 upgUNET 进行视觉比较（见 Fig. 9 (C-H) ）。在这个有限的测试中，我们发现 UNET 和 upgUNET 重建之间几乎没有差异。我们推测 upgUNET 比 B 模式图像呈现更真实的微血管结构。然而，很明显，深度学习图像提供的对比度和分辨率远优于使用 DAS 重建的等效图像。
In vivo reconstruction results. A human finger is scanned by the PAUS system. Two sagittal planes are tested. (A, B) US B-mode image. (C, D) Delay-and-sum result. The f-number is 0.5. (E, F) UNET deep-learning result. An input is a 2-D array using DAS. (G, H) upgUNET deep-learning result. An input is a 3-D multi-channel tensor. All images are displayed using a log-scale colormap. Mapping ranges for US and PD images are 50dB and 40dB, respectively.
体内重建结果。一个人的手指被 PAUS 系统扫描。测试了两个矢状面。（A, B）US B 模式图像。（C, D）延迟求和结果。F 数为 0.5。（E, F）UNET 深度学习结果。输入是一个使用 DAS 的二维数组。（G, H）upgUNET 深度学习结果。输入是一个三维多通道张量。所有图像均使用对数比例颜色图显示。US 和 PD 图像的映射范围分别为 50dB 和 40dB。
Discussion 讨论
As discussed in the introduction, the inverse problem of PA imaging can be solved exactly only when the detection surface represents a whole sphere, cylinder or infinite plane. These conditions can be obtained for small animal imaging [42] but are difficult to achieve in a clinical environment.
如引言中所讨论的，PA 成像的逆问题只有当检测表面代表一个完整的球体、圆柱体或无限平面时才能精确解决。这些条件可以在小动物成像 [42] 中获得，但在临床环境中难以实现。
Hand-held US probes with relatively narrow bandwidth yields serious image quality losses for PA images, and such techniques alone will unlikely be accepted for medical use [43].
手持US探头的相对狭窄带宽会严重损失 PA 图像的质量，仅凭这种技术不太可能被医疗使用接受 [43] 。
However interleaved spectroscopic PAUS has been recently shown to dramatically improve the capabilities of diagnostic US in monitoring interventional procedures even under limited view and bandwidth conditions [3]. Very recently, a fast-sweep PAUS approach has been developed to operate at clinically acceptable frame rates for both PA and US modalities [1].
然而，最近交错的光谱 PAUS 已被证明可以显著提高在有限视图和带宽条件下监测介入程序的诊断US能力。最近，一种快速扫描 PAUS 方法已被开发出来，在 PA 和US模式下以临床可接受的帧率运行 [1] 。
The goal of this paper was to investigate whether deep-learning algorithms can improve the quality of images obtained with a hand-held US probe. We note that large-scale absorbing heterogeneities are unlikely to be fully corrected with the proposed upgUNET algorithm because the low frequencies associated with these objects are not preserved within the limited bandwidth of the detector. In our opinion, however, microvascular networks can be improved with advanced reconstruction algorithms based on upgNET. Although not fully validated over a wide range of experimental conditions, the simulations, phantom measurements and in vivo results presented here strongly support this statement.
本文的目的是调查深度学习算法是否可以提高使用手持US探头获得的图像质量。我们注意到，大规模吸收异质性不太可能通过所提出的 upgUNET 算法完全纠正，因为这些对象的低频率在探测器的有限带宽内没有被保留。然而，在我们看来，微血管网络可以通过基于 upgNET 的高级重建算法得到改善。尽管没有在广泛的实验条件下得到充分验证，但这里呈现的模拟、幻影测量和体内结果强烈支持这一说法。
We have explored several different signal processing methods beyond traditional DAS reconstruction and found that deep learning has the potential for enhanced image quality and real-time implementation if input data are structured to match the network architecture of a reasonably sized net. Here we reformatted elemental signals from an US transducer array into multi-channel 3D data (tensor) using prior knowledge of propagation delays related to simple wave physics. Effective decomposition of this tensor by the neural network can capture latent structures and features.
我们探索了多种不同的信号处理方法，超越了传统的 DAS 重建，并发现如果输入数据结构与合理大小网络的网络架构相匹配，深度学习有可能提高图像质量并实现实时实施。在这里，我们使用了与简单波动物理相关的传播延迟的先验知识，将来自US换能器阵列的元素信号重新格式化为多通道 3D 数据（张量）。神经网络对这个张量的有效分解可以捕捉潜在的结构和特征。
Both experiments and simulations have demonstrated how the proposed neural network improves image quality for PA reconstruction. It produces images with stronger contrast, higher spatial resolution compared to DAS, and few structure losses even given the limited spatial and temporal bandwidth of the real-time system used for PA data acquisition. Additional image quality improvements were demonstrated using multi-channel data as the network input (upgUNET). The final advantage of this approach is that the network effectively learns filtering weights from training data while standard methods must explicitly impose parameter values, filtering dictionaries, or regularizers.
实验和模拟都证明了所提出的神经网络如何提高 PA 重建的图像质量。与 DAS 相比，它产生的图像对比度更强，空间分辨率更高，并且即使在实时系统用于 PA 数据采集的空间和时间带宽有限的情况下，结构损失也很少。使用多通道数据作为网络输入（upgUNET）进一步展示了图像质量的改进。这种方法的最终优势在于网络能够从训练数据中有效学习过滤权重，而标准方法必须显式施加参数值、过滤字典或正则化器。
Again, the novelty of our architecture is leveraging image reconstruction in the data-domain. It is less computational complex than end-to-end fully-connected networks that directly learn from raw-data. In addition, it preserves information better than image-to-image networks. Our architecture is based on U-net to efficiently cope with the non-stationary nature of off-axis signals, artifacts and noise, as well as to restore vascular structures. Compared to DL-based iterative schemes, it is much better suited to real-time imaging due to lower computations.
再次强调，我们架构的新颖之处在于利用数据域中的图像重建。它的计算复杂性低于直接从原始数据学习的端到端全连接网络。此外，它比图像到图像的网络更好地保留了信息。我们的架构基于 U-net，以有效应对非定常性的离轴信号、伪影和噪声，以及恢复血管结构。与基于 DL 的迭代方案相比，由于计算量较低，它更适合实时成像。
The computational cost of the matrix multiplication mapping a 2D data array (
将二维数据数组（
The average computation time for image reconstruction from raw data is 40 msec, representing a 25 Hz real-time frame rate. Our current system functions at 50 Hz and higher frame rates, so these computation times must be reduced by at least a factor of two for true real-time implementation. Given an optimized hardware architecture for our specific reconstruction approach, a real-time frame rate of 50 Hz and higher is very realistic in the short term for the specific deep learning algorithm presented here for PA image reconstruction.
从原始数据重建图像的平均计算时间为 40 毫秒，代表了 25Hz 的实时帧率。我们当前的系统功能在 50Hz 及更高帧率下运行，因此这些计算时间必须至少减少两倍以实现真正的实时实现。考虑到我们特定重建方法的优化硬件架构，短期内实现 50Hz 及更高的实时帧率对于这里呈现的特定深度学习算法来说是非常现实的。
The architecture combining a trainable neural network with a non-trainable preprocessor inputting raw data can be effective for diverse medical imaging fields and remote sensing applications. For example, some groups have recently investigated trainable apodization methods to replace slow MV beamforming for US B-mode imaging [45], [46]. They showed that a simple network enables fast imaging computation without reducing MV image quality. In addition, one group has applied the similar architecture to a semi-cylindrical PA tomographic system [47].
将可训练神经网络与非可训练预处理器结合的架构，输入原始数据，可有效应用于多种医学成像领域和远程感测应用。例如，一些团队最近研究了可训练的调制方法，以替代速度较慢的MV波束成形技术进行B模式成像。他们展示了一个简单的网络能够实现快速成像计算，而不降低MV图像质量。此外，有一个团队将类似的架构应用于半圆柱形PA断层扫描系统。
As noted in Sections IV and V, the deep learning approach presented here was trained over a wide range of SNR and demonstrated superior results compared to conventional DAS at all tested levels of SNR. This characteristic is especially attractive for PAUS imaging in vivo using a handheld probe with limited spatial and temporal bandwidth where depth-dependent light attenuation limits the penetration depth of PA images. To demonstrate potential gains in penetration depth with DL-based reconstruction, a test case (Fig. 10) was analyzed where a simple model of fluence variation with depth was applied to ground truth data.
如 Sections IV 和 V 中所述，这里介绍的深度学习方法经过了广泛的信噪比训练，并在所有测试的信噪比水平上显示出比传统 DAS 更优越的结果。这一特性对于使用手持探头进行的 PAUS 成像尤为吸引人，因为手持探头的空间和时间带宽有限，且深度依赖的光衰减限制了 PA 图像的穿透深度。为了展示基于 DL 重建的穿透深度潜在增益，分析了一个测试案例（ Fig. 10 ），其中应用了一个简单的随深度变化的光流模型到基础真实数据上。
![Fig. 10. - Reconstruction results using synthetic data to explore penetration depth when PA signal attenuates with depth. One particular absorption object is tested and all images are displayed using a log-scale colormap. (A) Sum of light fluences from 20 fibers in a scattering medium. MCX Studio is used for the Monte Carlo simulation [44]. The details of the transducer geometry can be found Ref [1]. The effective attenuation coefficient of the medium is 1.73 cm−1. (B) Ground-truth object. (C) Attenuated object due to the heterogeneous fluence. Measurements are obtained from the object using Eq. 3. (D-G) Reconstruction results using (D) DAS, (E) deep-learning, (F) DAS with fluence compensation, and (G) deep-learning with fluence compensation.](/mediastore_new/IEEE/content/media/42/9242349/9091172/kim10abcdefg-2993835-large.gif)
Reconstruction results using synthetic data to explore penetration depth when PA signal attenuates with depth. One particular absorption object is tested and all images are displayed using a log-scale colormap. (A) Sum of light fluences from 20 fibers in a scattering medium. MCX Studio is used for the Monte Carlo simulation [44]. The details of the transducer geometry can be found Ref [1]. The effective attenuation coefficient of the medium is 1.73 cm−1. (B) Ground-truth object. (C) Attenuated object due to the heterogeneous fluence. Measurements are obtained from the object using Eq. 3. (D-G) Reconstruction results using (D) DAS, (E) deep-learning, (F) DAS with fluence compensation, and (G) deep-learning with fluence compensation.
使用合成数据重建结果以探索 PA 信号随深度衰减时的穿透深度。测试了一个特定的吸收对象，所有图像都使用对数比例颜色图显示。（A）来自 20 个纤维在散射介质中的光流量总和。MCX Studio 用于蒙特卡洛模拟 [44] 。换能器几何的细节可以在参考文献 [1] 中找到。介质的有效衰减系数为 1.73 厘米 −1 。（B）基准对象。（C）由于异质流量而衰减的对象。从对象获取测量数据 Eq. 3 。（D-G）使用（D）DAS，（E）深度学习，（F）带流量补偿的 DAS，和（G）带流量补偿的深度学习的重建结果。
The sum of fluences from 20 fibers is presented in Fig. 10 (A) for a scattering medium with an effective attenuation coefficient of 1.73 cm−1. MCX Studio was used for the Monte Carlo simulation [44] producing this distribution. The PA signal for this case (Fig. 10 (C)) depends not only on the absorption of the vascular object (Fig. 10 (B)) but also on the depth-dependent fluence. Figs. 10 (D) and (E) compare conventional DAS and DL-based reconstructions with noise added. The average photoacoustic SNR for this example is about 25 dB at a depth of 5 mm and about 5 dB at a depth of 15 mm.
20 个纤维的光流量总和在 Fig. 10 (A) 中呈现，用于具有 1.73 厘米有效衰减系数 −1 的散射介质。MCX Studio 用于进行蒙特卡洛模拟 [44] ，产生这种分布。这种情况下的 PA 信号（ Fig. 10 (C) ）不仅取决于血管对象的吸收（ Fig. 10 (B) ），还取决于深度依赖的流量。 Figs. 10 (D) and (E) 比较了传统的 DAS 和基于 DL 的重建，并添加了噪声。这个例子的平均光声信噪比在 5 毫米深度约为 25 dB，在 15 毫米深度约为 5 dB。
Figs. 10 (F) and (G) compare DAS and DL-based reconstruction after fluence compensation has been applied. The specific compensation approach was similar to that applied in many PA systems in which the reconstructed images in Figs. 10 (D) and (E) are multiplied by a depth-dependent gain factor, similar to the time-gain control (TGC) function in a real-time US scanner, based on the expected fluence attenuation function. Clearly, DL-based reconstruction not only improves overall image quality as noted above, but also increases PA image penetration by about 50% for this particular case. In general, enhanced penetration for real-time, handheld PAUS imaging systems is another attractive feature of DL-based reconstruction.
Figs. 10 (F) and (G) 在应用了流量补偿后，比较 DAS 和基于 DL 的重建。具体的补偿方法与许多 PA 系统中应用的方法类似，在 Figs. 10 (D) and (E) 中重建的图像乘以一个深度依赖的增益因子，类似于实时US扫描仪中的时间增益控制（TGC）功能，基于预期的流量衰减函数。显然，基于 DL 的重建不仅如上所述改善了整体图像质量，而且还将 PA 图像的穿透力提高了约 50%。总的来说，对于实时手持 PAUS 成像系统，增强的穿透力是基于 DL 的重建的另一个吸引人的特性。
One limitation of the CNN identified in these studies is poor image quality for structures aligned almost vertically (i.e., nearly parallel to the normal to the 1D transducer array). As expected, restoration is challenging for this case since most of the PA signal spatial frequencies are not captured by the transducer array. As shown in the top images of Fig. 6, this loss is more serious if the PA signal power of the object is weak compared with the noise power. Part of our future work will explore alternative approaches leveraging additional information from ultrasound imaging, or better system conditions such as acquiring one additional view at a significant angle with respect to the array normal at the first position.
这些研究中发现的 CNN 的一个限制是对几乎垂直（即，几乎平行于 1D 换能器阵列法线）排列的结构的图像质量较差。如预期，由于大部分 PA 信号的空间频率未被换能器阵列捕获，这种情况的恢复是具有挑战性的。如 Fig. 6 顶部图像所示，如果物体的 PA 信号功率与噪声功率相比较弱，这种损失更为严重。我们未来的工作部分将探索利用来自超声成像的额外信息，或更好的系统条件，例如在第一个位置相对于阵列法线获取一个显著角度的额外视图的替代方法。
The second limitation of the method is the presence of some unexpected image artifacts for real data tests, which can reduce specificity. Obtaining real ground-truth maps is impracticable at large scale. Thus, reducing discrepancies between synthetic and real data is needed to guarantee that a trained CNN works best for real data.
该方法的第二个限制是在实际数据测试中出现了一些意外的图像伪影，这可能降低特异性。在大规模上获取真实的基准图是不切实际的。因此，需要减少合成数据和实际数据之间的差异，以保证训练有素的 CNN 对实际数据最有效。
One of the main difficulties in reconstruction arises from the object dimension. Although a transducer lens focuses to a 2D plane (imaging plane), it cannot fully limit the sensitivity to a selected plane in an object. In other words, signals still come from points outside the imaging plane. Thus, generating training data based on 3D structures is required for the next step. Experiments are ongoing to map the full 3D PSF of our PA system and to include these details into the forward model for synthetic data generation. We expect that new synthetic data accompanied with an extended detection view will further improve image quality.
重建中的一个主要困难来自于对象的尺寸。尽管换能器透镜聚焦于一个二维平面（成像平面），但它不能完全限制对选定平面的敏感度。换句话说，信号仍然来自成像平面之外的点。因此，基于三维结构生成训练数据是下一步所需的。我们正在进行实验，以绘制我们的 PA 系统的完整三维 PSF，并将这些细节纳入合成数据生成的前向模型中。我们期望新的合成数据配合扩展的检测视图将进一步提高图像质量。
There certainly are alternative learning frameworks. Our architecture is limited to supervised learning requiring ground-truth pairing with input data. Some architectures such as Generate Adversarial Networks (GANs) do not need pairs [48], [49]. They can provide realistic fake images from measurements using real vascular image datasets obtained by other medical imaging modalities. Since only real data and real images are required for training, they can prevent unexpected artifacts. Also, different learning frameworks can bring additional information from US images to bear on PA image reconstruction.
当然有其他的学习框架。我们的架构仅限于需要与输入数据配对的监督学习。一些架构，如生成对抗网络（GANs），不需要配对 [48] ， [49] 。它们可以使用通过其他医学成像方式获得的真实血管图像数据集，从测量中提供逼真的假图像。由于训练只需要真实数据和真实图像，它们可以防止意外的伪影。此外，不同的学习框架可以从US图像中带来额外的信息，以用于 PA 图像重建。
Overall, our future studies will address these two limitations and develop specific deep learning architectures for real-time implementation. The goal is to create a CNN tuned to the problem of PA image reconstruction using limited spatial and temporal bandwidth data for robust, high-quality spectroscopic PAUS imaging at frame rates of 50 Hz and higher.
总体而言，我们未来的研究将针对这两个局限性，并开发专门的深度学习架构以实现实时应用。我们的目标是创建一个针对PA图像重建问题调整过的卷积神经网络（CNN），使用有限的空间和时间带宽数据，以实现稳健、高质量的光谱光声超声成像（PAUS），帧率达到50 Hz及以上。
Conclusions 结论
In this paper, we described a deep convolutional network for real-time PAUS imaging. The PAUS platform has the potential for real-time clinical implementation but the limited view/bandwidth of clinical US arrays degrades PA image quality. Our target of interest is microvessels, which can be almost entirely restored by advanced reconstruction methods, in contrast to large homogeneous objects. We developed a deep learning network for this application trained with realistic simulation data synthesized from ground-truth vascular image sets using a deterministic PA forward operator and the measured impulse response of a real transducer. Reformatting raw channel data into a multi-channel array as a pre-processing step increased learning efficiency with respect to network complexity and imaging performance. The neural network is based on U-net, decomposing signals via multi-scale feature maps. Coupling the trainable network with the transformation method, we imaged structures mimicking vascular networks both in simulation and experiment. Overall, this approach reconstructs PA data with much higher image quality than conventional methods but loses some portions of complex absorber geometries and generates minor artifacts.
在本文中，我们描述了一个用于实时 PAUS 成像的深度卷积网络。PAUS 平台有实时临床实施的潜力，但临床US阵列的有限视角/带宽降低了 PA 图像质量。我们感兴趣的目标是微血管，与大型均质物体相比，可以通过先进的重建方法几乎完全恢复。我们为这一应用开发了一个深度学习网络，该网络使用从真实血管图像集合合成的现实模拟数据进行训练，使用确定性 PA 前向算子和真实换能器的测量脉冲响应。将原始通道数据重新格式化为多通道阵列作为预处理步骤，提高了学习效率，降低了网络复杂性和提高了成像性能。该神经网络基于 U-net，通过多尺度特征图解构信号。将可训练网络与转换方法结合，我们在模拟和实验中成像模拟血管网络的结构。总体而言，这种方法重建的 PA 数据比传统方法具有更高的图像质量，但丢失了一些复杂吸收体几何形状的部分并产生了轻微的伪影。
Appendix I 附录一Adjoint of Operator F 操作符 F 的伴随
Adjoint of Operator F 操作符 F 的伴随
The solution of the forward model in Eq. 2 can simply be written as
\begin{equation*} p(t, \mathbf {r}') = F(s(\mathbf {r})),\tag{19}\end{equation*}
where \begin{align*} < &F(s(\mathbf {r})), q(t, \mathbf {r}') > \\=&\int _\Xi \int _{R} \left[{ \beta \int _\Omega \frac {\frac {\partial }{\partial t}\delta \left({t-\frac {\vert \mathbf {r} - \mathbf {r}'\vert }{v_{s}}}\right)}{\vert \mathbf {r} - \mathbf {r}'\vert } s (\mathbf {r}) d \mathbf {r}}\right] q (t, \mathbf {r}') d \mathbf {r}'dt \\=&\int _\Omega s (\mathbf {r}) \left[{ \int _\Xi \frac {-\beta }{\vert \mathbf {r} - \mathbf {r}'\vert } \tilde {q} \left({\frac {\vert \mathbf {r} - \mathbf {r}'\vert }{v_{s}}, \mathbf {r}'}\right) d \mathbf {r}' }\right] d \mathbf {r} \\=& < s, F^{T}(\tilde {q})>\tag{20}\end{align*}
正向模型在 Eq. 2 中的解可以简单地写为
\begin{equation*} p(t, \mathbf {r}') = F(s(\mathbf {r})),\tag{19}\end{equation*}
\begin{align*} < &F(s(\mathbf {r})), q(t, \mathbf {r}') > \\=&\int _\Xi \int _{R} \left[{ \beta \int _\Omega \frac {\frac {\partial }{\partial t}\delta \left({t-\frac {\vert \mathbf {r} - \mathbf {r}'\vert }{v_{s}}}\right)}{\vert \mathbf {r} - \mathbf {r}'\vert } s (\mathbf {r}) d \mathbf {r}}\right] q (t, \mathbf {r}') d \mathbf {r}'dt \\=&\int _\Omega s (\mathbf {r}) \left[{ \int _\Xi \frac {-\beta }{\vert \mathbf {r} - \mathbf {r}'\vert } \tilde {q} \left({\frac {\vert \mathbf {r} - \mathbf {r}'\vert }{v_{s}}, \mathbf {r}'}\right) d \mathbf {r}' }\right] d \mathbf {r} \\=& < s, F^{T}(\tilde {q})>\tag{20}\end{align*}
Appendix II 附录二Framelet Expansion 小波框架展开
Framelet Expansion 小波框架展开
Preprocessed data can be represented as a third-order tensor
\begin{align*} H_{d_{1},d_{2}}(\mathbf {F})=&\begin{bmatrix} H_{d_{1},d_{2}}(\mathbf {F}_{1}) &\quad \cdots &\quad H_{d_{1},d_{2}}(\mathbf {F}_{J}) \end{bmatrix}, \\ H_{d_{1},d_{2}}(\mathbf {F}_{j})=&\begin{bmatrix} H_{d_{1}}(\mathbf {f}_{1,j}) &\quad H_{d_{1}}(\mathbf {f}_{2,j}) &\quad \cdots &\quad H_{d_{1}}(\mathbf {f}_{d_{2},j}) \\ H_{d_{1}}(\mathbf {f}_{2,j}) &\quad H_{d_{1}}(\mathbf {f}_{3,j}) &\quad \cdots &\quad H_{d_{1}}(\mathbf {f}_{d_{2}+1,j}) \\ \vdots &\quad \vdots &\quad \ddots &\quad \vdots \\ H_{d_{1}}(\mathbf {f}_{M,j}) &\quad H_{d_{1}}(\mathbf {f}_{1,j}) &\quad \cdots &\quad H_{d_{1}}(\mathbf {f}_{d_{2}-1,j}) \end{bmatrix}, \!\! \\\tag{21}\end{align*}
where \begin{align*} H_{d_{1}}(\mathbf {f}_{m,j}) \!=\! \begin{bmatrix} \mathbf {f}_{m,j}[{1}] &\quad \mathbf {f} _{m,j}[{2}] &\quad \cdots &\quad \mathbf {f} _{m,j}[d_{1}] \\ \mathbf {f}_{m,j}[{2}] &\quad \mathbf {f} _{m,j}[{3}] &\quad \cdots &\quad \mathbf {f} _{m,j}[d_{1}+1] \\ \vdots &\quad \vdots &\quad \ddots &\quad \vdots \\ \mathbf {f}_{m,j}[N] &\quad \mathbf {f} _{m,j}[{1}] &\quad \cdots &\quad \mathbf {f} _{m,j}[d_{1}-1] \end{bmatrix},\!\! \\\tag{22}\end{align*}
\begin{equation*} H_{d_{1} d_{2}} (\mathbf {F}) = \frac {1}{\alpha } \tilde {\mathbf {\Phi }}\mathbf {\Phi }^{T} H_{d_{1} d_{2}} (\mathbf {F}) {\mathbf {\Psi }}\tilde {\mathbf {\Psi }}^{T} = \frac {1}{\alpha } \tilde {\mathbf {\Phi }}\mathbf {C}\tilde {\mathbf {\Psi }}^{T}\tag{23}\end{equation*}
预处理数据可以表示为三阶张量
\begin{align*} H_{d_{1},d_{2}}(\mathbf {F})=&\begin{bmatrix} H_{d_{1},d_{2}}(\mathbf {F}_{1}) &\quad \cdots &\quad H_{d_{1},d_{2}}(\mathbf {F}_{J}) \end{bmatrix}, \\ H_{d_{1},d_{2}}(\mathbf {F}_{j})=&\begin{bmatrix} H_{d_{1}}(\mathbf {f}_{1,j}) &\quad H_{d_{1}}(\mathbf {f}_{2,j}) &\quad \cdots &\quad H_{d_{1}}(\mathbf {f}_{d_{2},j}) \\ H_{d_{1}}(\mathbf {f}_{2,j}) &\quad H_{d_{1}}(\mathbf {f}_{3,j}) &\quad \cdots &\quad H_{d_{1}}(\mathbf {f}_{d_{2}+1,j}) \\ \vdots &\quad \vdots &\quad \ddots &\quad \vdots \\ H_{d_{1}}(\mathbf {f}_{M,j}) &\quad H_{d_{1}}(\mathbf {f}_{1,j}) &\quad \cdots &\quad H_{d_{1}}(\mathbf {f}_{d_{2}-1,j}) \end{bmatrix}, \!\! \\\tag{21}\end{align*}
\begin{align*} H_{d_{1}}(\mathbf {f}_{m,j}) \!=\! \begin{bmatrix} \mathbf {f}_{m,j}[{1}] &\quad \mathbf {f} _{m,j}[{2}] &\quad \cdots &\quad \mathbf {f} _{m,j}[d_{1}] \\ \mathbf {f}_{m,j}[{2}] &\quad \mathbf {f} _{m,j}[{3}] &\quad \cdots &\quad \mathbf {f} _{m,j}[d_{1}+1] \\ \vdots &\quad \vdots &\quad \ddots &\quad \vdots \\ \mathbf {f}_{m,j}[N] &\quad \mathbf {f} _{m,j}[{1}] &\quad \cdots &\quad \mathbf {f} _{m,j}[d_{1}-1] \end{bmatrix},\!\! \\\tag{22}\end{align*}
\begin{equation*} H_{d_{1} d_{2}} (\mathbf {F}) = \frac {1}{\alpha } \tilde {\mathbf {\Phi }}\mathbf {\Phi }^{T} H_{d_{1} d_{2}} (\mathbf {F}) {\mathbf {\Psi }}\tilde {\mathbf {\Psi }}^{T} = \frac {1}{\alpha } \tilde {\mathbf {\Phi }}\mathbf {C}\tilde {\mathbf {\Psi }}^{T}\tag{23}\end{equation*}
