P波間期的心臟電流源重建及電活動磁成像*

周大方 蔣式勤? 趙晨 Peter van Leeuwen

1)(同濟大學電子與信息工程學院,上海 201804)

2)(維藤/黑爾德克大學健康學院,德國維藤 D-58448)

1 引 言

1970年,Cohen等[1]利用超導量子干涉儀(superconducting quantum interference device,SQUID)在人體胸腔表面無接觸、無創地測量到了心臟磁場圖(magnetocardiogram,MCG).心臟磁場十分微弱(約10–12Tesla),遠低于地球磁場(約10–6Tesla).SQUID的優點是[2,3]: 心磁測量無需外加激勵,無需接觸式電極,受試個體完全無損,操作方便.尤其是測量磁信號的時間分辨率高(≦1ms).目前,SQUID測量MCG信號的多通道系統已從早期的4通道發展到300通道.利用心磁圖無創獲取心臟功能信息的研究也取得了實質的進展[4?7].雖然心磁圖的研究比心電圖(electrocardiogram,ECG)推遲了80多年,但是,心臟電/磁信號同生共源,心磁圖研究可以作為心電圖的一種補充.MCG對心臟的切向電流(tangential current)和再入電流(reentrant current)較為敏感,可提供心臟的功能信息[3,8].

重建磁場電流源需要求解場–源逆問題[9,10].近20年腦磁(magnetoencephalogram,MEG)[11?13]和腦電(electroencephalogram,EEG)[14,15]的研究中發展了一種基于空間濾波器的重建分布電流源的方法.該方法給研究對象劃分網格,并在每個網格交點上構造一個空間濾波器.然后,重建濾波器位置上電流源的偶極矩,獲取分布源空間譜的信息.研究表明,重建分布源的精度與空間濾波器的輸出噪聲有關[16,17].雖然心磁測量是在屏蔽室內完成的,但是,在信號測量、數據處理和心臟電流源重建的環節中,仍需采用抑制噪聲的方法[18?22],尤其是當P波間期信噪比(signal-to-noise ratio,SNR)較低(≤ 30 dB)的時候.

最小方差波束成形(minimum variance beamforming,MVB)是一種改進的空間濾波方法[11].為了提高重建電流源的精度,MVB采用了空間濾波器的輸出總功率最小化和噪聲空間譜強度歸一化,以及在線調整權矩陣中測量磁信號的自適應技術[23,24].因此,MVB重建電流源的精度比非自適應的最小范數空間濾波(minimum norm spatial filtering,MN)方法高[25,26].Sekihara等[11,13]還提出了一種可降低MVB空間濾波器輸出噪聲功率的信號子空間投影方法,并將其用于腦功能的磁源定位.腦磁源定位的方法需要假設腦的有效電流源數目少于磁場測量通道的數目,以及綜合導聯場矩陣(composite lead-field matrix)是列滿秩的.然而,由于心臟結構復雜且有動態變化,一般情況下這兩個條件不能同時滿足[17].為此,我們曾經提出了一種抑制空間濾波器輸出噪聲功率增益(suppressing spatial filter output noise-power gain,SONG)的波束成形,可用于R波峰時刻心臟電活動成像[27].

對心臟P波間期電活動無創成像的難點是P波心臟磁信號比R波峰弱,信噪比較低.針對這一問題,本文提出了一種基于空間濾波器的可提高分布源空間譜估計強度對比度(improving intensity contrast of distributed source spatial spectrum estimation,IIC)的波束成形.因為分布源模型中導聯矩陣表示傳感器陣列對一個源所在位置的測量靈敏度,所以,IIC方法在空間濾波器的加權矩陣中引入導聯場矩陣,可以使濾波器輸出估計對磁場電流源及其分布比較敏感,從而改進分布源強度的對比度.再通過設置源強度閾值,提取每個時刻重建源中偶極矩強度極大的電流源,就可以消除其他位置上相對弱的重建分布源與超出心臟范圍的偽源,提高P波間期電流源重建的精度.文中給出了IIC方法單源重建的理論分析與仿真試驗,并比較了該方法與MVB和SONG電流源重建方法的性能.文中還用2個健康人的P波段61通道心臟磁場測量數據重建心臟電流源,并分析了IIC方法的成像結果.

