用于心臟電活動成像的空間濾波器輸出噪聲抑制方法?

周大方 張樹林 蔣式勤?

1)(同濟大學(xué)電子與信息工程學(xué)院,上海 201804)2)(中國科學(xué)院上海微系統(tǒng)與信息技術(shù)研究所,信息功能材料國家重點實驗室,上海 200050)(2018年2月6日收到;2018年4月8日收到修改稿)

1 引 言

波束成形方法是一種空間濾波技術(shù).它可以通過構(gòu)造一個空間濾波器,提取該位置上感興趣的源強度信息,同時抑制來自其他位置上源的影響[1,2].這種方法采用分布源模型估計分布電流源偶極矩的強度及其空間位置信息,即分布源的空間譜.相比采用單電流偶極子源模型,波束成形方法解決了估計多源的計算問題,可以細致地描述生物電活動,但是,面臨如何準(zhǔn)確重建多源的問題.近年來,最小方差波束成形(minimum variance beamforming,MVB)[1?7]被用于重建分布等效電流偶極子(equivalent current dipole,ECD)源的研究,例如,識別預(yù)激綜合征(wolf f-parkinson-white)病人心臟房室的旁路傳導(dǎo)[8,9],以及定位房顫(atrial fibrillation,AF)病人的AF病灶[10]和利用磁場圖(magnetocardiogram,MCG)仿真數(shù)據(jù)的心臟輪廓成像[7].在重建磁場分布電流源時,MVB通過自適應(yīng)的空間濾波技術(shù),最小化空間濾波器輸出總功率及歸一化噪聲空間譜強度,降低了測量噪聲對電流源空間譜估計的影響.可以發(fā)現(xiàn),如果在MVB方法的基礎(chǔ)上,再針對性地約束空間濾波器輸出的噪聲功率增益,可以進一步提高分布電流源空間譜估計的源分辨能力,從而增強心臟電活動磁成像的分辨率.

本文提出了一種可抑制空間濾波器輸出噪聲功率增益(suppressing spatial f i lter output noisepower gain,SONG)的波束成形方法.用一個對空間濾波器輸出功率有影響的低跡的半正定矩陣構(gòu)造一種新的濾波器權(quán)矩陣,可以約束空間濾波器的輸出噪聲功率的增益.該低跡的半正定矩陣滿足特征值不大于1,且矩陣的跡小于其階數(shù).這樣,可以提高分布電流源空間譜估計的分辨率.文中通過理論分析和仿真實驗比較了SONG和MVB方法.給出了用兩個健康人的36通道MCG數(shù)據(jù)得到的心臟電活動成像.結(jié)果表明,SONG方法分辨電流源的能力較強,能夠觀察到Rpeak時刻健康人的心室內(nèi)有較強電活動等明顯的電生理特征.

2 波束成形方法

2.1 問題的提出

假設(shè)第j個單位電流偶極子源的位置和方向分別是rj=(xj,yj,zj)(j=1,2,···,n)和[1,0,0]T,[0,1,0]T或[0,0,1]T.它產(chǎn)生的理想磁場列向量為lX,j,lY,j或 lZ,j,相應(yīng)的導(dǎo)聯(lián)場矩陣表示rj處電流源的測量靈敏度.已知心臟磁場的測量數(shù)據(jù)b(t)[11,12],求電流源偶極矩q(t,rj)的逆問題[13,14],可用線性方程表示[15,16]:

其中,b(t)=[b1(t),b2(t),···,bc(t)]T表示t時刻c個測量通道的磁場向量;v(t)是t時刻的測量噪聲向量;等效電流源的偶極矩q(t,rj)=q(t,rj)η(t,rj),j=1,2,···,n, 其中n是電流源的數(shù)目.源偶極矩的強度是標(biāo)量

表示單位向量.rj=(xj,yj,zj)是第j個電流源位置的坐標(biāo).下文中,q(t,rj),η(t,rj),q(t,rj),v(t)和b(t)將簡寫為qj,ηj,qj,v和b.

空間濾波是一種常用的分布源重建方法.將心臟磁場測量數(shù)據(jù)b(t)作為空間濾波器的輸入,估計的分布源的偶極矩作為輸出.空間濾波可用加權(quán)的線性運算表示:

其中,W(t,rj)是空間濾波的權(quán)矩陣.(2)式可簡寫為

