基于網絡計算的心電仿真計算研究

李 波，沈文楓，鄭衍衡，徐煒民

（上海大學 計算機工程和科學學院，上海200072）

0 引 言

心電仿真計算需要高性能計算環境，心電計算機仿真計算是心臟生物電現象研究中的一個重要課題［1］，隨著計算機計算能力的提高，心電活動的計算機仿真也隨著迅速發展。但是對于心臟模型的仿真計算，要求在空間上達到高分辨率、在計算上達到高精度，需要巨大的計算量，必須采用并行計算［2］。到目前為止，絕大部分的這類計算是在高性能計算機系統中采用（message passing interface，MPI）。GPU用于高性能計算，使得在高校和研究機構內部大量的PC機在計算性能上有很大的提高。利用這些PC機的空閑時間進行計算，可以代替昂貴的高性能計算機，廉價的、快速的實現心電仿真計算。

本文提出采用由PC機中的CPU和GPU組成異構網絡計算的方案，利用空閑的PC資源，實現心電仿真計算，并且提出了面向短的響應時間要求的調度算法，取得了很好的效果。

1 相關工作

Tusscher等進行心室傳導系統的模擬模型，采用C＋＋和MPI，是在由10臺Dell650Precision服務器，每臺服務器用雙核Intel Xeon 2.66GHz CPU組成的Beowolf集群式系統上運行的［3］。Reumann et al。在配置有512、1024、2048、4096個節點的IBM Blue Gene／L上，進行了大規模的心臟模擬［4］，這兩篇文章表明心電仿真主要在昂貴的高性能計算機系統上進行，在器件上都沒有采用GPU。本文研究在由CPU和GPU組成的網絡計算系統上運行心電仿真。Wenfeng Shen et al.研究了在多核和GPU上進行心電仿真的并行計算［1，5］，是在多核CPU和GPU組成的單機系統中的并行處理。本文研究的是CPU和GPU在多機系統中的并行處理。

在本文所采用的網絡計算系統的體系結構上，和著名的Bonic系統是相近的，但是有兩個不同點：①本文所研究的系統對于響應時間方面的要求比較高，這是和Bonic系統的主要差別。Bonic系統的調度程序以高吞吐率為出發點，追求公平性和高效率，在Bonic系統中，如果引入短的響應時間的要求，會影響公平性和高效率［6］。我們研究的心電仿真計算的要求是較短的響應時間，因此對于調度算法提出了不同的要求。②計算規模比較小，主要應用在高校和研究單位內部，對于提交的作業采取先提交先服務的策略，在運行時，獨占全部可用資源，要求每一個作業的運行時間盡量短，可以減少作業的平均等待時間，相當于提高了系統的計算速度［7－8］。同時，計算資源來自部門內部，可以獲得各個處理器的計算性能，為設計面向較短的響應時間的調度程序提供了條件。

2 心電仿真計算原理及其算法

心臟的計算機模型一般是通過對心臟斷面圖像離散化并進行空間插補，本論文所采用實驗模型為日本魏氏全心臟模型［9－10］，該模型是由一個56＊56＊90的斜交坐標系構成。模型上部為心房，下部為心室，模型細胞在空間上是等距的，直徑為1.5ms。模型的實際細胞數據為5萬個，通過使模型細胞具有方向性和具有和實際細胞相近的電生理特性，在解剖學和電生理學兩個方面對模型作進一步的修飾使得模型更加接近實際的心臟。

2.1 心電仿真算法

心電計算機仿真的核心便是心臟興奮序列的計算機仿真，所用算法是其關鍵所在，直接決定了計算的速度和仿真的效果，并影響到模型的實用性和仿真結果。在心電興奮傳播的過程中，算法的目標是在每一個計算過程中確定未興奮的心肌細胞時候接受臨近的已興奮心肌細胞的興奮傳播而開始其自身的興奮過程。

