西南石油大學(xué)計算機科學(xué)學(xué)院，四川 成都 610500

醫(yī)學(xué)超聲成像是利用超聲在人體組織內(nèi)的傳播特性而成像的方式。由于醫(yī)學(xué)超聲成像的便捷、廉價、安全、成像速度快等優(yōu)點[1]，使得該技術(shù)逐漸成為醫(yī)學(xué)上廣泛應(yīng)用的成像技術(shù)之一。超聲成像技術(shù)的發(fā)展離不開超聲成像模擬技術(shù)，超聲成像模擬技術(shù)在各種超聲圖像處理算法、超聲彈性成像算法、剪切波成像算法、超聲醫(yī)師的崗前培訓(xùn)等方面具有廣泛的運用，也是研究新成像算法和估計算法的重要第一步。模擬作為測試方法的一種基本方式，主要為了在可控制的環(huán)境下證實或者否定提出的假設(shè)。由于可以在模擬中控制所有的參數(shù)，可以首先構(gòu)建一個簡單的數(shù)字體模模型，然后逐步添加其它復(fù)雜的結(jié)構(gòu)，進而模擬出逼近真實人體組織的數(shù)字體模。構(gòu)建好模型后，由于模型里的幾何結(jié)構(gòu)或者內(nèi)部不同區(qū)域的物理屬性是已知的，通過超聲模擬器掃描數(shù)字體模得到的數(shù)據(jù)可以作為檢驗和評估各種算法的手段。

超聲成像模擬整個過程可以分為兩部分，第一部分是建立人體組織的數(shù)字體模，這個過程通常采用蒙特卡洛模擬法隨機產(chǎn)生大量的散射子。人體組織是通過這些隨機分布的散射子來表示，同時根據(jù)散射子的位置設(shè)置不同的反射系數(shù)；第二部分是模擬超聲成像的過程，許多研究小組設(shè)計了不同的模擬解決方案。如基于空間沖激響應(yīng)的模擬方案，代表性的研究方案為 Jorgen Arendt Jensen 所開發(fā)的 Field II[2]，也包括Sverre Holm 小組的 Ultrasim[3]和 D’hooge 小組的方案[4]；除此之外，還有基于有限元的模擬[5]；基于 K 空間的模擬[6]；以及基于點擴散函數(shù)與散射子的卷積方法[7]等。前面提到的方案的優(yōu)點是模擬成像的效果非常逼近于真實的超聲圖像，但缺點是計算時間太長。卷積方案的計算速度快，但是缺點是并不是完全符合真實超聲成像時空間變化和聲場的復(fù)雜性。

Field II[2]是由丹麥理工大學(xué)的 Jorgen Arendt Jensen 于 1999 年開發(fā)的程序，主要用于超聲領(lǐng)域的研究，F(xiàn)ield II[2]也是目前最廣泛使用的超聲模擬器，通常為了測試各種新的超聲成像技術(shù)和方法，第一步的工作就是通過 Field II 模擬超聲數(shù)據(jù)進行測試。它在無數(shù)的研究中顯示其是一個非常準確的超聲模擬解決方案。該程序可以仿真超聲掃描儀圖像的生成，也可以模擬不同類型的超聲換能器 (比如線陣、相控陣、2D 陣等陣列) 和相關(guān)圖像，并且已經(jīng)成功在很多研究中使用并且生成了精確的結(jié)果。Field II 作為公認的權(quán)威的超聲模擬器，相對于卷積方法來模擬超聲影像有其固有的優(yōu)勢，但是過長的計算時間又使它的應(yīng)用受到了一些限制。

由于受到 CPU 主頻提高的限制，當前 CPU 的技術(shù)朝著更寬的方向發(fā)展 (即一個物理 CPU 具有更多的CPU 核)。另外，隨著 Intel 公司提出的超線程技術(shù)的發(fā)展，使得單個 CPU 擁有處理多線程的能力，能讓處理器真正實現(xiàn)多線程工作，進而充分利用處理器的性能，例如雙核四線程或者四核八線程的 CPU。由于Field II 模擬超聲成像時天然的可并行性，在多 CPU核計算機或集群上，通過調(diào)用多個 Matlab 程序獨立計算一部分超聲信號掃描線是一種通常的做法。目前Field II pro[8]是 Field II的多線程版本，但是 Field II pro 需要付費且無分布式計算功能。