2 方 法

2.1 磁場分布源的模型及其空間譜估計

假設SQUID系統在人體胸腔表面測量到的心臟磁場由n個呈網格分布的電流源產生.測量面上第k通道的心磁測量信號用標量bk(t)表示,則t時刻c通道陣列信號的列向量可用b(t)=[b1(t),b2(t),...,bc(t)]T表示.磁場與分布源的關系模型可用如下線性方程表示[13,16]:

心臟電活動是一個隨機過程[3,17],故文中假設電流源的偶極矩是一個隨機變量,簡記為個分布源的空間譜可表示為

本文利用空間濾波器方法重建產生磁場的電流源.已知空間濾波器的輸入是測量信號向量輸出是需要計算的電流源偶極矩估計它們的關系可用線性加權運算表示[13,24]:

2.2 IIC波束成形方法

本文在MVB[12]和SONG[27]的基礎上,提出一種IIC波束成形方法.目的是通過提高分布源強度的對比度,減小P波間期信號噪聲的影響.

2.2.1 濾波器的加權矩陣與比較

MVB和SONG的濾波權矩陣分別為[12,27]:

為了改進分布源強度的對比度,本文在此基礎上,構造了IIC的濾波權矩陣:

根據實對稱矩陣有關定理[28,29]可以證明,(5)式中的實對稱陣U是一個“特征值不大于1,矩陣的跡小于其階數”的低跡的半正定陣.且有

如果空間濾波器的輸入v是高斯白噪聲,則分布源空間譜的強度估計由(7)式可見,當空間濾波器的輸入v是高斯白噪聲時,IIC方法降低了空間濾波器的輸出噪聲功率增益和輸出噪聲功率且降噪性能比SONG和MVB好.

2.2.2 空間譜估計強度對比度的分析

由于多源空間譜估計分析的困難,本文分析了用IIC方法的單電流源空間譜估計.首先定義波束成形方法的點擴散函數為任意空間位置上對電流源強度估計的歸一化.假設任意位置上分布源和上給定單電流源s的強度估計分別為并定義源強度的對比度與源強度的對比度成反比.

和

建立科學有效的評價機制。打破傳統的把考試分數作為唯一標準來衡量學生的教學評價方法，采用非標準化考試。從注重學習知識的考評向運用知識的考量轉變，主要考量學生運用知識分析解決問題的能力；從注重最終結果的考核向學生成長過程考核，主要考量課堂參與度、線上學習次數和時長、團隊匯報展示、團隊作業、隨堂實驗等，同時客觀記錄并科學評價學生成長歷程，學生的創新實驗、技術研發、發表論文、研究課題、獲得專利、競賽成績和自主創業等都要記錄在學生成長檔案；由教師單一評價主體向自評、互評、他評和教師評價等多元評價主體轉變，進行綜合評價。

可見,這3種波束成形方法中,IIC的單源空間譜估計強度的對比度最高,即

2.2.3 局部強度極大電流源提取

本文在改進空間譜估計對比度的基礎上,采用了設置源強度閾值,提取每個時刻重建分布源中局部強度極大電流源的方法.這是因為在重建磁場分布源時,需要消除環境噪聲和計算誤差引起的大量分布的重建弱源與重建的偽源(可能超出心臟范圍).本文假設源空間中每個分布源的最小鄰域中包括26個其他分布源,圖1中它們的位置分別用紅色和藍色表示.并用每個時刻空間譜估計的最大強度作為提取局部極大源的閾值的基數,其中比例系數可根據經驗確定.這樣,改進源強度對比度后,可以通過設定閾值,消除部分重建的偽源與大量較弱的重建分布源.最后,利用那些局部強度極大的電流源來研究心臟的電活動.