心電計算的仿真算法主要分為三步：第一步是以惠更斯原理模型為基礎，仿真心電興奮傳播的過程。第二步是根據傳播過程中記錄的偶極子數目，計算心電圖的容積導體和體表電勢。第三步是保存計算結果。其中第二步計算容積導體和體表電勢占用了整個計算機仿真過程的絕大部分時間。在體表軀干表面，心內膜和心外膜等四邊邊界大約3800三角形元素的電勢計算，產生大概50000個電偶極子源。容積導體中每個三角元素是異勢的。根據特定方程，從當前偶極子源采用邊界元方法計算的體表電位，心外膜和心內膜。如圖1所示。

電勢計算方程如下

圖1 人體軀干

式中：（rk，j）——在第k個表面上的第j個三角的體表電勢值，rk，j——中央位置向量。－號和＋號分別代表內部和外部信號。M則代表細胞模型的當前偶極子數目，J代表偶極子源。在整個容積導體中體表、體外膜表面、右心內膜表面和左心內膜表面分別被劃分為344，1002，307，278個節點和684，2002，610，552個三角形。心電仿真計算主要是對方程（1）進行計算機分解計算，求得計算結果。

2.2 心電仿真計算流程圖

心電仿真計算流程如圖2所示。

圖2 心電仿真計算流程

通過實驗驗證，對算法中每一步計算作時間分析，得到如下結果：在第一步的傳播過程中，所占時間為整個計算機仿真過程的0.1%，第二步的電勢計算占總仿真時間的99.7%，保存結果大概占總時間的0.2%。因此主要從第二步計算電勢來考慮，尋求最優并行化。第二步是根據傳播過程中記錄的偶極子數目，計算心電圖的容積導體和體表電位。由于心電計算的特殊性質，在每3ms的時間周期內所計算的偶極子數目不盡相同。每個周期的電勢計算只和這個周期產生的偶極子數目相關，一般情況下產生的偶極子數目越多，電勢計算所花費時間就越多。而偶極子數目在第一階段的興奮性傳播過程中便可以統計出來，通過合理的分析可以預測計算負載，為調度算法創造了條件。

3 異構網絡計算的結構

3.1 心電仿真計算在計算結構上的特點分析

心電仿真計算在計算結構上有以下的特點：高度并行，占總仿真時間99.7%的計算是高度并行的，計算模塊之間的相互通信量非常少，各個模塊完全可以獨立計算。非常適合利用網絡計算。各個模塊的計算量是可以預測的。為在網絡計算中實現定量調度提供了條件。

3.2 異構網絡結構介紹

異構網絡結構如圖3所示。

圖3 異構網絡結構

整個網絡由多臺的四核處理器、八核處理器、GPU處理器組成。這個網絡環境在高校的實驗室和研究所的研究室還是有比較大的代表性。

4 調度算法

調度算法的目標是達到按照各個服務器處理速度的比例分配任務，達到這個目標的方案很多，最簡單的方案是依照處理器的計算能力按比例進行分配，但是GPU用作高性能計算時，它的性能是隨著被處理的程序的并行程度變化的，因此在調度程序中要盡量使得計算量大（并行度高）的模塊分配給GPU計算。圖4是Geforce440GPU的處理速度隨著處理模塊的偶極子數量變化曲線。

調度算法的執行步驟如下：

圖4 Geforce440GPU的處理速度與處理模塊偶極子數量的關系

設不帶有GPU的處理器有m臺，帶有GPU的處理器有n臺。

（1）把總計算量按照CPU和GPU的計算能力分二大類，以保證異構系統中GPU的能力得到充分發揮。然后，在分類內部按照處理器能力分配計算量。步驟如下：

1）求出需要計算偶極子的總量，用Q總計算量表示。

2）對于各個CPU和GPU種類，用實驗方法獲得每類的計算能力（單位時間內可以計算的偶極子數量），CPUi計算能力和GPUi計算能力 表示分別表示第i個CPU和GPU的計算能力。按照降序排序，i＝1，表示計算能力最強。