本文的工作是利用 Matlab 提供的并行計算工具箱 PCT (Parallel Computing Toolbox) 和分布式計算服務(wù)器 MDCE (Matlab Distributed Computing Server) 開發(fā)一個可以利用具有多個 CPU 核的節(jié)點或集群加速Field II 模擬超聲成像過程的開源架構(gòu)，這樣的架構(gòu)可以實現(xiàn)多處理器、多核甚至集群化的高效并行計算。達到充分利用多 CPU 核計算機或集群的計算能力來減少目前最常用的超聲模擬器 Field II 的計算時間。

1 方法

1.1 Matlab 并行計算工具箱和分布式計算服務(wù)器

Matlab 分布式計算主要分為兩類：一類是由多核或者多處理器構(gòu)成的計算機平臺，一類是由多個計算結(jié)點構(gòu)成的集群系統(tǒng)。并行計算工具箱 PCT (Parallel Computing Toolbox) 提供了多種高級編程結(jié)構(gòu)，只需要簡單修改原來的代碼就可以使原來串行執(zhí)行的Matlab 代碼在多個 worker 上并行執(zhí)行，達到充分利用計算機資源和提高計算效率的目的[9]。配合 Matlab分布式計算服務(wù)器，可以在計算機集群中進行分布式并行計算。簡單來說，用戶只需通過客戶端用 PCT建立好作業(yè)任務(wù)，然后把任務(wù)提交給作業(yè)管理器 JM(Job Manager)，JM 自動劃分任務(wù)給 worker 進行計算，worker 計算完成后再將結(jié)果返回到 JM。

Matlab 的分布式計算服務(wù)器 MDCE (Matlab Distributed Computing Server) 用于支持基于集群的分布式計算[10]。該服務(wù)能夠啟動作業(yè)管理器來計算各個任務(wù)并負責把結(jié)果返回給客戶端。用戶幾乎不需要做任何改動，便可以讓并行計算工具箱 PCT 的程序得到擴展，使之運行在集群內(nèi)任意計算機的任意數(shù)量 worker 上。一個簡單的 Matlab 集群示意圖如圖1 所示。

圖1 Matlab 的集群結(jié)構(gòu)Fig. 1 Matlab cluster structure

該架構(gòu)中一共有三種角色，分別為提交 Job 并將Job 劃分為多個 task 的客戶端 (即桌面系統(tǒng))，管理作業(yè)運行的 Job Manager (JM) 或者 Scheduler，以及執(zhí)行task 具體運算的 worker。在一個分布式系統(tǒng)中通常包括多個 worker，具體數(shù)量取決于節(jié)點機核心數(shù)量，以便同時執(zhí)行多個 task。

1.2 循環(huán)并行 parfor

PCT 工具箱提供的 parfor 循環(huán)能夠在多核計算機或者集群上調(diào)用多個 worker 來處理計算任務(wù) 。使用 parfor 循環(huán)要求循環(huán)體的每次迭代都是獨立的，迭代之間沒有依賴性，執(zhí)行順序的變化對結(jié)果沒有影響。由于 parfor 使用了 MPI 作為底層通訊協(xié)議，parfor 可以跨機運行進而擴展到集群上運行，而用戶只需簡單修改原來的代碼。利用 parfor 并行 for 循環(huán)的步驟：

步驟一：使用 Matlab pool 函數(shù)配置和開啟Matalb 并行池；

步驟二：將無數(shù)據(jù)依賴的串行 for 循環(huán)中改為parfor 循環(huán)，并注意是否需要修改循環(huán)體；

步驟三：執(zhí)行完畢后使用 delete 關(guān)閉并行池。

Matlab 對于 parfor 關(guān)鍵字的處理是將任務(wù)動態(tài)分配到多個 worker 上并行執(zhí)行。假設(shè) n 為集群中 worker 的數(shù)量，m 為循環(huán)次數(shù)，如果 m/n 為整數(shù)，則循環(huán)被均勻分配到每個 worker；否則循環(huán)分配不均，某些 worker 則負載過大。所以，為了提高并發(fā)執(zhí)行的效率，應(yīng)該保證設(shè)置的 m 和 n 是整除關(guān)系。

