基于相位展开的磁敏感加权图像中相位伪影处理方法的研究

聂生东田晓迪钱黎俊陈增爱许建荣

1(上海理工大学医学影像工程研究所，上海 200093)2(上海交通大学附属仁济医院放射科，上海 200127)

基于相位展开的磁敏感加权图像中相位伪影处理方法的研究

聂生东1*田晓迪1钱黎俊2陈增爱2许建荣2

1(上海理工大学医学影像工程研究所，上海 200093)2(上海交通大学附属仁济医院放射科，上海 200127)

针对Hanning滤波器难以有效解决磁敏感加权成像过程中产生的相位伪影问题，提出一种基于相位展开技术的相位伪影处理方法。采用两个步骤有效消除磁敏感加权成像过程中产生的相位伪影。首先，采用传统的Hanning滤波器对相位伪影进行初步处理;然后，采用相位展开技术，对由Hanning滤波器处理后的残余伪影进行进一步处理，从而达到完全消除相位伪影的目的。对扫描层厚为0.6 mm的100层脑部SWI原始图像进行相位伪影的处理，结果表明，该方法能有效消除磁敏感加权图像中的相位伪影。采用该方法处理SWI图像中的相位伪影，可为后续采用最小密度投影方法显示清晰的血管影像奠定基础。

磁敏感加权成像;相位伪影;Hanning滤波器;相位展开

Abstract:The phase artifact，produced during susceptibility weighted imaging(SWI)，can not be eliminated efficiently by a Hanning filter.To solve this kind of problem，a phase artifact processing method was presented in this paper based on phase unwrapping technology.The method includes two steps.Firstly，a traditional Hanning filter was adopted to process the phase artifact.Then，the remnant phase artifact was processed using phase unwrapping technique.The proposed method was evaluated using 100 thin-section brain SWI slices with the interval of 0.6 mm.The results showed that the phase artifact could be completely eliminated from SWI images.The application of our method could make a basis for displaying clear vein blood vessel image by using minimum intensity projection technique.

Key words:susceptibility weighted imaging;phase artifact;Hanning filter;phase unwrapping

引言

磁敏感加权成像(susceptibilityweighted imaging，SWI)是近年来发展起来的一种全新的磁共振成像方法［1-2］。不同于以往的质子密度、T1或T2加权成像，这种方法可利用不同组织间磁敏感性的差异产生图像对比，进而可对各组织显影的新技术。与传统磁共振成像方式相比，SWI在静脉显影方面具有独特的优势，可应用于脑肿瘤、脑出血或其他有静脉参与的病灶研究中，从而有效地改善对这些疾病的诊断。静脉血的主要成分为顺磁性的去氧血红蛋白，动脉血则是抗磁性的氧合血红蛋白［3-5］。这种磁敏感性差异将最终导致两种血管信号的相位信息强度不同，在与幅度信息相加权后，使静脉能独立于动脉清晰成像，这是SWI成像的理论基础。在临床医学诊断中，由外伤所导致的微出血以及很多肿瘤病变都与静脉有关，目前的成像手段对静脉血管的成像效果都不甚理想，而利用SWI可获得清晰的静脉血管影像，因此SWI具有广泛的临床应用价值。

图2 相位图像中的相位伪影及处理效果。(a)未经处理的相位伪影;(b)经Hanning滤波器处理后的残留伪影;(c)伪影经相位展开处理后的效果Fig.2 Phase artifact and the processed effect.(a)phase artifact without process;(b)the remnant phase artifact after processed by hanning filter;(c)processed result after phase unwrapping