3）把異構系統中的計算機處理能力，按照CPU和GPU分為兩類，對兩類處理器能力分別進行統計，CPU總計算能力和GPU總計算能力分別代表CPU處理器和GPU處理器總的計算能力。

4）把計算總量為Q總計算量分配給CPU處理器和GPU處理器，用CPU總計算量和GPU總計算量表示CPU處理器和GPU處理器所承擔的總計算量。則GPU總計算量：CPU總計算量＝GPU總計算能力：CPU總計算能力，GPU總計算量＋CPU總計算量＝Q總計算量

5）按照計算能力的大小計算每個處理器所分配的計算量，用CPUi計算量表示分配給第i類CPUi的計算量，i＝1，2…m。GPUi計算量表示分配給第i類GPUi的計算量，i＝1，2…n。則CPUi計算量＝CPU總計算量＊CPUi計算能力／CPU總計算能力。同樣，GPUi計算量＝GPU總計算量＊GPUi計算能力／GPU總計算量。

6）對偶極子計算任務按照計算量進行降序排序，并且劃分成Q1，Q2…QN，其中N＝n（GPU處理器數）＋m（CPU處理器數）。使得 Qi＝GPUi計算量，i＝1，2…n，Qn＋I＝CPUi計算量，I＝1，2…m。

（2）對于GPU1處理器，按照偶極子計算量的降序逐個模塊提交，對于CPUm處理器，按照偶極子計算量的升序逐個模塊提交。對于其它的GPUi處理器（i＝2，3…n）和CPUi處理器（i＝1，2…m－1）按照偶極子計算量的中間優先，逐步向二邊交叉延伸。這樣的提交次序，便于在計算的最后階段，把剩余計算的任務向鄰居節點轉移。

（3）網絡計算的特點之一，是在運行過程中會發生增加或者減少處理器的情況，此時，按照以下步驟進行：

1）計算出還沒有提交的任務（包括在減少的處理機中被中斷并且取消的任務）的剩余總計算量，設定為Q剩余總計算量。并且按照新的系統配置情況計算出到全部任務執行完成還需要多少時間，設定為T剩余總計算時間。

2）計算出正在正常運行的各個處理器i當前的任務結束還需要多少時間，設定為T i剩余計算時間，修改處理器i的計算能力為CPUi新計算能力＝CPUi計算能力×（1－T i剩余計算時間／T剩余總計算時間），對于GPU處理器也同樣處理。

3）按照新的系統配置情況，使用調度程序計算各個處理器需要承擔的任務。由于各個處理器的剩余計算量是不同的，因此各個處理器的新計算能力也是各不相同的，原有的分類將取消，每個處理器都單獨成為一個類。

（4）當某一處理器k結束自身計算任務，同時k－1或者k＋1處理器還有多個任務等待，可以把相應的計算任務遷移到處理器k上。直至整個系統計算任務結束為止。

（5）在開始計算過程以后，需要密切注意GPUn計算速度，如果發現GPUn速度低于期望值很多時，則表示在GPUn中，計算了大量的小任務，使得GPU計算速度下降。從我們分析偶極子計算量的分布情況來看，70－80%的偶極子模塊計算時，GPU處理器的計算速度基本上保持著最佳狀態，也就是說，只要在網絡計算系統中，CPU處理器的計算量超過20%，我們所提出的調度算法是可以獲得最佳結果的。我們認為在高性能計算中，這種情況是有典型性的。當然如果在極端情況下，計算任務中包含了大量的小任務，迫使GPU處理器承擔了很多小任務，由于我們的調度算法在分配GPU處理器和CPU處理器承擔的計算量時，是按照GPU處理器在處理大任務情況下的性能。當GPU處理器需要處理大量小任務時，完成所分配的任務需要的時間將會超過預定時間，適當把分配給GPU處理器的任務遷移到CPU處理器。在遷移的過程中，可以參照步驟（3）的方法。

