Field II 超声模拟器的并行与分布式计算研究

西南石油大学计算机科学学院，四川 成都 610500

医学超声成像是利用超声在人体组织内的传播特性而成像的方式。由于医学超声成像的便捷、廉价、安全、成像速度快等优点[1]，使得该技术逐渐成为医学上广泛应用的成像技术之一。超声成像技术的发展离不开超声成像模拟技术，超声成像模拟技术在各种超声图像处理算法、超声弹性成像算法、剪切波成像算法、超声医师的岗前培训等方面具有广泛的运用，也是研究新成像算法和估计算法的重要第一步。模拟作为测试方法的一种基本方式，主要为了在可控制的环境下证实或者否定提出的假设。由于可以在模拟中控制所有的参数，可以首先构建一个简单的数字体模模型，然后逐步添加其它复杂的结构，进而模拟出逼近真实人体组织的数字体模。构建好模型后，由于模型里的几何结构或者内部不同区域的物理属性是已知的，通过超声模拟器扫描数字体模得到的数据可以作为检验和评估各种算法的手段。

超声成像模拟整个过程可以分为两部分，第一部分是建立人体组织的数字体模，这个过程通常采用蒙特卡洛模拟法随机产生大量的散射子。人体组织是通过这些随机分布的散射子来表示，同时根据散射子的位置设置不同的反射系数；第二部分是模拟超声成像的过程，许多研究小组设计了不同的模拟解决方案。如基于空间冲激响应的模拟方案，代表性的研究方案为 Jorgen Arendt Jensen 所开发的 Field II[2]，也包括Sverre Holm 小组的 Ultrasim[3]和 D’hooge 小组的方案[4]；除此之外，还有基于有限元的模拟[5]；基于 K 空间的模拟[6]；以及基于点扩散函数与散射子的卷积方法[7]等。前面提到的方案的优点是模拟成像的效果非常逼近于真实的超声图像，但缺点是计算时间太长。卷积方案的计算速度快，但是缺点是并不是完全符合真实超声成像时空间变化和声场的复杂性。

Field II[2]是由丹麦理工大学的 Jorgen Arendt Jensen 于 1999 年开发的程序，主要用于超声领域的研究，Field II[2]也是目前最广泛使用的超声模拟器，通常为了测试各种新的超声成像技术和方法，第一步的工作就是通过 Field II 模拟超声数据进行测试。它在无数的研究中显示其是一个非常准确的超声模拟解决方案。该程序可以仿真超声扫描仪图像的生成，也可以模拟不同类型的超声换能器 (比如线阵、相控阵、2D 阵等阵列) 和相关图像，并且已经成功在很多研究中使用并且生成了精确的结果。Field II 作为公认的权威的超声模拟器，相对于卷积方法来模拟超声影像有其固有的优势，但是过长的计算时间又使它的应用受到了一些限制。

由于受到 CPU 主频提高的限制，当前 CPU 的技术朝着更宽的方向发展 (即一个物理 CPU 具有更多的CPU 核)。另外，随着 Intel 公司提出的超线程技术的发展，使得单个 CPU 拥有处理多线程的能力，能让处理器真正实现多线程工作，进而充分利用处理器的性能，例如双核四线程或者四核八线程的 CPU。由于Field II 模拟超声成像时天然的可并行性，在多 CPU核计算机或集群上，通过调用多个 Matlab 程序独立计算一部分超声信号扫描线是一种通常的做法。目前Field II pro[8]是 Field II的多线程版本，但是 Field II pro 需要付费且无分布式计算功能。

本文的工作是利用 Matlab 提供的并行计算工具箱 PCT (Parallel Computing Toolbox) 和分布式计算服务器 MDCE (Matlab Distributed Computing Server) 开发一个可以利用具有多个 CPU 核的节点或集群加速Field II 模拟超声成像过程的开源架构，这样的架构可以实现多处理器、多核甚至集群化的高效并行计算。达到充分利用多 CPU 核计算机或集群的计算能力来减少目前最常用的超声模拟器 Field II 的计算时间。

1 方法

1.1 Matlab 并行计算工具箱和分布式计算服务器

Matlab 分布式计算主要分为两类：一类是由多核或者多处理器构成的计算机平台，一类是由多个计算结点构成的集群系统。并行计算工具箱 PCT (Parallel Computing Toolbox) 提供了多种高级编程结构，只需要简单修改原来的代码就可以使原来串行执行的Matlab 代码在多个 worker 上并行执行，达到充分利用计算机资源和提高计算效率的目的[9]。配合 Matlab分布式计算服务器，可以在计算机集群中进行分布式并行计算。简单来说，用户只需通过客户端用 PCT建立好作业任务，然后把任务提交给作业管理器 JM(Job Manager)，JM 自动划分任务给 worker 进行计算，worker 计算完成后再将结果返回到 JM。