1.3 超聲模擬成像過程

利用 Field II 進行人體組織模擬成像的流程主要包括系統(tǒng)的初始化、探頭的設(shè)置、陣元的設(shè)定、散射體模型的建立、聲場的計算五個步驟。超聲成像模擬的第一步在于計算機數(shù)字體模 (即散射體模型) 的建立，通過在文件中定義一系列散射子的位置和相應(yīng)的反射系數(shù)，進而計算掃描線來模擬生成B-mode 圖像，調(diào)整相關(guān)散射子的反射系數(shù)可以模擬人體中不同位置的組織。如果組織具有病變，則該組織的反射系數(shù)將大于周圍的其它組織；如果某一部分主要由液體構(gòu)成，那么這部分的反射系數(shù)將較小。由于這種基于蒙特卡洛方法會產(chǎn)生大量的隨機分布的散射子，導(dǎo)致整個模擬 RF 掃描線的過程會花費非常長的時間，即使在當前最快的 CPU，通常也需要數(shù)小時完成一幅 B-Mode 圖像所需的所有掃描線的模擬工作。如果模擬三維的超聲數(shù)據(jù)，其計算時間將會更長。

1.4 算法并行性分析與實現(xiàn)

Field II 超聲模擬器模擬超聲 RF 掃描線的過程根據(jù)探頭類型、散射子數(shù)量、掃描線數(shù)量等參數(shù)的變化，計算時間從十幾分鐘到數(shù)小時不等。但是不同掃描線之間的計算是獨立的、沒有數(shù)據(jù)依賴。為了加速這個過程并充分利用現(xiàn)有計算機的計算資源，通過PCT 工具箱和分布式計算服務(wù)器 MDCE 是一種簡單快捷的方法。主要思路是由 PCT 工具箱根據(jù)計算機計算資源的數(shù)量分配相應(yīng)的 worker，每一個 worker在某一時刻執(zhí)行一條 RF 掃描線的模擬任務(wù)，每一個worker 內(nèi)部需要初始化一個 Field II 程序的實例并執(zhí)行相應(yīng)的計算。

圖2 所示的流程圖演示了整個計算過程。首先讀入計算機數(shù)字體模數(shù)據(jù)，初始化探頭基本參數(shù)，然后打開 Matlab 并行池，并設(shè)置啟動的 worker 數(shù)目，parfor 循環(huán)將計算 RF 線的任務(wù)分配到集群上不同worker 上執(zhí)行。待 RF 線全部計算完成后將并行池關(guān)閉，并取出計算結(jié)果生成 B 模式圖像。

圖2 超聲模擬的并行計算架構(gòu)Fig. 2 Parallel computing architecture of ultrasonic simulating

2 實驗與討論

實驗中操作系統(tǒng)為 Windows Server 2008 R2，集群由三臺計算節(jié)點組成，節(jié)點的配置為：Intel Core Xeon CPU E5-4607 2.2GHz (2 處理器 12 核心 24 線程，支持超線程) 內(nèi)存 32G。為保證程序的兼容性以及節(jié)點機之間正常通信，均安裝 Matlab R2012b 版本，且安裝在同一路徑下。安裝完成后，首先在各個節(jié)點機上安裝 MDCE 服務(wù)，通過命令中心加入各個節(jié)點機，然后啟動 Job Manager 并完成集群計算節(jié)點的配置。該集群每臺節(jié)點機配置 24 個 worker，集群總共擁有 72 個 worker。實驗中為了各個 worker 能調(diào)用 Field II 程序，需要將 Field II 程序安裝到同一路徑下。

2.1 實驗結(jié)果