5 實驗結果

實驗環境：

CPU處理器：Intel（R）Core（TM）i7CPU 920＠2.67GHz雙核。

GPU處理器：NVIDIA GeForce 9800GT。

實驗內容：

（1）采用3種不同配置，如下：

配置1：CPU處理器10臺，GPU處理器30臺。

配置2：CPU處理器20臺，GPU處理器20臺。

配置3：CPU處理器30臺，GPU處理器10臺。

（2）在每個配置上，采用二種不同算法，如下：

算法1：本調度算法。

算法2：先到先服務算法，有多個空閑處理器時，選用性能最高的處理器。

（3）輸入負荷：室上性心動過速心電仿真的4.8s偶極子數。

實驗結果：

如圖5所示，對于三證不同的配置平臺，分別使用算法1和算法2。

圖5 不同配置下的實驗結果

6 結束語

本論文在分析了心電仿真計算特點的基礎上，提出了在網絡計算系統中的一個新的調度算法，在系統的響應時間方面取得了很好的效果。

在本文心電仿真計算中，除了在并行性方面適合網絡計算的要求以外，其中的主要前提是計算工作量的可預測性，這個要求在高性能計算方面還是有一定的可行性的，例如，本調度算法可以應用到N－Body問題Fmm算法中，因為在每個分塊中近距離粒子作用力的計算量和粒子數的平方成正比，它的計算量也是可以預測的。

進一步工作，將在總結其它類型的高性能計算的算法特點的基礎上，研究動態調度算法。

［1］Shen W，Wei D，Xu W，et al.Parallelized computation for computer simulation of electrocardiograms using personal computers with multi－core CPU and general－purpose GPU［J］.Computer Methods and Programs in Biomedicine，2010，100（1）：87－96.

［2］McFarlane R，Biktasheva I V.High performance computing for the simulation of cardiac electrophysiology［C］／／Sliema：Proceedings of the Third International Conference on Software Engineering Advances，2008：13－18.

［3］Tusscher K，Panfilov A.Modelling of the ventricular conduction system［J］.Progress in Biophysics and Molecular Biology，2008，96（1－3）：152－170.

［4］Reumann M，Fitch BG，Rayshubskiy A，et al.Large scale cardiac modeling on the Blue Gene supercomputer［C］／／Vancouver，BC：Proceedings of IEEE Engineering in Medicine and Biology Society，2008：577－580.

［5］Reumann M，Fitch B G，Rayshubskiy A，et al.Strong scaling and speedup to 16，384processors in cardiac electro－mechanical simulations［C］／／Minneapolis，MN：Annual International Conference of the IEEE Engineering in Medicine and Biology Society，2009：2795－2798.

［6］Donassolo B，Legrand A，Geyer C.Non－cooperative scheduling considered harmful in collaborative volunteer computing environments cluster，cloud and grid computing（CCGrid）［C］／／Newport Beach，CA：11th IEEE／ACM International Symposium on，2011：144－153.

［7］Sato D，Xie Y，Weiss J N，et al.Acceleration of cardiac tissue simulation with graphic processing units［J］.Medical and Biological Engineering and Computing，2009，47（9）：1011－1015.

［8］Zhu H，Sun Y，Rajagopal G，et al.Facilitating arrhythmia simulation：The method of quantitative cellular automata modeling and parallel running［OL］.Biomedical Engineering Online，2004－03－29.

［9］Zhu Xin，Wei Daming，Wang Hui.Simulation of intracardial potentials with anisotropic computer heart models［C］／／Shanghai：Proceedings of the IEEE Engineering in Medicine and Biology，2005：1071－1074.

［10］Zhu Xin，Wei Daming.Computer simulation of intracardiac potential with whole－heart model［J］.International Journal of Bioinformatics Research and Applications，2007，3（1）：100－122.