基于貪婪稀疏方法的心臟磁場(chǎng)源重構(gòu)*

邴 璐 王偉遠(yuǎn) 王永良 蔣式勤?

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

2)(中國科學(xué)院上海微系統(tǒng)與信息技術(shù)研究所，上海 200050)

(2012年12月29日收到)

1 引言

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

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

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

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

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

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

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

2 方法

分布電流源產(chǎn)生的磁場(chǎng)可用如下線性方程表示：

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

2.1 正交匹配追蹤方法

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

2.2 近似正交化過程

對(duì)于大規(guī)模的線性優(yōu)化問題，上述內(nèi)積計(jì)算非常耗時(shí).我們采用梯度追蹤方法替代精確的正交化[16]，即采用近似正交映射(方向優(yōu)化)的方法.OMP算法中，更新電流源強(qiáng)度向量可簡(jiǎn)單表示為

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

2.3 不同的原子選擇過程

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

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

(8)式中，參數(shù)t的選擇會(huì)影響算法的性能，但是沒有具體的選擇方法，甚至當(dāng)所有內(nèi)積遠(yuǎn)低于閾值時(shí)，算法可能失效.

本文中，由于待重構(gòu)的電流源的稀疏度是未知的，所以，我們?cè)贏中選擇所有使上述內(nèi)積最大的原子

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

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

3 心磁數(shù)據(jù)的電流源重構(gòu)

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

圖1 一組正常人的心臟磁場(chǎng)強(qiáng)度數(shù)據(jù)(QRS段(234—333 ms);ST-T段(457—556 ms))

3.1 用快速貪婪稀疏方法的計(jì)算結(jié)果

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

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

圖2 QRS段中各時(shí)刻的重構(gòu)電流源分布與磁場(chǎng)等高線圖

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

3.2 源重構(gòu)的精度

通常用兩個(gè)指標(biāo)評(píng)價(jià)源重構(gòu)精度[5].指標(biāo)1是測(cè)量磁場(chǎng)信號(hào)與重構(gòu)磁場(chǎng)信號(hào)差的平均值.定義每個(gè)時(shí)間間隔中重構(gòu)電流源產(chǎn)生的磁場(chǎng)信號(hào)b′(t)與測(cè)量信號(hào)b(t)差的平均值為

指標(biāo)2是測(cè)量磁場(chǎng)信號(hào)與重構(gòu)磁場(chǎng)信號(hào)差的絕對(duì)值與測(cè)量信號(hào)最大值的比率.定義測(cè)量磁場(chǎng)和由重構(gòu)電流源產(chǎn)生的磁場(chǎng)信號(hào)差的絕對(duì)值與測(cè)量信號(hào)最大值bmax的比率為

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

圖4 ST-T段各時(shí)刻的重構(gòu)電流源分布與磁場(chǎng)等高線圖

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

圖6 QRS段源重構(gòu)精度的分析 (a)指標(biāo)1;(b)指標(biāo)2

圖7 QRS段測(cè)量和重構(gòu)的磁場(chǎng)等高線圖 (a)測(cè)量的;(b)重構(gòu)的

圖8 ST-T段源重構(gòu)精度的分析 (a)指標(biāo)1;(b)指標(biāo)2

圖9 ST-T段測(cè)量和重構(gòu)的磁場(chǎng)等高線圖 (a)測(cè)量的;(b)重構(gòu)的

4 討論

4.1 運(yùn)算速度比較

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

圖10 QRS段和ST-T段測(cè)量和重構(gòu)磁場(chǎng)等高線圖的擬合優(yōu)度 (a)QRS;(b)ST-T

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

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

4.2 重構(gòu)電流源的時(shí)空分布特點(diǎn)

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

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

4.3 研究限制

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

5 結(jié)論

我們用快速貪婪稀疏方法得到了心臟磁場(chǎng)的電流源分布.與我們前期所用的改進(jìn)FOCUSS方法相比，在源估計(jì)精度大致相當(dāng)?shù)那闆r下，這種新方法的計(jì)算速度有了明顯的提高.我們根據(jù)定性分析QRS段和ST-T段中稀疏源與MRI圖的對(duì)應(yīng)關(guān)系，稀疏源與磁場(chǎng)等高線圖的對(duì)應(yīng)關(guān)系，以及根據(jù)測(cè)量心臟磁場(chǎng)和用重構(gòu)電流源產(chǎn)生的磁場(chǎng)等高線圖的擬合優(yōu)度和兩個(gè)源定位精度指標(biāo)的定量分析，驗(yàn)證了這種方法及相關(guān)結(jié)果的有效性.

因?yàn)樾呐K磁場(chǎng)測(cè)量數(shù)據(jù)的質(zhì)量直接影響到源重構(gòu)的結(jié)果，所以，在評(píng)判一組心磁數(shù)據(jù)是否可利用時(shí)，我們的經(jīng)驗(yàn)是，先用多種方法檢驗(yàn)該數(shù)據(jù)是否能夠得到相似的源分布結(jié)果，心磁圖對(duì)應(yīng)的梯度圖是否存在邊界誤差等.

本文討論了正常人QRS段和ST-T段稀疏電流源的時(shí)空分布特點(diǎn)，分析了心臟除極和復(fù)極時(shí)最強(qiáng)稀疏源移動(dòng)軌跡的相關(guān)信息.但是對(duì)于典型心臟病人的心臟電活動(dòng)的研究和統(tǒng)計(jì)分析，尚需一定數(shù)量的臨床數(shù)據(jù)和病人資料.

本文作者感謝張懿教授，謝曉明教授及其團(tuán)隊(duì)，王樂民教授，羅明教授，權(quán)微微醫(yī)生和李文生醫(yī)生為本研究提供的幫助.感謝中科院上海微系統(tǒng)所和上海瑞金醫(yī)院對(duì)本研究的大力支持.

[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心臟傳導(dǎo)系統(tǒng)基礎(chǔ)與臨床(北京：清華大學(xué)出版社)第84頁]

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