基于信號子空間主特征向量的心磁源重構

李浩然，朱俊杰

河南理工大學電氣工程與自動化學院，河南焦作 454000

前言

心臟磁場研究的發展離不開超導量子干涉儀的成功研制，利用測量到的心臟磁場數據重構心臟的電流源，可視化心臟電流空間分布一直是心臟電生理活動逆向問題研究的重點和難點［1-5］。在心磁源重構的研究中有很多重構心臟電流源的方法，利用高精度的源重構算法重構心臟電流源，分析心臟電興奮時電流分布的情況，可視化心臟電活動的傳導過程，這一研究對心臟疾病的早期預防和診斷有重要的意義［6-7］。

波束形成器是一種常被用于心磁圖和腦磁圖源重構的成像方法。1997年，Van Veen等［8］在腦電研究的過程中利用線性約束最小方差波束形成（Minimum Variance Beamforming,MVB）方法對腦電信號源進行了定位；2001年，Sekihara 等［9］提出信號子空間特征投影的方法，并將該方法用于腦磁源定位的研究中；2006年，Kim 等［10］提出一種獨立成分分析的預處理方法，利用該方法對預激綜合征患者的傳導旁路進行定位；2011年，王偉遠等［11-13］針對QRS 波段心臟電活動變化快的特點，研究時間窗寬度、源移動速度以及信號噪聲對源重構精度的影響；2014年，趙晨等［14］提出一種用于等效源重構的磁場極值信號法，該方法可以用于在非均勻介質的情況下求解心磁逆問題；2018年，周大方等［15］提出一種抑制空間濾波器噪聲輸出功率增益的波束形成方法，該方法具有很好的噪聲抑制能力，提高了重建分布電流源偶極矩強度估計的源分辨率。

本研究提出一種改進的基于信號子空間主特征向量的波束形成（Improved Signal-subspace Principal Eigenvector based Beamforming,ISPEB）方法，此方法約束了空間濾波器輸出的噪聲功率和增益，補償了信號主特征向量對噪聲空間分布均勻性的影響。仿真實驗中，通過SPEB 和ISPEB 方法對心磁仿真數據進行源重構，估計其對應的分布電流源空間譜。最后采用MVB、SPEB 和ISPEB 方法對2 例沿z 軸正向測量的健康受試的單周期61通道心臟磁場數據進行源重構成像，并比較不同兩種方法的電活動成像結果。

1 波束形成方法

1.1 問題的提出

假設SQUID 系統在人體胸腔表面測量到的心臟磁場由n個呈網格狀分布的電流源偶極子源產生［16-20］，則檢測平面對應的電流源模型可以用如下線性方程表示:

其中，b(t)=[b1(t),b2(t),…,bc(t)]T是t時刻c通道陣列信號的列向量；q(t,rj)=q(t,rj)η(t,rj)(j= 1,2,…,n)表示n個分布電流源的偶極矩；rj=(xj,yj,zj)(j= 1,2,…,n)表示電流源空間的任意網格點位置；q(t,rj)=表示源偶極矩的強度；η(t,rj)=[ηX(t,rj),ηY(t,rj),ηZ(t,rj)]T是源偶極矩的方向；v(t)是t時刻測量的噪聲向量。下文中，q(t,rj)、η(t,rj)、q(t,rj)、v(t)和b(t)將簡寫為qj、ηj、qj、v和b。

本研究使用空間濾波的方法重建心臟電流源，將心臟磁場測量數據b(t) 作為空間濾波器的輸入，估計的分布源的偶極矩(t,rj)(j= 1,2,…,n)為輸出，空間濾波可用加權的線性運算表示：

其中，W(t,rj)T是空間濾波的權矩陣，式（2）可簡寫為：

MVB方法的基本原理是先用空間濾波技術重建心臟的分布電流偶極子源，根據可描述電流源偶極矩平均強度的分布電流源空間譜估計，對心臟電流源成像。電流源偶極矩與該位置上的電流源強度有關。

帶入式（3），得到空間濾波器的總輸出功率為：