圖1 分布源位置及指定的26個鄰域位置,分別用紅色和藍色表示Fig.1.The position of the source and 26 specified adjacent positions are shown in red and blue colors,respectively.

3 心臟電活動的仿真與成像

電活動的仿真與成像是指用仿真的磁場數據重建電流源,并記錄電流源在空間中的位移,給出源位移隨時間變化的圖像.目的是檢驗用局部強度極大電流源近似磁場電流源,以及重建電流源在空間中的移動模擬電活動成像的效果.本文將根據重建電流源中偽源的個數、重建源的誤差、以及成像結果中電活動的方向,分析比較IIC,SONG,MVB以及一種非空間濾波算法即信賴域反射(trust region reflective,TRR)方法的仿真與成像的結果.TRR是一種單電流源重建方法[30,31],它通過最小化磁場殘差的范數和迭代更新,求解等效的單電流源[32].

3.1 仿真數據

假設包含心臟的人體軀干是沿坐標系Z軸水平分層的導體(horizontally layered conductor,HLC),心臟–軀干模型如圖2(a)所示.層間電導率之差與該處靜電勢的乘積被稱作二次電流密度,其方向與沿Z方向測量的心臟磁場平行[17].根據畢奧薩法爾定律,導體中二次電流密度產生的磁場可以忽略[10].對圖2(a)中包含心臟在內的導體

劃分網格,由(1)式可得,間距為1 cm的網格交點將構成一個n=12168的分布源空間.

3.1.1 磁場導聯場矩陣

參照Magnes 1300C生物磁強計系統[2],假設仿真的1 kHz磁場數據是在胸腔表面用61通道環狀陣列傳感器沿Z方向測量的.如圖2(a)所示,是通道的磁信號測量位置,故(1)式中的磁場導聯矩陣可計算如下[17]:

圖2 (a)水平分層導體.其中黑色直線表示分層軀干的邊界,黃色橢圓體表示心臟;(b),(c)測量通道1和61的磁場導聯,(ζ=X,Y,Z)曲線.其中紅、綠和藍色分別表示X,Y和Z方向導聯的曲線Fig.2.(a)The horizontally layered conductor(HLC),where the black line indicates the boundary of HLC,and the yellow spheroid the heart;(b),(c)The X,Y and Z lead-field based plots of channels 1 and 61 are expressed in red,green and blue,respectively.

3.1.2 仿真的61通道磁場數據

首先,用已知偶極矩強度的單或雙電流源產生模擬心臟的磁場.為簡單起見,假定該單或雙電流源通過了源空間中的直線參考路徑,用于仿真沿參考路徑的電活動.因為P波間期心臟電活動的強度先增后減,所以用半周期正弦函數近似描述該單或雙電流源偶極矩強度的變化[17]:

其中rp是參考路徑上某個分布源的位置.τrp是電流源出現在位置rp的時刻.ω是弧度頻率,正弦函數的周期T=2π/ω.根據P波間期人體心臟的心房和心室肌中電活動的最大速度約1 m/s[33],令直線參考路徑上相鄰位置間電流源移動的時間?t=τrp+1?τrp=0.01s,以及移動速度(moving velocity)MV=|rp+1?rp|/?t≈1m/s.其中rp與rp+1表示直線參考路徑上相鄰位置的前后.當τrp≦t<τrp+1時,電流源出現在位置rp,其偶極矩強度q(t,rp)按照半周期正弦函數先增后減.當t=τrp+1時,該電流源到達位置rp+1,其偶極矩強度q(t,rp+1)將先增加后減小至0.

然后,在給定單或雙電流源偶極矩和計算相應磁場導聯矩陣的基礎上,61通道的仿真磁場數據可用下式計算:

其中b(t)=[b1(t),b2(t),...,b61(t)]T是61通道磁場測量信號.q(t,rp)=q(t,rp)η(t,rp)表示單或雙電流源的偶極矩.Lp是單或雙電流源的磁場導聯矩陣.v為高斯白噪聲,根據文中P波間期心磁測量信號的信噪比給定.

如第2節所述,利用61通道的仿真磁場數據求解逆問題,即可得到分布源的空間譜.再從中提取局部偶極矩強度極大的電流源作為重建電流源,并給出其沿給定參考路徑位移的圖像.

3.2 仿真結果1

假設參考路徑1是一條起止點坐標為Ps=(7,5,6)cm和Pe=(–3,–5,16)cm的空間直線.如圖3(a)—圖3(e)中黑色圓圈所示,該路徑上相鄰分布源的間距為單電流源通過參考路徑上相鄰位置的時間速度圖3(b)—圖3(e)給出了SNR=20 dB時,用上述4種方法得到的重建電流源和電活動成像結果.圖中可見,用IIC每個時刻重建的電流源沿著參考路徑1上的黑圈移動.黑圈中的綠色由淺到深表示時間變化.在該參考路徑以外,IIC出現的偽源最少,SONG與TRR的結果次之.

表1中比較了上述4種方法的仿真結果中的分布電流源估計誤差,分別考慮了無噪聲以及信噪比SNR為30與20 dB的情況.表1中分別表示X,Y,Z方向分布電流源估計的根均方誤差和總誤差.無噪聲時,IIC,SONG和TRR的源估計誤差較小.MVB的誤差相對較大.SNR=30 dB時,相比其他方法,IIC可以降低誤差E.SNR=20 dB時,IIC的單源估計誤差比其他方法小.如圖3所示,IIC能夠仿真單源沿參考路徑的電活動,偽源數量最少,單源重建的精度相對其他方法要好.

3.3 仿真結果2

假設參考路徑2包括兩條起止點坐標分別為[Ps1,Pe1]=[(7,5,6),(7,–5,6)]cm和[Ps2,Pe2]=[(–3,5,16),(–3,–5,16)]cm的空間直線.如圖4(a)—圖4(e)中黑圈所示,這兩條路徑上相鄰分布源的間距均為1 cm.仿真中,令兩個電流源分別通過各自參考路徑上相鄰位置的時間速度其中一個電流源比另一個電流源提前0.005 s開始移動.圖4給出了SNR=20 dB時,4種方法的仿真結果和電活動成像.圖中可見,用IIC每個時刻重建的電流源沿參考路徑2中的黑圈移動.黑圈中的綠色由淺到深表示時間變化.在該參考路徑以外,IIC出現的偽源最少,明顯優于其他3種方法.

表2比較了上述3種空間濾波方法仿真結果中的分布電流源估計誤差,分別考慮了無噪聲以及信噪比SNR為30與20 dB的情況.TRR是一種非空間濾波的等效單源重建方法,表中未列出.從表2可見,信噪比對雙電流源估計的誤差有影響.表中用IIC方法同時估計兩個電流源的誤差較小,而用SONG和MVB的估計誤差相對較大.SNR=20 dB時,IIC的雙電流源估計誤差最大值E=2.44cm,比其他方法小.當SNR=30 dB時,IIC的雙電流源估計誤差相對SNR=20 dB時小.如圖4所示,IIC亦能夠仿真雙電流源沿參考路徑的電活動,偽源數量較少,其雙電流源重建的精度相對其他方法要高.

圖3 (a)黑色圓圈○和箭頭表示沿直線方向的參考路徑1;(b)?(e)是SNR=20 dB時,IIC,SONG,MVB和TRR方法的源重建結果.黑圈中的淺綠色到深綠色表示重建電流源位置隨時間的變化Fig.3.(a)The black circles ○ and arrow indicate the reference path 1 along the straight line;(b)?(e)The results of the current source reconstruction using IIC,SONG,MVB and TRR methods when SNR is 20 dB.The green points at different color levels denote the reconstructed source locations at different times.

