基于貪婪稀疏方法的心臟磁場源重構*

邴 璐 王偉遠 王永良 蔣式勤?

1)(同濟大學電子與信息工程學院，控制科學與控制工程系，上海 201804)

2)(中國科學院上海微系統與信息技術研究所，上海 200050)

(2012年12月29日收到)

1 引言

研究心臟電活動的方法之一是用人體胸腔表面測量到的心臟磁場數據重構產生該磁場的電流源.這種方法也被稱作心臟電流源重構(cardiac source reconstruction)或心電活動的三維磁成像[1].目的是通過心臟電活動的可視圖像，無創地研究其心臟功能與診斷心臟疾病.

心臟磁場與腦磁場的源重構具有相似的理論基礎.一些重要的源重構方法，比如，空間濾波法[2,3]和稀疏分解方法[4]，最早均用于腦電/磁(EEG/MEG)研究.腦的外部磁場信號是由頭顱內特殊部位的細胞電流產生的.這些特殊部位數量不多，體積很小.然而，心臟內部的細胞電流較為復雜，電活動的傳導機理表明，它是在心臟內快速傳播的[5].因此，如何建立心臟磁場的數學模型，如何優化其逆問題求解方法，如何面向臨床應用，實現心臟電活動可視化，是當前心磁研究的重點.

理論上，通常用點狀分布的電偶極子模擬心臟磁場的激勵源.本文采用固定位置的分布源模型，將計算磁場源參數轉化為求解偶極矩的逆問題.但是，這種基于分布源模型的線性方程中，未知源的數量往往超過磁場測點的個數，理論上沒有唯一解.因此，需要施加各種約束條件，并采用優化方法求解[6,7].

稀疏信號重構是一種通過有限數據求局部能量的優化方法[8].被廣泛應用于各種信號估計，信號分類，信號處理，圖像壓縮和生物醫學成像.簡單說，稀疏信號(sparse signal)是指一個信號向量中只有有限個非零元素.稀疏分解方法(sparse decomposition)可以將一組信號在給定變換域(過完備字典)中，通過其中向量的線性組合，得到信號最稀疏的表達.

1995年，Gorodnitsky等將一種聚焦欠定系統求解器 (focalunderdetermined system solver，FOCUSS)方法引入腦磁的研究[4].他們在加權最小范數解的基礎上，提出了一種逐步加權的迭代優化算法.即利用上一步估計的源強度，作為權矩陣的對角線元素，從而在給定的分布源中挑選出所謂的局部能量最大解或稱之為稀疏解.2005年，Liu等在此基礎上提出了具有標準收縮空間的LORETA-FOCUSS方法.他們將Standard LORETA(sLORETA)的平滑估計結果作為初值，再采用FOCUSS的迭代權最小化范數，自動調整源搜索空間，進行腦電的時空源定位[9].該方法在復雜源和無噪聲的情況下獲得了很好的定位精度.2007年，Xu等，提出了腦電源重構的模迭代稀疏解[10]，這種方法對于分散的稀疏源有較好的定位能力.同年，Manoharan等研究了腦磁場源重構的多尺度FOCUSS方法，討論了動態初值選取對源重構的影響[11].最近，我們也研究了FOCUSS方法，并提出了一種加入稀疏性模約束條件的方法，對FOCUSS算法作了改進[12].仿真結果表明，它可以有效地收斂于真實的稀疏源的活動位置，實現磁場信號的電流源重構.

目前，在腦源重構中所用的稀疏方法主要有“凸”優化方法和貪婪方法，其各有優缺點.“凸”優化方法的基礎是l1范數最小化約束[8].貪婪方法包括匹配追蹤(MP)[13]或正交匹配追蹤(OMP)[14]等算法.凸優化算法的源重構精度較高，貪婪算法的計算復雜度較低[15].但尚未發現該方法被用于求解心磁逆問題或用于心磁源重構.

本文的工作是用貪婪向量(或稱作原子)選擇的思想和改進后的稀疏分解算法進行心臟磁場信號的源重構.目的是提高源重構的計算速度和精度，通過重構稀疏源在三維空間中的移動軌跡，使心臟電活動可視化，并揭示心臟除極和復極時電活動的傳導.

2 方法

分布電流源產生的磁場可用如下線性方程表示：

其中，b是一個M×1的向量，表示用超導量子干涉器記錄的M個陣列測點上的磁場信號.x是N×1的向量，為待定的N個分布電流源的強度.A是一M×N的導聯矩陣，表示分布電流源與磁場測點的關系.這里A的每一列向量被稱作原子.通常，測量點個數遠遠低于源搜索空間的維數，所以，M?N，方程(1)是欠定的.本文通過一種快速貪婪方法得到了該問題的稀疏解，即采用改進的正交匹配追蹤方法，快速搜索對應源空間中的局部能量最大解.