本文中所實現(xiàn)的基于 Matlab 的并行與分布式計算的 Field II 超聲模擬器的性能通過兩組實驗測試完成。兩組實驗均通過對 Field II 程序開發(fā)者提供的模擬示例改寫完成，采用兩個節(jié)點機進行計算。第一組實驗的仿真體模為一個 50mm × 60mm × 10 mm (寬度× 深度 × 厚度) 的三維結(jié)構(gòu)。由 10 萬個散射子構(gòu)成。實驗?zāi)M采用線陣探頭進行掃描，總共生成 64 條 RF線。第二組實驗的仿真體模為一個 100 mm × 100 mm ×15 mm (寬度 × 深度 × 厚度) 的三維結(jié)構(gòu)。由 100 萬個散射子構(gòu)成，實驗?zāi)M采用相控陣探頭進行掃描，總共生成 128 條 RF 線。兩次實驗均通過蒙特卡洛法創(chuàng)建隨機數(shù)用于模擬人體組織細小結(jié)構(gòu)的散射子，其強度均符合高斯分布。模擬所需時間的測定通過 Matlab自帶的 tic/toc 完成，采用 10 次計算平均值的方式獲得每一次不同線程模擬計算的耗時結(jié)果。表 1 中列出了兩組實驗時探頭的主要參數(shù)。

表1 Field II 模擬的探頭參數(shù)Table 1 Field II simulating parameter of transducer

圖3 所示第一組實驗的模擬結(jié)果。圖中，x 軸表示的是本文程序中啟動的線程，y 軸是原始的 Field II 程序相比的加速比。對比于原始的 Field II 程序，本文所開發(fā)的并行與分布式版本的加速比大致上與所使用線程成正比。實驗中啟動線程的數(shù)量依次為 2，4，8，16，32。由于生成的 RF 線為 64 條，為了使計算任務(wù)均勻分配，防止部分 worker 負載過大，線程數(shù)量設(shè)置能被 RF 線數(shù)量整除。當采用 32 個線程 (對應(yīng)為 32 個 CPU 核) 時，對于 10 萬個散射子的模擬所需要的時間為 145 秒，加速比為 22。

表2 為整個超聲成像模擬流程各階段所需時間?？梢詮谋?2 中明顯看出 RF 線的計算占據(jù)了模擬計算流程的主體，加載體模和初始化探頭參數(shù)的時間幾乎可以忽略不計。由于集群各節(jié)點間通信延遲以及網(wǎng)絡(luò)帶寬的影響，啟動 32 個 worker 需要消耗 10 秒左右的時間。所以，對于 RF 線的并行加速大幅度的提高了超聲成像模擬的整體速度。

圖3 本文實現(xiàn)的并行與分布式計算 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 所示為第二組實驗的結(jié)果。與圖 3 相同，x 軸和 y 軸分別表示啟動的線程和加速比。模擬采用的探頭為具有 128 陣列的相控陣探頭，中心頻率為 7 MHZ。成像的范圍為 -45 至 45 度的范圍。實驗中啟動線程的數(shù)量與第一組實驗相同。當采用 32 個線程(對應(yīng)為 32 個 CPU 核) 時，對于 100 萬個散射子的模擬所需要的時間為 4910 秒。加速比為 25.5。

圖3 和圖 4 中紅色的線顯示的為理想狀態(tài)下的加速比，即完全線性相關(guān)模擬計算所用 CPU 核的數(shù)量。本文所開發(fā)的并行與分布式版本的加速比低于理想狀態(tài)的加速比的主要原因是當采用的線程數(shù)量大于單個節(jié)點所能提供的 CPU 核數(shù)量時，程序會自動分配剩余的計算任務(wù)到其它計算節(jié)點上進行計算，這個過程會產(chǎn)生數(shù)據(jù)傳輸及任務(wù)管理上的開銷。由于本文的系統(tǒng)節(jié)點之間的數(shù)據(jù)傳輸采用的還是 100M 以太網(wǎng)結(jié)構(gòu)，這也是導(dǎo)致任務(wù)分發(fā)和結(jié)果回收時影響。當分配的線程數(shù)量如果大于整個集群所能提供的物理 CPU核數(shù)量時，這時性能還將會進一步的下降?？傮w上來說，當分配的線程數(shù)量不超過單個節(jié)點的物理 CPU核時，隨著線程數(shù)量的增加，加速比會逐漸增大并接近于理想的加速比，但趨于逐漸下降趨勢。

圖5 演示了不同散射子規(guī)模對模擬時間的影響。如圖所示，散射子數(shù)量對于模擬計算的時間開銷具有明顯的影響。當計算線程設(shè)置為 32 時，模擬 10 萬和 100 萬個散射子所構(gòu)成的體模所需要的時間分別為 145 秒和 4910 秒，由此可以發(fā)現(xiàn)，模擬計算并不是完全隨著散射子的規(guī)模進行倍數(shù)的增長，而是與探頭、散射子數(shù)量這個因素綜合的結(jié)果。