其中，tr[ ·] 表示矩陣的跡。假設式（1）中測量信號的噪聲v(t)為高斯白噪聲，滿足E(v)= 0 與E(vvT)=σ02I，σ02是濾波器噪聲輸入功率，I是單位陣。將式（1）和E(v)= 0，E(vvT)=σ02I帶入式（4），可以表示為：

其中，前一項為空間濾波器輸出的所有電流源偶極矩qi對應功率和；后一項σ02tr[WjTWj]是濾波器輸出噪聲功率，其中tr[WjTWj]是輸出噪聲功率的增益。對式（5）開方，可得估計的電流源偶極矩平均強度空間譜：

式（6）中估計的電流源空間譜也可以估計為：

最小化空間濾波器輸出的總功率，歸一化空間濾波器輸出的噪聲空間譜強度，得到MVB 方法的空間濾波權矩陣：

SPEB 重建分布源方法的濾波權矩陣以MVB 濾波權矩陣為基礎。將MVB 權矩陣Wj,MVB的列向量投影到信號子空間主特征向量x1為基的一維子空間，可得到新的矩陣：

其中，x1是對稱陣E(bbT)最大特征值λ1(λ1∈(0,1))對應的主特征向量。滿足x1Tx1= 1，屬于信號子空間。將式（9）帶入式（7），令E(bbT)=I，對MVB的權矩陣Wj,MVB使用信號子空間主特征向量投影策略后，新的噪聲空間譜強度為：

用式（10）對應的參數標準化權矩陣Wj,1，歸一化每個分布源位置的空間濾波器輸出的噪聲空間譜強度，使其重新分布均勻，歸一化后的SPEB 空間濾波器權矩陣為：

1.2 ISPEB方法

本研究提出ISPEB方法，令空間濾波權矩陣為：

可以證明，其中的實對稱矩陣V=E(bbT)(x1x1T)-1是一個對空間濾波器輸出功率有影響的矩陣。由式（9）可知，因為最大特征值λ1相比其他特征值，受噪聲影響最小，故可以采用該功率較小的濾波權矩陣來增強分布源重建與空間譜估計的抗噪能力。因此，ISPEB方法中輸出噪聲功率增益可以表示為：

對式（13）使用施瓦茲不等式可以得知，ISPEB方法可以更好地降低空間濾波器的輸出噪聲增益，將E(bbT)=I帶入式（8）、（11）和（12）可知，ISPEB、SPEB 和MVB 的噪聲空間譜強度都相等，等于1。綜上所述，ISPEB方法可以約束噪聲空間譜對分布源空間譜的影響和濾波器輸出功率的增益，相對于MVB和SPEB 方法，ISPEB 方法可以提高估計空間譜的分辨率。

1.3 ISPEB方法的分析與比較

由于多電流源重建問題比較復雜，采用波束形成方法進行電流源的空間譜估計決定了電活動成像的分辨能力［21-24］。本研究比較了ISPEB和SPEB單電流源重建的源分辨率。

由式（2）可知，當單電流源偶極矩的方向已知后，空間濾波器的源偶極矩估計退化為源偶極矩的強度估計，其中，和wj是退化后空間濾波器的輸出和權向量。由式（7）可知，任意位置的單電流源空間譜強度估計為假設給定單電流源S 在單源rs位置上的強度估計為。定義任意位置rj(j= 1,2,…,n) 上的點擴散函數為歸一化后的估計空間譜可以表示為

由式（12）可知，單源重建時，采用的空間濾波權向量為：

由式（7）和（14）可得ISPEB方法的單源重構的空間譜強度估計：

其中，lj是電流源產生的磁場列向量。將式（15）進行化簡，得到ISPEB的空間譜估計強度表示如下：

歸一化后的點擴散函數為：

同理可得SPEB的點擴散函數：

比較式（17）和（18），有：

通常，由于所有測量通道上的理想磁場信號平均功率大于噪聲平均功率，所以可以令α=根據施瓦茲不等式(ljT,f)2≤可知，由此，根據式（18）和（19）可得：