2.1 正交匹配追蹤方法

正交匹配追蹤算法(orthogonalmatching pursuit，OMP)[14]，是一種貪婪的迭代算法，是對匹配追蹤算法(matching pursuit，MP)[13]的改進.在OMP算法中，每步迭代都對所選的原子進行Gram-Schmidt正交化處理.即將信號投影到正交原子構成的空間，并將信號投影后得到的分量和殘余分量再次用同樣的方法分解.經M次分解后，該信號被分解為M個原子的線性組合.這樣可以對信號的逼近始終保持一個較高的收斂速度，從而克服了匹配追蹤算法收斂較慢的問題[15].

2.2 近似正交化過程

對于大規模的線性優化問題，上述內積計算非常耗時.我們采用梯度追蹤方法替代精確的正交化[16]，即采用近似正交映射(方向優化)的方法.OMP算法中，更新電流源強度向量可簡單表示為

可以驗證，本文采用的近似共軛梯度更新方法在計算復雜度方面優于梯度更新，而且，其源重構精度與OMP方法差別很小.也就是說，梯度追蹤算法中，不需要保證余量rn與被選擇的A的列向量的正交性.這樣，算法每次迭代時可以重新在A中選擇原子.

2.3 不同的原子選擇過程

在一般的OMP算法或標準的梯度追蹤算法中，每次迭代只選擇一個原子，需要多次迭代后才能知道是否出現非零原子.文獻[17]中的Stagewise OMP(StOMP)，通過給定一個閾值λ，使得一次迭代可以選擇一個以上的原子：

即選擇矩陣A中所有內積大于該閾值的原子：

(8)式中，參數t的選擇會影響算法的性能，但是沒有具體的選擇方法，甚至當所有內積遠低于閾值時，算法可能失效.

本文中，由于待重構的電流源的稀疏度是未知的，所以，我們在A中選擇所有使上述內積最大的原子

其中，γ為給定參數.被選擇的原子集合定義為

在MCG電流源的稀疏度未知的情況下，這種原子選擇方式也能保證比較精確的源重構.

3 心磁數據的電流源重構

本文用4通道超導量子干涉器(SQUID)測量到的一組正常人的心臟磁場強度數據，舉例說明了其QRS段(100 ms)和ST-T段(100 ms)源重構的計算結果.圖1是在人體胸腔表面20 cm×20 cm的測量平面上，經過9次測量，以及數據同步和信號處理后，得到的36通道的單周期心磁數據.分別從該數據中截取QRS段和ST-T段100 ms心磁數據.

圖1 一組正常人的心臟磁場強度數據(QRS段(234—333 ms);ST-T段(457—556 ms))

3.1 用快速貪婪稀疏方法的計算結果

假設在測量平面上有21×21個等效的分布電流源.人體胸腔表面向下沿Z方向從3 cm到12.5 cm分為20層，使整個源空間的分辨率為1.0×1.0×0.5 cm.用快速貪婪稀疏方法計算得到的電流源分布如圖2和圖4所示.圖中分別顯示了QRS段和ST-T段中幾個時刻的磁場等高線圖與稀疏源的對應關系.其中，源深度用數字表示，單位為cm.圖3和圖5分別是強度為65%以上的電流源(稱之為主導電流源)沿X，Y，Z方向隨時間變化的運動軌跡和強度變化曲線.將主導電流源顯示在MRI圖像的X-Y平面，Z-Y平面，X-Z平面中以便說明主導電流源的相對位置.圖3中，QRS段用實線劃分為begin-Q，Q-R，R-S和S-end四部分，以便觀察該時間段內心臟的電活動，圖5中標出了T波段的心臟電活動.

由圖2和圖4可見，重構的電流源多數分布在零磁場線附近，尤其是紅色的最強的那個等效電流源.如圖 3(a)，(b)，(c)和圖 5(a)，(b)，(c)所示，稀疏電流源的空間移動軌跡基本位于源空間范圍內，QRS段電流源的深度變化在4—10 cm左右，ST-T段的變化在 6—9 cm 左右.由圖 3(d)，(e)，(f)和圖 5(d)，(e)，(f)可見，強度為65%以上的稀疏電流源的移動軌跡基本位于心臟中.由于QRS段和ST-T段兩端的測量信號的信噪較低，重構的稀疏源中有強度較低的偽電流源.此外，由于X-Y平面的邊界上存在信號插值引起的誤差，所以，圖2和圖4中邊界上有重構的偽源[18].