MVB方法的基本原理是先用空間濾波技術(shù)重建心臟的分布電流偶極子源.然后,根據(jù)可描述電流源偶極矩平均強度的分布電流源空間譜估計,對心臟電流源成像.電流源偶極矩是源電流密度在鄰域的體積分,它與該位置上的源電流密度大小,也就是電流源強度有關(guān)[13,14].

式中,前一項為空間濾波器輸出的所有電流源偶極矩qi的功率和;后一項為空間濾波器輸出的噪聲功率,其中為噪聲功率的增益.對開方,可得估計的電流源偶極矩平均強度空間譜:

2.2 SONG波束成形方法

本文提出了一種SONG波束成形方法.令空間濾波權(quán)矩陣

可以證明,其中的c階實對稱矩陣V=E(bbT)是一個對空間濾波器輸出功率有影響的半正定矩陣,滿足“矩陣的跡低于其階數(shù),且其特征值不大于1”,(以下簡稱V是低跡半正定陣).因此,SONG方法中輸出噪聲功率增益可表示為

令任意矩陣Γ =Wj,MVB,也可以證明,存在不等式

當(dāng)V=E(bbT)時,由(10)和(11)式可知

也就是說,SONG方法可以或更好地降低空間濾波器的輸出噪聲功率增益.將E(bbT)=I代入(8)和(9)式可知,SONG的噪聲空間譜強度與MVB相同,均等于1.綜上,SONG方法不僅可以約束噪聲空間譜對分布源空間譜估計的影響,還可約束空間濾波器輸出噪聲功率的增益.相比MVB,SONG方法可以提高估計空間譜的分辨率.

2.3 SONG方法的分析與比較

采用波束成形方法,電流源的空間譜估計決定了電活動成像的分辨能力.由于多電流源重建問題比較復(fù)雜,本文比較了SONG和MVB單電流源重建的源分辨率.

由(2)式可知,當(dāng)單電流源偶極矩方向已知時,空間濾波器的源偶極矩估計可以退化為源偶極矩的強度估計是退化后空間濾波器的權(quán)向量.由(7)式可知,單源在任意位置上產(chǎn)生的估計空間譜強度為假設(shè)單源S位置rs上的估計空間譜強度為b ?s.為了分析估計空間譜的單源分辨率,定義任意位置rj(j=1,2,···,n)上的點擴散函數(shù)(point spread function)為(rj)[4],可用歸一化后的估計空間譜強度簡單表示為

由(9)式可知,單源重建時,SONG方法的權(quán)向量為

式中,

其中,lj是電流源產(chǎn)生的磁場列向量.

由(7)和(13)式可得SONG方法估計的空間譜

其中

將(15)和(17)式代入(14)式,可得SONG方法在任意位置rj(j=1,2,···,n)的估計空間譜強度

歸一化后的點擴散函數(shù)為

同理可得,MVB單源重建的估計空間譜強度和歸一化后的點擴散函數(shù)

比較(19)和(20)式,有

通常,所有測量通道上的理想磁場信號平均功率大于噪聲平均功率[11,12].所以,可令α=根據(jù)施瓦茲不等式性質(zhì),可知由此,從(20)和(21)式可以得到

可見,當(dāng)rj=rs時,單源位置上點擴散函數(shù)和相等,為最大值1. 當(dāng)rj?=rs時,其他空間位置上,SONG方法的點擴散函數(shù)比MVB的小.點擴散函數(shù)?j可以反映單源對空間其他位置的估計空間譜強度影響大小(rj?=rs)的取值小,說明單源對鄰域的影響擴散小.因此,相比MVB,SONG方法對單源的空間分辨率較高.

3 仿真與實驗

3.1 仿真數(shù)據(jù)和結(jié)果

相關(guān)電流源(correlated sources)是比較難分辨的,所以,文中利用仿真的磁場數(shù)據(jù),比較了SONG與MVB方法估計相關(guān)電流源的能力.

