基于網(wǎng)絡(luò)計算的心電仿真計算研究

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

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

0 引 言

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

本文提出采用由PC機中的CPU和GPU組成異構(gòu)網(wǎng)絡(luò)計算的方案，利用空閑的PC資源，實現(xiàn)心電仿真計算，并且提出了面向短的響應(yīng)時間要求的調(diào)度算法，取得了很好的效果。

1 相關(guān)工作

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

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

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

心臟的計算機模型一般是通過對心臟斷面圖像離散化并進行空間插補，本論文所采用實驗?zāi)Ｐ蜑槿毡疚菏先呐K模型［9－10］，該模型是由一個56＊56＊90的斜交坐標(biāo)系構(gòu)成。模型上部為心房，下部為心室，模型細(xì)胞在空間上是等距的，直徑為1.5ms。模型的實際細(xì)胞數(shù)據(jù)為5萬個，通過使模型細(xì)胞具有方向性和具有和實際細(xì)胞相近的電生理特性，在解剖學(xué)和電生理學(xué)兩個方面對模型作進一步的修飾使得模型更加接近實際的心臟。

2.1 心電仿真算法

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

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

電勢計算方程如下

圖1 人體軀干

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

2.2 心電仿真計算流程圖

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

圖2 心電仿真計算流程

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

3 異構(gòu)網(wǎng)絡(luò)計算的結(jié)構(gòu)

3.1 心電仿真計算在計算結(jié)構(gòu)上的特點分析

心電仿真計算在計算結(jié)構(gòu)上有以下的特點：高度并行，占總仿真時間99.7%的計算是高度并行的，計算模塊之間的相互通信量非常少，各個模塊完全可以獨立計算。非常適合利用網(wǎng)絡(luò)計算。各個模塊的計算量是可以預(yù)測的。為在網(wǎng)絡(luò)計算中實現(xiàn)定量調(diào)度提供了條件。

3.2 異構(gòu)網(wǎng)絡(luò)結(jié)構(gòu)介紹

異構(gòu)網(wǎng)絡(luò)結(jié)構(gòu)如圖3所示。

圖3 異構(gòu)網(wǎng)絡(luò)結(jié)構(gòu)

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

4 調(diào)度算法

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

調(diào)度算法的執(zhí)行步驟如下：

圖4 Geforce440GPU的處理速度與處理模塊偶極子數(shù)量的關(guān)系

設(shè)不帶有GPU的處理器有m臺，帶有GPU的處理器有n臺。

（1）把總計算量按照CPU和GPU的計算能力分二大類，以保證異構(gòu)系統(tǒng)中GPU的能力得到充分發(fā)揮。然后，在分類內(nèi)部按照處理器能力分配計算量。步驟如下：

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

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

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

4）把計算總量為Q總計算量分配給CPU處理器和GPU處理器，用CPU總計算量和GPU總計算量表示CPU處理器和GPU處理器所承擔(dān)的總計算量。則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）對偶極子計算任務(wù)按照計算量進行降序排序，并且劃分成Q1，Q2…QN，其中N＝n（GPU處理器數(shù)）＋m（CPU處理器數(shù)）。使得 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）按照偶極子計算量的中間優(yōu)先，逐步向二邊交叉延伸。這樣的提交次序，便于在計算的最后階段，把剩余計算的任務(wù)向鄰居節(jié)點轉(zhuǎn)移。

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

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

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

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

（4）當(dāng)某一處理器k結(jié)束自身計算任務(wù)，同時k－1或者k＋1處理器還有多個任務(wù)等待，可以把相應(yīng)的計算任務(wù)遷移到處理器k上。直至整個系統(tǒng)計算任務(wù)結(jié)束為止。

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

5 實驗結(jié)果

實驗環(huán)境：

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

GPU處理器：NVIDIA GeForce 9800GT。

實驗內(nèi)容：

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

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

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

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

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

算法1：本調(diào)度算法。

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

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

實驗結(jié)果：

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

圖5 不同配置下的實驗結(jié)果

6 結(jié)束語

本論文在分析了心電仿真計算特點的基礎(chǔ)上，提出了在網(wǎng)絡(luò)計算系統(tǒng)中的一個新的調(diào)度算法，在系統(tǒng)的響應(yīng)時間方面取得了很好的效果。

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

進一步工作，將在總結(jié)其它類型的高性能計算的算法特點的基礎(chǔ)上，研究動態(tài)調(diào)度算法。

［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，F(xiàn)itch 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，F(xiàn)itch 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－h(huán)eart model［J］.International Journal of Bioinformatics Research and Applications，2007，3（1）：100－122.