Matlab 的分布式计算服务器 MDCE (Matlab Distributed Computing Server) 用于支持基于集群的分布式计算[10]。该服务能够启动作业管理器来计算各个任务并负责把结果返回给客户端。用户几乎不需要做任何改动，便可以让并行计算工具箱 PCT 的程序得到扩展，使之运行在集群内任意计算机的任意数量 worker 上。一个简单的 Matlab 集群示意图如图1 所示。

图1 Matlab 的集群结构Fig. 1 Matlab cluster structure

该架构中一共有三种角色，分别为提交 Job 并将Job 划分为多个 task 的客户端 (即桌面系统)，管理作业运行的 Job Manager (JM) 或者 Scheduler，以及执行task 具体运算的 worker。在一个分布式系统中通常包括多个 worker，具体数量取决于节点机核心数量，以便同时执行多个 task。

1.2 循环并行 parfor

PCT 工具箱提供的 parfor 循环能够在多核计算机或者集群上调用多个 worker 来处理计算任务 。使用 parfor 循环要求循环体的每次迭代都是独立的，迭代之间没有依赖性，执行顺序的变化对结果没有影响。由于 parfor 使用了 MPI 作为底层通讯协议，parfor 可以跨机运行进而扩展到集群上运行，而用户只需简单修改原来的代码。利用 parfor 并行 for 循环的步骤：

步骤一：使用 Matlab pool 函数配置和开启Matalb 并行池；

步骤二：将无数据依赖的串行 for 循环中改为parfor 循环，并注意是否需要修改循环体；

步骤三：执行完毕后使用 delete 关闭并行池。

Matlab 对于 parfor 关键字的处理是将任务动态分配到多个 worker 上并行执行。假设 n 为集群中 worker 的数量，m 为循环次数，如果 m/n 为整数，则循环被均匀分配到每个 worker；否则循环分配不均，某些 worker 则负载过大。所以，为了提高并发执行的效率，应该保证设置的 m 和 n 是整除关系。

1.3 超声模拟成像过程

利用 Field II 进行人体组织模拟成像的流程主要包括系统的初始化、探头的设置、阵元的设定、散射体模型的建立、声场的计算五个步骤。超声成像模拟的第一步在于计算机数字体模 (即散射体模型) 的建立，通过在文件中定义一系列散射子的位置和相应的反射系数，进而计算扫描线来模拟生成B-mode 图像，调整相关散射子的反射系数可以模拟人体中不同位置的组织。如果组织具有病变，则该组织的反射系数将大于周围的其它组织；如果某一部分主要由液体构成，那么这部分的反射系数将较小。由于这种基于蒙特卡洛方法会产生大量的随机分布的散射子，导致整个模拟 RF 扫描线的过程会花费非常长的时间，即使在当前最快的 CPU，通常也需要数小时完成一幅 B-Mode 图像所需的所有扫描线的模拟工作。如果模拟三维的超声数据，其计算时间将会更长。

1.4 算法并行性分析与实现

Field II 超声模拟器模拟超声 RF 扫描线的过程根据探头类型、散射子数量、扫描线数量等参数的变化，计算时间从十几分钟到数小时不等。但是不同扫描线之间的计算是独立的、没有数据依赖。为了加速这个过程并充分利用现有计算机的计算资源，通过PCT 工具箱和分布式计算服务器 MDCE 是一种简单快捷的方法。主要思路是由 PCT 工具箱根据计算机计算资源的数量分配相应的 worker，每一个 worker在某一时刻执行一条 RF 扫描线的模拟任务，每一个worker 内部需要初始化一个 Field II 程序的实例并执行相应的计算。

图2 所示的流程图演示了整个计算过程。首先读入计算机数字体模数据，初始化探头基本参数，然后打开 Matlab 并行池，并设置启动的 worker 数目，parfor 循环将计算 RF 线的任务分配到集群上不同worker 上执行。待 RF 线全部计算完成后将并行池关闭，并取出计算结果生成 B 模式图像。

图2 超声模拟的并行计算架构Fig. 2 Parallel computing architecture of ultrasonic simulating

2 实验与讨论

实验中操作系统为 Windows Server 2008 R2，集群由三台计算节点组成，节点的配置为：Intel Core Xeon CPU E5-4607 2.2GHz (2 处理器 12 核心 24 线程，支持超线程) 内存 32G。为保证程序的兼容性以及节点机之间正常通信，均安装 Matlab R2012b 版本，且安装在同一路径下。安装完成后，首先在各个节点机上安装 MDCE 服务，通过命令中心加入各个节点机，然后启动 Job Manager 并完成集群计算节点的配置。该集群每台节点机配置 24 个 worker，集群总共拥有 72 个 worker。实验中为了各个 worker 能调用 Field II 程序，需要将 Field II 程序安装到同一路径下。