圖4 本文實現(xiàn)的并行與分布式計算 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 模擬不同數(shù)量散射子時本文實現(xiàn)的并行與分布式計算Field II 的執(zhí)行時間F i g. 5 Executing time of realized parallel and distributed computing of fi eld II on different amount scatters

為進一步的比較集群規(guī)模對于模擬時間的影響，本文在三節(jié)點平臺下對第一個實驗體模采用 36 核和72 核模擬生成 144 條 RF 線，所需時間分別為 341.8秒和 237.5 秒，從它們所需要的計算時間來看，加速比與所期望的線性或近似于線性的增長趨勢相差較遠。這也說明，當單條 RF 線的模擬時間較少的時候，采用更多的節(jié)點或線程進行模擬計算，計算數(shù)據(jù)的發(fā)送和結(jié)果的回收將影響到加速比的提升。

圖6 演示了一個復(fù)雜的數(shù)字乳腺體模[11]的模擬結(jié)果。這個數(shù)字乳腺體模為一個40 mm × 23 mm × 8 mm (寬度 × 深度 × 厚度) 的三維結(jié)構(gòu)，由 5 萬個散射子模擬構(gòu)成。模擬實驗采用線陣探頭進行掃描，主要參數(shù)設(shè)置可參考[11]，總共生成 256 條 RF 線。計算線程設(shè)置為 32，模擬總耗時為 145 秒。由于這個復(fù)雜的數(shù)字乳腺體模的內(nèi)部結(jié)構(gòu)是已知的，通過掃描這樣的一個復(fù)雜的數(shù)字乳腺體模生成的RF信號或者生成 B-mode 影像，可以幫助對現(xiàn)有的如醫(yī)學(xué)圖像分割算法或者其它醫(yī)學(xué)圖像處理算法進行評估與改進，同時結(jié)合有限元模擬，也可以進一步用于超聲彈性成像及剪切波彈性算法的評估。

圖6 復(fù)雜的計算機乳腺數(shù)字體膜的模擬結(jié)果F i g. 6 Simulating results of complicated computer breast digital phantom

4 結(jié)論

超聲成像的模擬通常需要較長的時間，由于超聲成像的模擬過程具有可并行計算的特性，因此結(jié)合現(xiàn)在個人計算機、工作站及計算機集群的多 CPU 架構(gòu)，有必要開發(fā)可并行和分布式計算的開放獲取的超聲模擬程序。Field II 作為當前學(xué)術(shù)界最廣泛使用的超聲模擬器，雖然其開發(fā)者也提供了一個可并行計算的新版本 Field II pro，但需要付費使用且無分布式計算功能。通過利用 Matlab 的并行計算工具箱和分布式計算服務(wù)器，只需要對原始的代碼進行簡單的修改就可以進行并行與分布式模擬計算。本文所實現(xiàn)的方法的性能通過模擬線陣和相控陣探頭的相關(guān)參數(shù)，結(jié)合用于產(chǎn)生計算機數(shù)字體模所需要的大量的隨機散射子被測試。測試過程中，主要測試了不同的線程數(shù)量下，線陣和相控陣探頭在不同的散射子規(guī)模下的計算時間。從測試的結(jié)果來看，本文所實現(xiàn)的方法的加速比大致上與所采用的 CPU 核成正比。因此，本文所實現(xiàn)的方法在性能上的提升主要是靠多線程執(zhí)行。本文所開發(fā)的方法將會是一個完全開源的框架，整個計算與原來的 Field II 代碼兼容，只需要進行簡單的修改和并行計算與分布式計算部署就可以實現(xiàn)。程序中可以根據(jù)個人或集群計算機物理CPU 配置，進行線程配置，計算任務(wù)將自動的進行動態(tài)劃分。本文所開發(fā)的方法將有助于加速超聲成像模擬的過程而不需要購買付費的 Field II pro，加速超聲成像新方法的研究。

致謝

感謝美國密歇根理工大學(xué)生物醫(yī)學(xué)工程系博士生王禹提供用于生成圖 6 的數(shù)字乳腺體模及單線程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.