圖2 QRS段中各時刻的重構電流源分布與磁場等高線圖

圖3 QRS段主導稀疏電流源的空間移動軌跡 (a)沿X方向位置軌跡;(b)沿Y方向位置軌跡;(c)沿Z方向位置軌跡(深度);(d)MRI圖像上X-Y平面的電流源軌跡;(e)MRI圖像上Z-Y平面的電流源軌跡;(f)MRI圖像上X-Z平面電流源軌跡

3.2 源重構的精度

通常用兩個指標評價源重構精度[5].指標1是測量磁場信號與重構磁場信號差的平均值.定義每個時間間隔中重構電流源產生的磁場信號b′(t)與測量信號b(t)差的平均值為

指標2是測量磁場信號與重構磁場信號差的絕對值與測量信號最大值的比率.定義測量磁場和由重構電流源產生的磁場信號差的絕對值與測量信號最大值bmax的比率為

圖6和圖8中給出了這兩個指標的分析結果.圖7和圖9分別是QRS段和ST-T段中每間隔4 s的測量值和重構的磁場等高線圖.可以看出，重構的磁場等高線圖與測量的等高線圖相似.圖10是它們的擬合優度(goodness of fit，GOF)隨時間變化曲線，平均擬合優度分別為99.36%和99.78%.

圖4 ST-T段各時刻的重構電流源分布與磁場等高線圖

圖5 ST-T段主導稀疏電流源的空間移動軌跡 (a)沿X方向位置軌跡;(b)沿Y方向位置軌跡;(c)沿Z方向位置軌跡(深度);(d)MRI圖像上X-Y平面的電流源軌跡;(e)MRI圖像上Z-Y平面的電流源軌跡;(f)MRI圖像上X-Z平面電流源軌跡

圖6 QRS段源重構精度的分析 (a)指標1;(b)指標2

圖7 QRS段測量和重構的磁場等高線圖 (a)測量的;(b)重構的

圖8 ST-T段源重構精度的分析 (a)指標1;(b)指標2

圖9 ST-T段測量和重構的磁場等高線圖 (a)測量的;(b)重構的

4 討論

4.1 運算速度比較

我們比較了兩種源重構算法，即本文給出的快速貪婪稀疏方法和之前的改進FOCUSS(MFOCUSS)方法[12]，計算一組心磁數據所需要的運算時間，如圖11所示.雖然后者導聯矩陣A的列數大約是前者的5倍，但是，用快速貪婪稀疏方法比MFOCUSS方法的運算速度快得多，所用時間約為MFOCUSS的1/50.如果計算同樣規模的導聯矩陣A，MFOCUSS算法可能失效.

圖10 QRS段和ST-T段測量和重構磁場等高線圖的擬合優度 (a)QRS;(b)ST-T

圖11 兩種算法的運算時間比較 (a)MFOCUSS算法;(b)快速貪婪稀疏算法

FOCUSS也是一種用來解決線性逆問題的優化方法[4].它通過稀疏約束優化一個lp(p≤1)的正則化代價函數，是一種使得解的能量局部化的迭代過程.該算法需要選擇正則化參數及給出停止條件.對于特定的應用對象來說這是不容易的.然而，本文的貪婪方法的參數可以通過信號稀疏度和計算性能要求綜合考慮[19].在選擇停止標準時，兩種算法完全不同.貪婪算法在每次迭代中選擇新原子，可控制稀疏度，減少誤差，直到滿足停止條件.然而，FOCUSS算法中，停止條件僅與代價函數的最小值有關.

4.2 重構電流源的時空分布特點

從圖2和圖4中可見，QRS段和ST-T段每個時刻的心磁圖中基本都有2個較強的電流偶極子，它們方向不同.從圖3中QRS段電流源沿X，Y和Z方向隨時間變化的曲線可見，Q點(大約254 ms)到R峰(大約284 ms)約30 ms.在除極過程中，等效電流源從被測量對象胸腔的右上方逐漸向左下方移動，與左心室興奮由內向外擴散相對應.然后，快速地從左下方逐漸返回右上方，與興奮由心尖向心底擴散，最晚興奮處在心臟的基底部相對應[20,21].

該正常人QRS段等效電流源深度的變化具有明顯特征，即先是從測量平面以下約4 cm處(對應Q點，約254 ms)沿Z方向逐漸向下，達到約10 cm(對應S點，約304 ms)后，迅速返回約6 cm處.從圖5中ST-T段電流源沿X和Y方向隨時間變化的曲線可見，大約到T峰時刻等效電流源開始向下移動.該正常人ST-T段等效電流源的深度變化在6 cm到9 cm左右.