當rj=rs時，點擴散函數可以取到最大值1，但此時仍然小于1，所以永遠小于1，即比小。當rj≠rs時，其他空間位置上，ISPEB 方法的點擴散函數仍然小于SPEB 方法的點擴散函數，點擴散函數可以反映空間其他位置的估計強度影響大小，取值小，說明單源對鄰域的影響擴散小。定義空間譜估計的單源分辨率為ψj=將式（20）帶入該式，得到：

因此，相比于SPEB 方法，ISPEB 方法有著更好的源分辨率。

2 仿真與實驗

2.1 仿真數據和相關源重構

相關電流源時比較難分辨，本研究利用仿真的磁場數據，比較了MVB、SPEB 和ISPEB 方法估計相關電流源的能力。

在仿真實驗中，采用雙電流偶極子仿真產生61通道的磁場數據（圖1），在3 組仿真數據中分別加入信噪比為10 dB 和20 dB 高斯白噪聲，源模型均在檢測平面下8 cm 處，檢測平面范圍為20 cm×20 cm。分別用3 種方法對雙電源偶極子進行源重構，并分析3種源重構的精度，采用重構出的電流源極矩強度較大的位置與給定電流源的位置誤差和重構電流源的分布面積作為評價源重構精度指標。

圖2是3 種方法對雙電流源偶極子重構得到的空間譜強度估計經過歸一化處理后的等高線表示，圖中的“×”表示給定的源位置?？梢钥闯?，在20 dB的信噪比下，SPEB 和ISEPB 方法明顯比MVB 方法有著較好的源重構能力；在10 dB的信噪比下，ISPEB方法有著最好的源重構能力。在低信噪比情況下，ISPEB方法的源重構能力最強。

圖3是SPEB 和ISPEB 方法對雙電流源偶極子重構的噪聲輸出功率增益抑制比。當信噪比為20 dB 時，SPEB 方法的噪聲輸出功率增益抑制比為627～657，ISPEB 方法的噪聲輸出功率增益抑制比為1625～1 780；當信噪比為10 dB 時，SPEB 方法的噪聲功率增益抑制比為635～657，ISPEB 方法的噪聲輸出功率抑制比為1 662～1 780。通過對比分析可知，當信號的信噪比從20 dB 下降到10 dB 時，ISPEB 方法的抑制能力更突出，因此，在多電流偶極子源重構中ISPEB 方法的噪聲抑制能力和源分辨能力更好。

2.2 心臟電活動成像

本文使用的心磁圖數據采集協議已得到當地倫理委員會批準，同時，受檢者給出了書面知情同意書。對兩名健康受試的61通道心臟磁場數據進行源重構分析。圖4是兩組健康者的心臟磁場數據圖，在整個心臟活動中，心臟磁場信號的R峰時刻心臟電活動較為劇烈，此時心室內電活動比心房的電活動強，易于觀察。分別用MVB、SPEB 和ISPEB 這3 種方法對兩名健康受試的心臟磁場數據R 峰時刻進行源重構。

圖5給出了兩名健康受試實測心臟磁場數據的源重構結果，圖中方框表示給定的源位置?？梢钥闯?，MVB方法重構出的心臟電流源分布效果最差，心臟磁場信號中的噪聲嚴重影響了源重構的結果。在信噪比較低的情況下，SPEB 和ISPEB 方法都有較好的成像結果。在R峰時刻，健康受試的心房除極已結束，心室正處于除極期，心室內黃色表示的電活動強度比心房內紅色表示的電活動強度高，心室與心房的源強度對比度增強了，心室電流源在心臟范圍內的分布及其電活動成像的清晰度也改善了。MVB方法成像結果比較模糊，上述電活動特征不明顯；而SPEB方法和ISPEB方法都有較好的清晰度。

3 結論

本研究在MVB 方法和SPEB 方法的基礎上，在SPEB 方法的權矩陣上添加了一個約束矩陣，使其在信噪比較低的情況下有著較好的成像結果，ISPEB方法利用一種低跡半正定矩陣構造了一個濾波器權矩陣，可以有效地降低噪聲功率增益，提高重建分布電流源偶極矩強度，即提高電流源空間譜估計的源分辨能力。實驗結果和分析比較表明，當信噪比較低時，ISPEB 方法在3 種方法中有著最好的成像結果，當信噪比較低時，可以采用ISPEB方法進行成像。