2.1 实验结果

本文中所实现的基于 Matlab 的并行与分布式计算的 Field II 超声模拟器的性能通过两组实验测试完成。两组实验均通过对 Field II 程序开发者提供的模拟示例改写完成，采用两个节点机进行计算。第一组实验的仿真体模为一个 50mm × 60mm × 10 mm (宽度× 深度 × 厚度) 的三维结构。由 10 万个散射子构成。实验模拟采用线阵探头进行扫描，总共生成 64 条 RF线。第二组实验的仿真体模为一个 100 mm × 100 mm ×15 mm (宽度 × 深度 × 厚度) 的三维结构。由 100 万个散射子构成，实验模拟采用相控阵探头进行扫描，总共生成 128 条 RF 线。两次实验均通过蒙特卡洛法创建随机数用于模拟人体组织细小结构的散射子，其强度均符合高斯分布。模拟所需时间的测定通过 Matlab自带的 tic/toc 完成，采用 10 次计算平均值的方式获得每一次不同线程模拟计算的耗时结果。表 1 中列出了两组实验时探头的主要参数。

表1 Field II 模拟的探头参数Table 1 Field II simulating parameter of transducer

图3 所示第一组实验的模拟结果。图中，x 轴表示的是本文程序中启动的线程，y 轴是原始的 Field II 程序相比的加速比。对比于原始的 Field II 程序，本文所开发的并行与分布式版本的加速比大致上与所使用线程成正比。实验中启动线程的数量依次为 2，4，8，16，32。由于生成的 RF 线为 64 条，为了使计算任务均匀分配，防止部分 worker 负载过大，线程数量设置能被 RF 线数量整除。当采用 32 个线程 (对应为 32 个 CPU 核) 时，对于 10 万个散射子的模拟所需要的时间为 145 秒，加速比为 22。

表2 为整个超声成像模拟流程各阶段所需时间。可以从表 2 中明显看出 RF 线的计算占据了模拟计算流程的主体，加载体模和初始化探头参数的时间几乎可以忽略不计。由于集群各节点间通信延迟以及网络带宽的影响，启动 32 个 worker 需要消耗 10 秒左右的时间。所以，对于 RF 线的并行加速大幅度的提高了超声成像模拟的整体速度。

图3 本文实现的并行与分布式计算 Field II 与原始的 Field II 的加速比 (模拟线阵探头)F i g. 3 Realized parallel and distributed computing of fi eld II compared with original fi eld II on ratio (simulating linear array transducer)

表2 超声成像模拟时间开销Table 2 Ultrasound imaging simulating time

图4 所示为第二组实验的结果。与图 3 相同，x 轴和 y 轴分别表示启动的线程和加速比。模拟采用的探头为具有 128 阵列的相控阵探头，中心频率为 7 MHZ。成像的范围为 -45 至 45 度的范围。实验中启动线程的数量与第一组实验相同。当采用 32 个线程(对应为 32 个 CPU 核) 时，对于 100 万个散射子的模拟所需要的时间为 4910 秒。加速比为 25.5。

图3 和图 4 中红色的线显示的为理想状态下的加速比，即完全线性相关模拟计算所用 CPU 核的数量。本文所开发的并行与分布式版本的加速比低于理想状态的加速比的主要原因是当采用的线程数量大于单个节点所能提供的 CPU 核数量时，程序会自动分配剩余的计算任务到其它计算节点上进行计算，这个过程会产生数据传输及任务管理上的开销。由于本文的系统节点之间的数据传输采用的还是 100M 以太网结构，这也是导致任务分发和结果回收时影响。当分配的线程数量如果大于整个集群所能提供的物理 CPU核数量时，这时性能还将会进一步的下降。总体上来说，当分配的线程数量不超过单个节点的物理 CPU核时，随着线程数量的增加，加速比会逐渐增大并接近于理想的加速比，但趋于逐渐下降趋势。

图5 演示了不同散射子规模对模拟时间的影响。如图所示，散射子数量对于模拟计算的时间开销具有明显的影响。当计算线程设置为 32 时，模拟 10 万和 100 万个散射子所构成的体模所需要的时间分别为 145 秒和 4910 秒，由此可以发现，模拟计算并不是完全随着散射子的规模进行倍数的增长，而是与探头、散射子数量这个因素综合的结果。

图4 本文实现的并行与分布式计算 Field II 与原始的 Field II 的加速比 (模拟相控阵探头)F i g. 4 Realized parallel and distributed computing of fi eld II compared with original fi eld II on ratio (simulating linear array transduce)

图5 模拟不同数量散射子时本文实现的并行与分布式计算Field II 的执行时间F i g. 5 Executing time of realized parallel and distributed computing of fi eld II on different amount scatters