4 心臟電活動的磁成像

本文利用P波間期采集的磁場信號,分析了健康人A和B的心臟電活動及其成像結果.這些數據是在磁屏蔽室內用Magnes 1300C 61通道生物磁強計系統沿Z方向測量的.信號采樣頻率1 kHz[2].假設P波間期心磁測量信號及其噪聲的平均功率為Ps和Pv,則信噪比 SNR=10lg[Ps/Pv].計算可得健康人A和B測量信號的P波間期信噪比 SNR∈[20,30]dB.

圖5中給出了健康人A和B的P波峰值(Ppeak)時刻心磁信號的源空間譜估計強度等高線圖.經歸一化處理后的等高線值用色標表示.圖中比較了用IIC,SONG和MVB三種空間濾波方法重建電流源的結果.為了顯示人體心臟及其房室的相對位置,圖中給出了一個健康人的心臟核磁共振影像(magnetic resonance imaging,MRI)冠狀位(coronal)和水平位(transverse)視圖,將其作為圖5中的XY與XZ面視圖,并分別與A和B的測量坐標系配準.圖5(c)和圖5(d)中A和B的心房內,用綠色三角形表示了冠狀位與水平位面的交點(0,?9,10)和(0,?6,10)cm.

表1 參考路徑1對應的單源估計誤差Table 1.The error of single source estimation correspond to the reference path 1.

圖5(c)和圖5(d)中給出了Ppeak時刻重建電流源的結果.IIC圖中,黃紅色標和緊密環繞的等高線表示較強電流源均集中在A和B的心房.健康人心房內的電活動較心室內強.這與Ppeak時刻健康人的心房正在除極,而心室尚未開始除極有關[33,34].然而,SONG的等高線擴散到了心臟以外,不像IIC那樣集中,說明SONG的空間譜強度估計的對比度相對較低,影響了重建電流源的精度.MVB的等高線圖顯示,心臟內外有很多無規律和強弱不同的偽源.重建電流源誤差相對較大是MVB電活動成像結果不理想的主要原因.

圖6(c)和圖6(d)是健康人A和B的P波前半段(Ponset—Ppeak,Ponset是P波起點)電活動成像結果.圖中用圓點表示重建源的位置,其中顏色由淺紅到深紅表示重建源位置隨時間的變化.圖6(c)和圖6(d)中IIC的結果顯示,P波前期健康人右心房除極時電活動從右心房的右上方向左移動.這與右心房的右上方靠近竇房結,右心房的左側靠近右心室有關[33].由于P波后半段(Ppeak—Pend,Pend是P波終點)的磁信號比較復雜,心房內電活動成像的特征不明顯,故文中沒有圖示.圖6(c)和圖6(d)中,SONG與IIC的電活動成像結果類似,A和B的心房外各有少量重建的偽源.圖6(c)中IIC和SONG顯示,雖然心臟P波前期磁信號呈單調上升,但是Ponset—Ppeak間期電活動的方向有變化.從起始時刻Ponset到結束時刻Ppeak,重建源移動方向如圖中綠色箭頭所示.MVB的成像結果比較模糊,所以A和B的心房內電活動特征不明顯.TRR是一種重建等效單源的方法,所以圖6(c)和圖6(d)中顯示,等效單電流源是從右上向左下移動的[32].實驗研究結果顯示,IIC圖6中P波間期健康人A和B心臟內分別有14個和26個極大電流源.相比其他方法超出心臟范圍的偽源較少.IIC和SONG的成像結果顯示了健康人右心房中電活動的方向.

5 討 論

用MVB,SONG和IIC研究源重建的結果表明,MVB采用的空間濾波器噪聲空間譜強度歸一化方法,SONG采用的抑制空間濾波器輸出噪聲功率增益方法,以及IIC采用的上述兩種方法,可以不同程度地抑制噪聲的影響,改進重建電流源的精度.表1、表2和圖3、圖4的仿真結果表明,IIC方法優于MVB和SONG.由(14)式和(15)式可知,IIC空間譜估計的源強度對比度比MVB源強度對比度大.

