周大方 張樹林 蔣式勤?

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

1 引 言

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

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

2 波束成形方法

2.1 問(wèn)題的提出

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

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

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

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

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

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

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

2.2 SONG波束成形方法

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

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

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

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

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

2.3 SONG方法的分析與比較

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

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

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

式中,

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

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

其中

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

歸一化后的點(diǎn)擴(kuò)散函數(shù)為

同理可得,MVB單源重建的估計(jì)空間譜強(qiáng)度和歸一化后的點(diǎn)擴(kuò)散函數(shù)

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

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

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

3 仿真與實(shí)驗(yàn)

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

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

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

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

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

圖1 (a),(b)兩種SNR情況下,兩個(gè)相關(guān)電流源產(chǎn)生的仿真磁場(chǎng)數(shù)據(jù)(紅色虛線表示重建相關(guān)源的時(shí)刻);(c),(d)SONG和MVB在XY平面(z=11 cm)上估計(jì)空間譜強(qiáng)度的等高線圖(×號(hào)表示相關(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 心臟電活動(dòng)成像

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

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

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

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

圖2 (a),(b)兩個(gè)健康人的36通道單周期心臟磁場(chǎng)數(shù)據(jù)及對(duì)應(yīng)的Rpeak時(shí)刻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)矩陣對(duì)應(yīng)的噪聲空間譜強(qiáng)度也恒等于1,并可以約束空間濾波輸出噪聲功率增益,但是,仿真結(jié)果表明,其他濾波權(quán)矩陣會(huì)使兩個(gè)相關(guān)電流源定位到它們的中間位置.因此,必須利用測(cè)量信號(hào)的二階特征矩陣E(bbT)構(gòu)造空間濾波權(quán)矩陣[2,7].因?yàn)闇y(cè)量信號(hào)矩陣E(bbT)可以反映電流源的相關(guān)性.圖1中用SONG方法,兩個(gè)給定的相關(guān)源可以準(zhǔn)確定位.

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

圖3 (a),(b)兩個(gè)健康人的分布源空間譜估計(jì)強(qiáng)度的等高線圖;(c),(d)兩個(gè)健康人的空間濾波器輸出噪聲功率增益比β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ǔ)上,提出了一種用于心臟電活動(dòng)成像的可抑制空間濾波器輸出噪聲功率增益的波束成形方法.該方法利用一種低跡半正定矩陣構(gòu)造了一個(gè)濾波器權(quán)矩陣,可以降低空間濾波器輸出的噪聲功率增益,提高重建分布電流源偶極矩強(qiáng)度分辨率即分布電流源空間譜估計(jì)的源分辨能力,從而增強(qiáng)心臟電活動(dòng)磁成像的分辨率.仿真實(shí)驗(yàn)和分析比較的結(jié)果表明,SONG方法優(yōu)于MVB方法.當(dāng)心磁信號(hào)的SNR不低于10 dB時(shí),采用該方法可以提高心臟電活動(dòng)成像的效果,將有助于相關(guān)的醫(yī)學(xué)研究和應(yīng)用.

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