为进一步的比较集群规模对于模拟时间的影响，本文在三节点平台下对第一个实验体模采用 36 核和72 核模拟生成 144 条 RF 线，所需时间分别为 341.8秒和 237.5 秒，从它们所需要的计算时间来看，加速比与所期望的线性或近似于线性的增长趋势相差较远。这也说明，当单条 RF 线的模拟时间较少的时候，采用更多的节点或线程进行模拟计算，计算数据的发送和结果的回收将影响到加速比的提升。

图6 演示了一个复杂的数字乳腺体模[11]的模拟结果。这个数字乳腺体模为一个40 mm × 23 mm × 8 mm (宽度 × 深度 × 厚度) 的三维结构，由 5 万个散射子模拟构成。模拟实验采用线阵探头进行扫描，主要参数设置可参考[11]，总共生成 256 条 RF 线。计算线程设置为 32，模拟总耗时为 145 秒。由于这个复杂的数字乳腺体模的内部结构是已知的，通过扫描这样的一个复杂的数字乳腺体模生成的RF信号或者生成 B-mode 影像，可以帮助对现有的如医学图像分割算法或者其它医学图像处理算法进行评估与改进，同时结合有限元模拟，也可以进一步用于超声弹性成像及剪切波弹性算法的评估。

图6 复杂的计算机乳腺数字体膜的模拟结果F i g. 6 Simulating results of complicated computer breast digital phantom

4 结论

超声成像的模拟通常需要较长的时间，由于超声成像的模拟过程具有可并行计算的特性，因此结合现在个人计算机、工作站及计算机集群的多 CPU 架构，有必要开发可并行和分布式计算的开放获取的超声模拟程序。Field II 作为当前学术界最广泛使用的超声模拟器，虽然其开发者也提供了一个可并行计算的新版本 Field II pro，但需要付费使用且无分布式计算功能。通过利用 Matlab 的并行计算工具箱和分布式计算服务器，只需要对原始的代码进行简单的修改就可以进行并行与分布式模拟计算。本文所实现的方法的性能通过模拟线阵和相控阵探头的相关参数，结合用于产生计算机数字体模所需要的大量的随机散射子被测试。测试过程中，主要测试了不同的线程数量下，线阵和相控阵探头在不同的散射子规模下的计算时间。从测试的结果来看，本文所实现的方法的加速比大致上与所采用的 CPU 核成正比。因此，本文所实现的方法在性能上的提升主要是靠多线程执行。本文所开发的方法将会是一个完全开源的框架，整个计算与原来的 Field II 代码兼容，只需要进行简单的修改和并行计算与分布式计算部署就可以实现。程序中可以根据个人或集群计算机物理CPU 配置，进行线程配置，计算任务将自动的进行动态划分。本文所开发的方法将有助于加速超声成像模拟的过程而不需要购买付费的 Field II pro，加速超声成像新方法的研究。

致谢

感谢美国密歇根理工大学生物医学工程系博士生王禹提供用于生成图 6 的数字乳腺体模及单线程FIELD 模拟代码。

[1]Jerrold T. Bushberg, John M. Boone. The essential physics of medical imaging: Lippincott Williams & Wilkins, 2011.

[2]Jørgen Arendt Jensen. Field: A program for simulating ultrasound systems: Citeseer, 1996.

[3]Jon Petter Åsen, Sverre Holm. Huygens on speed: Interactive simulation of ultrasound pressure fi elds: IEEE, 2012.

[4]Jan D Hooge, Johan Nuyts, Bart Bijnenset al. The calculation of the transient near and far fi eld of a baf fl ed piston using low sampling frequencies. The Journal of the Acoustical Society of America, 1997, 102 (1): 78-86.

[5]Gianmarco Pinton, Gregg Trahey. 3C-2 Full-Wave Simulation of Finite-Amplitude Ultrasound in Heterogeneous Media: IEEE, 2007.

[6]Bradley E. Treeby, Jiri Jaros, Alistair P. Rendellet al. Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using ak-space pseudospectral method. The Journal of the Acoustical Society of America, 2012, 131 (6): 4324-4336.

[7]Sjur Urdson Gjerald, Reidar Brekken, Torbjorn Hergumet al. Real-time ultrasound simulation using the GPU. IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 2012, 59 (5): 885-892.

[8]Jørgen Arendt Jensen. A multi-threaded version of Field II: IEEE, 2014.

[9]Vinay K. Ingle, John G. Proakis. Digital Signal Processing Using MATLAB: A Problem Solving Companion:Cengage Learning, 2016.

[10]Xianhua Ding, Xiaoming Wang, Yu Duet al. A reliable platform using matlab distributed computing server integrated with AIM: IEEE, 2014.

[11]Yu Wang, Emily Helminen, Jingfeng Jiang. Building a virtual simulation platform for quasistatic breast ultrasound elastography using open source software: A preliminary investigation. Medical physics, 2015, 42 (9): 5453-5466.