由于SWI是一种新的磁共振成像技术，现有的磁共振成像仪得到的仅仅是SWI原始图像，其所配备的图像处理工作站没有对SWI原始图像进行处理的相关软件，目前国内外也没有成熟的软件实现对SWI原始图像的后处理［5］。所谓的 SWI原始图像是指采用梯度回波序列扫描获得的幅值图像(magnitude image)和相位图像(phase image)，如图1中的(a)和(b)所示。常规磁共振成像中所使用的都是幅值图像，它描述弛豫过程中质子发出的信号强度，相位图像则描述质子在该过程中行经的角度。幅值图像中包含了绝大部分的组织对比信息，而相位图像则从磁敏感性的角度反映组织对比，特别是磁化率差异较大的组织。这两种图像是在扫描过程中同时获得的，所以，它们总是成对出现，每一对图像所对应的解剖位置都完全一致。但是，从图1中的(a)和(b)很难获得血管的诊断信息。要利用原始SWI图像获得连续清晰的血管影像，需要对原始图像进行一系列复杂的后处理过程。

在对SWI原始图像的处理中，需要解决的关键问题就是如何有效处理SWI相位图像中的相位伪影(或 field inhomogeneity artifact，磁场不均匀性伪影)。相位伪影是由于磁场的不均匀性造成的，通常出现在具有不同磁化率物质的交界面周围，如副鼻窦、蝶鞍和颅底等位置［5］。在成像过程中由于磁场的不均匀性，导致了在SWI的相位图像中出现相位模糊(或相位缠绕)问题，在相位图上，表现为环绕在交界面周围的扰动波纹，如图2(a)所示(图中圈出的部分)。近年来，对 SWI原始图像进行处理的文献报道非常少，且基本上都在国外文献上发表。目前，已报道的对SWI原始图像处理的文献中通常采用 4 个步骤［1，5］:

1)幅值和相位图像的重组;

2)消除相位图像中磁场不均一性伪影;

3)制作相位蒙片并与幅值图像加权获得磁敏感加权图像;

4)通过最小密度投影显示连续层面的静脉血管结构。

其中，第2步是一个关键步骤，其对伪影的处理结果直接影响到后续的采用最小密度投影方法显示血管影像的效果。在目前报道的文献中，都是采用Hanning滤波器来处理SWI相位图像中的伪影问题。利用Hanning滤波器处理SWI相位图像中的相位伪影，存在的主要问题是:经Hanning滤波器处理后，相位伪影不能完全被消除，仍然残留部分伪影，如图2(b)所示。如果不能完全消除伪影，残留伪影将会保留在后续生成的蒙片上(如图3(a)所示)，将会对最终血管显示结果造成不利影响。为此，本研究提出基于相位展开技术的SWI图像中相位伪影的处理方法，首次将广泛应用在合成孔径雷达干涉测量、合成孔径声纳、光学干涉、微波干涉等领域的相位展开技术应用于处理 SWI图像中的伪影［6，7］，成功解决了 SWI图像中的相位伪影问题，如图2(c)所示。由相位展开消除伪影后生成的蒙片，残留的伪影明显被消除，如图3(b)所示。

图3 后续生成的相位蒙片。(a)由图2(b)生成的相位蒙片;(b)由图2(c)生成的相位蒙片Fig.3 The different phase mask produced by Fig.2(b)and Fig.2(c).(a)the phase mask produced by Fig.2(b);(b)the phase mask produced by Fig.2(c)

1 相位展开技术及其应用

1.1 关于相位缠绕与相位展开问题

关于相位缠绕(wrapping)问题，可以用图4进行简单说明。一个相位本应在0～4π的信号因为某种原因，落在了-π到+π之间，这就是所谓的相位缠绕问题。将此缠绕信号重新恢复到原相位范围，称为相位展开(phase unwrapping)，有些文献称为相位解缠。

图4 真实相位与缠绕相位Fig.4 The sketch map of the true phase and the intertwist phase

在合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)中存在相位缠绕问题［6-13］。合成孔径雷达干涉测量技术，是以合成孔径雷达复数据提取的相位信息为信息源，获取地表的三维信息和变化信息［6］。InSAR通过两副天线，同时观测(单轨模式)或两次近平行观测(重复轨道模式)获取地面同一景观的两幅包含幅值与相位信息的二维复数据图像［6-7，13］。因此，与 SWI 成像一样，InSAR一次测量也得到同一景观的一对图像，即幅值图像和相位图像。由于相位展开技术能有效解决合成孔径雷达干涉测量中的相位缠绕问题，因此应用相位展开技术解决 SWI图像中的相位缠绕(相位伪影)问题也是一个自然的选择。