假定軀干G0={(x,y,z)|x∈[?12.5,12.5],y∈[?12.5,12.5],z∈[3,20]}(cm).隨機給定兩個相關(guān)源的位置坐標(biāo)為(6.5,?2.5,11)和(?2.5,6.5,11)(cm),它們的相關(guān)系數(shù)為0.9824.G0被劃分為間距1 cm的10625個體素(voxel).并用ECD源模型產(chǎn)生兩組仿真的36通道Z軸測量磁場數(shù)據(jù),采樣頻率為1 kHz[7,12],如圖1(a)和圖1(b)所示.假設(shè)其中測量噪聲分別為根均方(root-mean-square,rms)信噪比(signal-to-noise ratio,SNR)20 dBrms和10 dBrms的高斯白噪聲[7,12].產(chǎn)生仿真的磁場數(shù)據(jù)和進行電流源空間譜估計,都要用到等效電流源的導(dǎo)聯(lián)場矩陣.文中采用水平分層導(dǎo)體作為軀干模型,并通過解磁場的正問題求導(dǎo)聯(lián)場矩陣[7,13,14].

圖1(c)和圖1(d)是SONG和MVB在XY平面(z=11 cm)上歸一化后的空間譜估計強度的等高線圖,等高線的步長為1%.結(jié)果表明,SNR對相關(guān)電流源的空間譜估計有影響.與MVB相比,用SONG方法估計的相關(guān)源強度比鄰域的估計源強度有明顯的增強.MVB方法中,表示濾波器的輸出噪聲功率增益,本文在此基礎(chǔ)上給出了SONG與MVB方法的輸出噪聲功率增益比

如圖1(e)和圖1(f)所示.其中,橫坐標(biāo)是分布源位置的索引,縱坐標(biāo)是增益比β.將(12)式代入β可知β>0.當(dāng)β>0時,說明SONG方法抑制噪聲的效果比MVB好.圖1表明,測量信噪比SNR分別為20和10 dB時,用SONG重建相關(guān)源時,濾波器輸出噪聲增益比β分別為584—622和592—623.在SNR=10 dB時抑制噪聲增益的效果較SNR=20 dB時更明顯.因此,圖1(c)和圖1(d)表明,SONG方法有較強的相關(guān)源分辨能力,并且在SNR較低時,其源分辨能力更加明顯.

圖1 (a),(b)兩種SNR情況下,兩個相關(guān)電流源產(chǎn)生的仿真磁場數(shù)據(jù)(紅色虛線表示重建相關(guān)源的時刻);(c),(d)SONG和MVB在XY平面(z=11 cm)上估計空間譜強度的等高線圖(×號表示相關(guān)源的給定位置);(e),(f)用SONG方法的空間濾波器輸出噪聲功率增益比βFig.1.(a),(b)Simulated magnetic data generated by two correlated current sources with noise,respectively.The red dashed line denotes the time of source reconstruction.(c),(d)Estimated spatial spectrum intensity contour on XY plane(z=11 cm)using SONG and MVB.The black cross signs×denote the given locations of two correlated sources.(e),(f)Ratio β of noise-power gain of spatial f i lter output using SONG method.

3.2 心臟電活動成像

分布電流偶極子源的偶極矩強度可以反映分布電流的強度,因此,可利用分布源空間譜估計的強度對心臟電活動成像.文中還用兩個健康人的心臟磁場數(shù)據(jù)比較了SONG和MVB成像的效果.在求導(dǎo)聯(lián)場矩陣時,采用水平分層導(dǎo)體作為軀干模型[7,13,14].用這種模型時,導(dǎo)體邊界上的單位法向量均平行于心臟測量磁場的單位向量,所以,求解導(dǎo)聯(lián)場矩陣時,軀干體電導(dǎo)的影響可以忽略[7,13,14].

圖2是沿Z軸測量的兩個健康人的36通道單周期心臟磁場曲線.測量平面為25 cm×25 cm.36通道的SNR約為10—20 dB,帶寬0.01—100 Hz,采樣頻率1 kHz[12].

圖2中心臟磁場信號的Rpeak時刻對應(yīng)心室的除極期.這時健康人心室電活動較強,心房的相對較弱[17],比較容易識別[11,12].因此,文中比較了SONG和MVB兩種方法的成像結(jié)果.將測量平面與健康人核磁共振影像(magnetic resonance imaging,MRI)的心臟冠狀位(coronal view)、水平位(transverse view)和矢狀位(sagittal view)圖的坐標(biāo)配準(zhǔn),然后,用MRI中心臟的位置作為房室位置的參考.圖3(a)和圖3(b)給出了歸一化后的分布源空間譜估計強度的等高線圖.其中,冠狀位視角的XY面(z=10 cm)、水平位視角的XZ面(y=2.5 cm)和矢狀位視角的Y Z面(x=0.5 cm)的交點(0.5,2.5,10)cm位于心室內(nèi),用藍色的小方框表示,譜強度等高線的步長為1%.圖3(a)和圖3(b)表明,SONG方法的成像結(jié)果能夠反映健康人Rpeak時刻電活動的特征.因為Rpeak時刻,健康人的心房除極已結(jié)束,心室正處于除極期.根據(jù)圖中色標(biāo),可以觀察到心室內(nèi)黃色表示的電活動強度明顯比心房內(nèi)紅色表示的電活動強度高,以及心室鄰域的電活動相對較強[7,11,12,17,18].由于SONG增加了心臟內(nèi)外分布源空間譜估計強度的差異,心臟分布電流源的分辨率提高了.用MVB的成像結(jié)果相對模糊,特征不明顯.