本文采用在濾波權矩陣中引入導聯場矩陣的方法,研究了有多個電流源的情況.圖3和圖4的仿真結果說明,IIC可以較好地消除重建源中的弱源與偽源,改進了2個源重建的精度.然而,當多個電流源的偶極矩方向變化時,多源重建的精度還需要深入研究.

表2 參考路徑2對應的雙電流源估計誤差Table 2.The error of two sources estimation correspond to the reference path 2.

圖4 (a)黑色圓圈○和箭頭表示沿直線方向的參考路徑2;(b)?(e)SNR=20 dB時,IIC,SONG,MVB和TRR方法的源重建結果.黑圈中的淺綠色到深綠色表示重建電流源位置隨時間的變化Fig.4.(a)The black circles ○ and arrow indicate the reference path 2 along the straight line;(b)?(e)The results of the current source reconstruction using IIC,SONG,MVB and TRR methods when SNR is 20 dB.The green points at different color levels denote the reconstructed source locations at different times.

IIC采用提高空間譜估計的源強度對比度和提取偶極矩極大源的方法去除重建的偽源.從圖3—圖6的源重建和電活動成像結果可見,這種方法比其他方法要好.仿真和實測成像的結果中,根據經驗,重建源強度的閾值分別等于當前時刻電流源空間譜的最大強度估計的0.01%和1%.圖6中用IIC時,健康人A的瞬時極大源數最多是2個,健康人B的最多是4個.

圖5中2個健康人的IIC的電活動成像結果顯示,參照MRI中人體心臟與其左右房室的位置,用黃色表示的Ppeak時刻較強電流源均集中在他們的右心房.圖6中紅色圓點表示Ponset—Ppeak間期強度極大源也集中在右心房.他們右心房的電活動均比左心房強.這與右心房靠近竇房結,即心臟電興奮的起搏點,以及健康人右心房的電興奮早于左心房相符合[33].圖6(c)中IIC和SONG的電活動成像結果顯示,雖然Ponset—Ppeak間期心磁信號單調上升,但是,這時電活動的方向有變化.可能這與參考文獻[33]指出的“心臟肌纖維的走向呈螺旋狀,以及在纖維的走向上,肌電阻率較低,肌電活動較強”有關.圖6中未能給出左心房電活動的成像結果,可能與心臟磁場只有Z軸方向的信號測量有關.

我們還用IIC方法分析了2個病人的數據[4].他們心臟左、右冠狀動脈嚴重狹窄,P波間期心房電活動與健康人明顯不同.可能是心臟冠脈嚴重狹窄引起心肌缺血,他們的心房電活動相對較弱,未顯示出心房中電活動的方向.因此,文中沒有給出病人的電活動成像結果,有關問題還需要深入研究.

6 結 論

本文針對心臟P波間期電流源重建及心臟電活動成像的問題,提出了一種可提高分布源空間譜估計強度對比度(IIC)的波束成形方法.該方法在空間濾波權矩陣中引入導聯場矩陣,可以改進分布源空間譜強度的對比度,有降低濾波器輸出噪聲功率增益的作用.再根據經驗設置源強度閾值,提取重建分布源中偶極矩強度極大的電流源,可以去除重建源中相對較弱的分布源和心臟范圍以外的偽源.可以提高P波間期電流源重建的精度.文中對4種方法的理論分析,源重建仿真和電活動成像結果的比較可見,IIC方法將有助于P波間期的心臟電流源重建和電活動的成像.對有局部心肌缺血的病人,心臟電活動成像的問題有待深入研究.

感謝中國科學院上海微系統與信息技術研究所的張懿教授,謝曉明教授和孔祥燕教授及張樹林博士對本研究的支持和技術交流.