1.2 纳入标准 （1）年龄≥60岁。（2）本次住院行冠脉造影，且造影结果显示冠脉血管1支及（或）以上血管病变≥50%。（3）患者冠心病诊断符合我国2010年ST段抬高型心肌梗死诊疗指南［9］中或我国2012版非ST段抬高型急性冠脉综合征［10］的标准。（4）患者糖尿病诊断均符合我国2型糖尿病防治指南（2013年版）［11］中的相关标准。（5）临床资料及随访资料完整。（5）对本研究知情同意，且愿意按时接受随访。

与非相干信号处理只关心信号的幅度与强度不同，相干信号处理还同时涉及信号的时间和空间的相位信息［8］。SWI成像和干涉孔径雷达成像是相干信号处理的典型代表，然而因为系统原因，所得到的与感兴趣的物理量相联系的相位均被缠绕在(-π，π］之间，存在着相位模糊(相位缠绕)的现象。恢复与磁场一致的真实相位，就是本研究所要讨论的相位展开问题。

1.2 相位展开技术简述

从数学上讲［6，13］，给定一个二维相位矩阵 Ψi，j，为了对该矩阵进行相位展开，需要对每个像素点(i，j)加上或减去2π的整数倍，从而得到一个连续的函数:

在解决相位缠绕的问题中，目前有很多方法［6，13］，但其中主要采用的方法有路径跟踪法和最小范数法。路径跟踪法(path-following method)是建立在残差定理的基础上，通过建立枝切线，阻止积分路径穿过相位不连续区域，从而实现与路径无关的相位积分，得到没有误差传递的相位展开结果。其代表算法包括Goldstein的枝切法(branch cut)［9］、Xu 和 Cumming的区域生长法(region growing)等。最小范数法则基于最小均方意义，其思想是最小化测量相位梯度与展开相位梯度的均方差，主要工作集中在 Romero、Franceschetti及Lanari的格林函数法，Ghigilia、Romero的先验条件共轭梯度法(PCG)［10］，Pritt于 1996 年提出的加权多重网格算法［11］。此外，通过求相位数据的二次偏导数并设置阈值，从而生成二值质量图，这就是1991年Bone提出的利用质量评价来指导相位展开的质量图引导法［14］。

为了从众多相位展开方法中选择能够有效处理SWI图像中相位伪影的方法，分别采用Goldstein枝切法、质量图引导法和最小二乘法等不同算法，对SWI相位数据进行对比处理实验，结果发现，以Goldstein的枝切法为代表的路径积分法对SWI相位数据的处理效率最高［15］。因此，下面主要介绍Goldstein的枝切法在SWI相位伪影处理中的应用。

1.3 基于路径跟踪的相位展开算法

Goldstein等提出的枝切法是基于路径跟踪的相位展开方法的典型代表［9］。其基本思想是:首先，通过识别正负残差点，并连接邻近的正负残差点，实现残差点的电荷平衡;然后，生成最优枝切线来确定积分路径，防止误差沿积分路径传递;最后，使相位展开的积分路径不穿过枝切线，从而得到真实的相位信息。算法基本步骤包括识别残差点、生成枝切线和绕过枝切线积分，以下对算法进行详细描述。

首先，遍历图像寻找残差点，当扫描到残差点时，将一个3像素×3像素的窗口放置在该残差点上，在窗口内寻找其他的残差点，当寻找到新的残差点时将它们连接起来;接着，判断极性，如果极性相反，则表明通过连接枝切线实现了残差点的极性平衡，如果极性相同且此点未与其他残差点相连，则将其加入到极性表中记录其极性，并将其电荷值累计;然后，将窗口移到新的残差点上继续搜索，直到达到极性平衡，否则将窗口放大到5像素×5像素、7像素×7像素或更大范围来搜索。当搜索到边界时，将残差点与边界相连。生成枝切线后，通过积分相邻像素上的相位差分进行相位展开。其中的关键是:无论在窗口中找到的残差点是否与其他残差点相连，都将它与窗口中心的残差点相连。另外，当搜索窗口到达图像边界时，将残差点与边界相连，阻止积分路径穿过。图5是路径跟踪法框图，图6描述了枝切线设置的基本过程。

图5 路径跟踪法框图Fig.5 The flow chart of the path-following method

2 实验结果及讨论

对100层SWI相位图片中的第1层运用Goldstein方法进行搜索，产生残差点和枝切线结果。图7是第一次对图2(a)搜索时产生的残差点结果图像，其中被圈起及划线的区域是被消去的对偶极子较多的区域;图8是消除对偶极子后所剩的残差点;图9为消除对偶极子后产生的枝切线。

图2(c)是对图2(a)进行识别残差点、消除对偶极子、设置枝切线后相位展开的结果图像。图10(a)是利用Hanning滤波器对相位伪影处理后，再利用最小密度投影获得的图像。图10(b)是利用相位展开技术，对经Hanning滤波器处理后的残留伪影进行处理，再利用最小密度投影获得的图像。从图10可以明显看出，经相位展开技术处理后，可以获得更清晰的血管显示影像。

图11～图13是利用枝切法对40层图像处理的结果。

对比图11和图12可以看出，利用本方法可以有效地消除SWI图像中的相位伪影，经最小密度投影处理后，可得到清晰的血管影像(见图13)。

另外，本研究从算法的处理速度和处理效果两个方面，对常用的两种算法——枝切法和质量图引导法的处理结果进行了对比实验。从视觉效果上看，两种算法对SWI图像中相位伪影的处理效果基本相同。在处理速度方面，利用枝切法处理每幅SWI图像的平均时间是476 ms，而利用质量图引导法处理每幅SWI图像的平均处理时间为736 ms。因此，对SWI图像中相位伪影的处理速度，枝切法比质量图引导法要高得多。

图6 枝切线设置的基本过程。(a)正负残差点;(b)连接正负残差点;(c)正负残差点;(d)3像素×3像素搜索窗;(e)连接残差点;(f)搜索窗扩大至3像素×3像素;(g)3个残差点;(h)连接3个残差点;(i)将残差点接地Fig.6 The basic process of set up branch cut line.(a)the negative＆positive residues;(b)to connect the negative＆positive residues;(c)the negative＆positive residues;(d)the 3 pixels×3 pixels searching window;(e)to connect the negative＆positive residues;(f)extending the searching window to 5 pixels×5 pixels;(g)three residues;(h)to connect the three residues;(i)grounding the residues.

图7 第1次生成的残差点Fig.7 The residues produced at the first time

图8 消去对偶极子后的残差点Fig.8 The residues after the dipole was deleted

最小二乘相位展开算法是求解二维相位展开问题最稳健的方法之一，也是一种在雷达干涉孔径成像中常用的相位展开技术［16］，笔者也试图在 SWI图像的相位伪影处理中应用这一技术。但通过实验发现，使用最小二乘法进行相位展开时，算法平均运行时间在60 s以上，运行耗时太高。又因所处理的SWI图像为180层连续图像，如果只进行相位展开处理就需要60×180 s的时间，而 SWI后处理还包含其他步骤。因此，在处理效果基本相同的情况下，运行时耗的多少是决定SWI的后处理算法选取的一个重要因素。由于最小二乘法相对较低的处理效率，研究中没有选取最小二乘法对SWI图像进行相位展开处理。

图9 消除对偶极子后的枝切线Fig.9 The branch cut line after the dipole was deleted

图10 经最小密度投影处理后的图像显示结果。(a)有残留伪影的血管显示图像;(b)消除残留伪影后的血管显示图像Fig.10 The blood vessel image display results after processed by minimum intensity projection.(a)blood vessel display image with remnant phase artifact;(b)blood vessel display image after remnant phase artifact was eliminated

图11 处理前的相位图像Fig.11 The phase image before processed

图12 消除伪影后的相位图像Fig.12 The phase image after artifact be eliminated

图13 经最小密度投影处理后的血管图像Fig.13 The blood vessel image display results after processed by minimum intensity projection

3 结论

SWI是一种新型的磁共振成像方式，其主要特点就是利用不同物质间的磁敏感性差异进行静脉的显影成像。目前，SWI在临床上应用的瓶颈问题是图像后处理技术，特别是如何消除磁场不均一性伪影问题。近年来，有关处理SWI图像的文献都是采用Hanning滤波器，滤除部分磁场不均一性伪影。由于经Hanning滤波器处理后仍残留部分伪影，严重影响了最终的血管显示结果。因此，本研究首次在SWI图像后处理中引入了在雷达干涉孔径成像中常用的相位展开技术，专门来解决SWI相位图像中的磁场不均一性伪影问题。实验结果表明，相位展开技术可有效处理磁场不均一性伪影，为最终显示清晰的血管影像奠定了基础。

［1］Thomas B，Somasundaram S，Thamburaj K.Clinical applications of susceptibility weighted MR imaging of the brain-a pictorial review［J］.Neuroradiology，2008，50:105-116.

［2］Shrinivas B.SWI，a new MRI sequence—How useful it is?［J］.Ind J Radil Imag，2006，16(1):13-14.

［3］Ogawa S，Lee TM，Nayak AS，et al.，Oxygenation sensitive contrast in magnetic resonance image of rodent brain at high magnetic fields［J］.Magn Reson Med，1990，14(1):68-78.

［4］Ogawa S，Lee TM，Kay AR，et al.Brain magnetic resonance imaging with contrast dependent on blood oxygenation［J］.Proc Natl Acad Sci，1990，87(24):9868-9872.

［5］钱黎俊.磁敏感加权成像后处理算法实现与临床应用［D］.上海:上海交通大学，2007.

［6］孟智强，夏雨人.相位质量图在二维相位解缠中的应用探讨［J］.计算机应用与软件，2006，23(2):63-65.

［7］岑小林.InSAR相位解缠算法研究［D］.长沙:湖南大学，2008.

［8］李笑郁，毛士艺.干涉SAR与MRI中相位展开算法研究［J］.中国体视学与图像分析，2001，4(6):193-197.

［9］Goldstein RM，ZebkerHA，WernerCL.Satellite radar interferometr:two-dimensional phase unwrapping［J］.Radio Sci，1988，23(4):713-720.

［10］Ghiglia DC，Romero LA.Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods［J］.J Opt Soc Amer A，1994，11(1):107-117.

［11］Pritt MD.Phase unwrapping by means of multigrid techniques for Interferometric SAR［J］.IEEE Trans on GRS，1996，34(3):728-738.

［12］孙龙.偏微分方程在合成孔径雷达干涉测量中的应用研究［D］.合肥:安徽大学，2005.

［13］左臣婕.分布式小卫星SAR地面目标检测及速度估计［D］.北京:中国科学院电子学研究所，2007.

［14］Bone DJ.Fourier fringe analysis:The two dimensional phase unwrapping problem［J］.Applied Optics，1991，30(25):3627-3632.

［15］田晓迪.基于相位展开的磁敏感加权图像后处理方法的研究［D］.上海:上海理工大学，2009.

［16］Xu Wei，Cumming.I A region growing algorithm for InSAR phase unwrapping［J］.IEEE Transactions on Geoscience and Remote Sensing，1999，37(1):124-134.

Study on Processing Method for Phase Artifact in SWI Based on Phase Unwrapping Technique

NIE Sheng-Dong1*TIAN Xiao-Di1QIAN Li-Jun2CHEN Zeng-Ai2XU Jian-Rong2
2(Institute of Medical Imaging Engineering，University of Shanghai for Science and Technology，Shanghai 200093，China)2(Renji Hospital，Shanghai Jiaotong University，Shanghai 200127，China)

TP391

A

0258-8021(2010)04-0517-07

10.3969/j.issn.0258-8021.2010.04.007

2009-12-31，

2010-04-12

国家自然科学基金资助项目(60972122);上海市教委科研创新项目(09YZ216)

*通讯作者。 E-mail:nsd4647@sohu.com