圖3(c)和圖3(d)的結(jié)果表明,Rpeak時刻用SONG方法,兩個健康人的濾波器輸出噪聲增益比β分別為206—228和215—232.也就是說,當(dāng)實測心磁數(shù)據(jù)SNR為10—20 dB時,相比MVB,SONG能夠更好地抑制空間濾波器輸出噪聲.如圖3(a)和圖3(b)所示,SONG方法的源分辨能力以及成像效果相對較好.

圖2 (a),(b)兩個健康人的36通道單周期心臟磁場數(shù)據(jù)及對應(yīng)的Rpeak時刻Fig.2.(a),(b)The 36-channel MCG data of single-cycle from two healthy subjects,as well as the time-points of Rpeakfor MCG imaging.

4 討 論

由(12)式可見,SONG波束成形方法可以約束空間濾波器輸出噪聲的功率增益.用構(gòu)造濾波權(quán)矩陣Wj,SONG類似的方法,還可以構(gòu)造其他的濾波權(quán)矩陣.雖然理論上其他空間濾波權(quán)矩陣對應(yīng)的噪聲空間譜強度也恒等于1,并可以約束空間濾波輸出噪聲功率增益,但是,仿真結(jié)果表明,其他濾波權(quán)矩陣會使兩個相關(guān)電流源定位到它們的中間位置.因此,必須利用測量信號的二階特征矩陣E(bbT)構(gòu)造空間濾波權(quán)矩陣[2,7].因為測量信號矩陣E(bbT)可以反映電流源的相關(guān)性.圖1中用SONG方法,兩個給定的相關(guān)源可以準(zhǔn)確定位.

我們還利用MCG仿真數(shù)據(jù),研究了重建分布源的最小范數(shù)空間濾波(minimum norm spatial f i ltering,MN)方法[19,20],這種非自適應(yīng)的空間濾波方法,采用了Tikhonov正則化技術(shù).仿真結(jié)果表明,有可利用的心臟三維輪廓時,MN的分布電流源成像結(jié)果較好[19].參考文獻[19]的圖3和圖4給出了重建的16個等效電流源,與本文圖3中SONG方法的成像結(jié)果類似.

圖3 (a),(b)兩個健康人的分布源空間譜估計強度的等高線圖;(c),(d)兩個健康人的空間濾波器輸出噪聲功率增益比βFig.3.(a),(b)Contour of the estimated intensity of distributed source spatial spectrum of two healthy subjects;(c),(d)spatial f i lter output noise-power gain ratio β of two healthy subjects.

5 結(jié) 論

本文在研究MVB方法的基礎(chǔ)上,提出了一種用于心臟電活動成像的可抑制空間濾波器輸出噪聲功率增益的波束成形方法.該方法利用一種低跡半正定矩陣構(gòu)造了一個濾波器權(quán)矩陣,可以降低空間濾波器輸出的噪聲功率增益,提高重建分布電流源偶極矩強度分辨率即分布電流源空間譜估計的源分辨能力,從而增強心臟電活動磁成像的分辨率.仿真實驗和分析比較的結(jié)果表明,SONG方法優(yōu)于MVB方法.當(dāng)心磁信號的SNR不低于10 dB時,采用該方法可以提高心臟電活動成像的效果,將有助于相關(guān)的醫(yī)學(xué)研究和應(yīng)用.

感謝中國科學(xué)院上海微系統(tǒng)與信息技術(shù)研究所的張懿教授、謝曉明教授和孔祥燕教授及其團隊為本研究提供可用的心臟磁場數(shù)據(jù)與核磁共振影像數(shù)據(jù),以及有關(guān)的技術(shù)交流.