4.3 研究限制

本研究中，所用的導聯矩陣A尚不能精確的反映不同個體的胸腔和心臟組織情況，有待進一步研究與比較.我們用邊界元方法構建的體電導模型進行了數值仿真，結果說明，心臟和軀干邊界的影響會引起源定位誤差，體電導率的大小也會影響重構電流源的精度.可能引起的源重構誤差有待深入研究.

5 結論

我們用快速貪婪稀疏方法得到了心臟磁場的電流源分布.與我們前期所用的改進FOCUSS方法相比，在源估計精度大致相當的情況下，這種新方法的計算速度有了明顯的提高.我們根據定性分析QRS段和ST-T段中稀疏源與MRI圖的對應關系，稀疏源與磁場等高線圖的對應關系，以及根據測量心臟磁場和用重構電流源產生的磁場等高線圖的擬合優度和兩個源定位精度指標的定量分析，驗證了這種方法及相關結果的有效性.

因為心臟磁場測量數據的質量直接影響到源重構的結果，所以，在評判一組心磁數據是否可利用時，我們的經驗是，先用多種方法檢驗該數據是否能夠得到相似的源分布結果，心磁圖對應的梯度圖是否存在邊界誤差等.

本文討論了正常人QRS段和ST-T段稀疏電流源的時空分布特點，分析了心臟除極和復極時最強稀疏源移動軌跡的相關信息.但是對于典型心臟病人的心臟電活動的研究和統計分析，尚需一定數量的臨床數據和病人資料.

本文作者感謝張懿教授，謝曉明教授及其團隊，王樂民教授，羅明教授，權微微醫生和李文生醫生為本研究提供的幫助.感謝中科院上海微系統所和上海瑞金醫院對本研究的大力支持.

[1]He B 2003 Modeling and Imaging of BioelectricalActivity：Principles and applications(Minnesota：Kluwer Academic/Plenum Press)

[2]Gutierrez D，NehoraiA，Dogandzic A 2006 IEEE Transactions on BiomedicalEngineering 53 840

[3]Moiseev A，Gaspar J M，Schneider J A，Herdman A T 2011 Neuroimage 58 481

[4]Gorodnitsky I F，George J S，RaoB D 1995 Electroencephalography and clinicalNeurophysiology 95 231

[5]Melis M D，Tanaka K，Uchikawa Y 2010 IEEE Transactions on magnetics 38 1203

[6]MichelC M，Murray M 2004 ClinicalNeurophysiology 115 2195

[7]Davis G 1994 Ph.D.Dissertation(U.S.A：New York University)

[8]Chen S S，DonohoD L，Saunders M A 1998 SIAM Journalof Scientific Computing 20 33

[9]Liu H S，Schimpf P H，Dong G Y，GaoX R，Yang F S，GaoS K 2005 IEEE Transactions on BiomedicalEngineering 52 1681

[10]Xu P，Tian Y，Chen H F，YaoD Z 2007 IEEE Transactions on BiomedicalEngineering 54 400

[11]Manoharana A，Morana J E，Bowyera S M，Masona K M，Tepleya N，Smitha B J，Barkleya G L 2007 InternationalCongress Series.Eng.1300 665

[12]Bing L，Jiang S Q，2012 InternationalConference on System Simulation Shanghai，China，April4—6，2012 p298

[13]Mallat S，Zhang Z 1993 IEEE Transactions on SignalProcessing 41 3397

[14]PatiY C，Rezaiifar R，Krishnaprasad S 1993 27th Annum Asilomar Conference on Signals Systems and Computers California，U.S.A，November 1—3，1993 p101

[15]Mallat S，Davis G，Zhang Z 1994 SPIE Journalof OpticalEngineering 33 2183

[16]Blumensath T，Davies M 2008 IEEE Transactions on SignalProcessing 56 2370

[17]DonohoD，Tsaig Y，DroriI，Starck J 2006 Tech.Rep.(U.S.A：Stanford University Press)

[18]Zhu J J，Jiang S Q 2012 InternationalConference on System Simulation Shanghai，China，April4—6，2012 p302

[19]Davies M E，Blumensath T 2008 Proceedings of the third InternationalSymposium on Communications,Controland SignalProcessing Malta，March 12—14，2008 p774

[20]Wang C S 2007 Cardiac conduction system(Beijing：Qinghua University Press)p84(in Chinese)[王成山2007心臟傳導系統基礎與臨床(北京：清華大學出版社)第84頁]

[21]Yasunaga H，Yasunaga H，Harumizu S，NaomiI，SeikoK，Tohru S，Masayasu H 2002 Journalof Electrocardiology 35 p105