本发明涉及一种医学图像处理方法，特别涉及一种基于数据一致性和图像伪影分解技术的运动伪影消除方法。
背景技术：
：：磁共振成像技术是利用氢质子的核磁共振现象进行成像的一种技术。人体内包含单数质子的原子核，例如广泛存在的氢原子核，其质子具有自旋运动。带电原子核的自旋运动，在物理上类似于单独的小磁体，而且在没有外部条件影响下这些小磁体的方向性分布是随机的。当人体置于外部磁场中时，这些小磁体将按照外部磁场的磁力线重新排列具体为在平行于或反平行于外在磁场磁力线的两个方向排列，将上述平行于外在磁场磁力线的方向称为正纵向轴，将上述反平行于外在磁场磁力线的方向称为负纵向轴，原子核只具有纵向磁化分量，该纵向磁化分量既具有方向又具有幅度。用特定频率的射频(rf，radiofrequency)脉冲激发处于外在磁场中的原子核，使这些原子核的自旋轴偏离正纵向轴或负纵向轴，产生共振，这就是磁共振现象。上述被激发原子核的自旋轴偏离正纵向轴或负纵向轴之后，原子核具有了横向磁化分量。停止发射射频脉冲后，被激发的原子核发射回波信号，将吸收的能量逐步以电磁波的形式释放出来，其相位和能级都恢复到激发前的状态，将原子核发射的回波信号经过空间编码等进一步处理即可重建图像。在磁共振扫描成像过程中，由于检测时间比较长，病人往往会有自主的或无意识的运动。这种运动会造成图像变模糊、产生伪影或者图像错配，影响医生对病人的诊断。文章“handbookofmripulsesequences”(bernstein,m.a,kingk.f,xiaohongjoezhou)指出消除运动伪影是磁共振成像领域研究的重点和技术难题之一。为了规避、矫正或者消除运动引起的图像质量下降，基本上有两种主要策略：由“reducingmotionartifactsintwo-dimensionalfouriertransformimaging.(magnresonimaging1986；4:359–376)”公开的前瞻性运动矫正和回顾性运动矫正技术。在前瞻性运动矫正技术，通过使用一定的方法，例如导航回波、呼吸触发等方法，实时的监控成像部位的空间位置，根据检测到的空间位置信息，实时的改变脉冲序列形状，尽可能的规避物体运动导致到图像误差。前瞻性运动矫正方法需要专门的硬件支持特别的位置监控信息，这些位置监控信息的精准度会对运动矫正效果产生重要影响；同时，对于非刚性物体运动，特别是高阶运动成分矫正效果不够明显。回顾性运动矫正技术中，物体的运动位置信息也会通过一定的技术，例如导航回波或者门控触发方法，记录下来。但并不会实时的改变扫描序列的形状。由于物体的图像以及对应的位置信息都被记录下来，因此，可以对生成的图像进行后处理，将物体的运动信息用于矫正初步生成的图像，进而消除由于运动导致的图像质量下降。1999年j.g.pipe提出了propeller(periodicallyrotatedoverlappingparallellineswithenhancedreconstruction)技术。该技术基于k空间非等间隔的数据采集，减少数据采集的时间，见图1；利用k空间中心数据带重叠区域的过采样数据估计运动形式并矫正运动伪影；运用网格化算法将非笛卡尔数据转化为笛卡尔数据，再进行傅立叶逆变换最终生成重建图像。这种方法对于刚性运动伪影的消除效果非常显著，已经在头颅磁共振成像中获得了成功应用。但是，这种方法主要是针对刚性运动物体，对于非刚性运动物体的运动，例如腹部、血管波动，吞咽等的矫正效果并不明显，而且由于要求需要必须过采样，延长了扫描时间。1993年，qing-sanxiang等人提出通过两次扫描，利用运动信息的相干性消除运动伪影的方法。这种方法中，包含运动伪影的图像(以下称为采集图像)可以被分解成不包含运动的复数图像(以下称为目标图像)和仅仅包含运动伪影的复数图像(以下称为伪影)。两次独立的扫描过程中，在扫描对象运动状态一致的情况下，修改图像采集的k空间(即数据空间)顺序，导致该两次独立扫描后，由于运动所导致的仅仅包含运动伪影的复数图像的相位发生了变化，而不包含运动伪影的目标图像没有影响，因此，只要能够估计出由于运动导致的伪影相位变化，即可通过解析的方法，直接求出伪影的大小和目标图像的大小。对于每个像素而言，由于运动导致的伪影的相位变化取值范围在0～2*pi之间。将整个成像范围划分为多个小的区域，不断调节该小区域内每个像素的相位，这时，每个小的区域对应的梯度能量就会发生变化，当该区域的梯度能量为最小时，就能够确认这个小区域内的相位最为准确。根据此时的伪影相位值，就能够直接求出伪影和目标图像的大小。这种基于“图像-伪影分解”的方法的优势在于，不需要专门的硬件检测成像物体的运动情况，而且不需要对物体的运动形式建模，刚性和非刚性运动都能适用；但是，这种方法对运动的形式要求太严格，必须保证两次独立采集中成像对象的运动模式基本一致，而这一点在实际是很难保证的；同时，这种方法需要改变k空间(数据空间)的填充顺序，而这种做法会明显的改变基于回波链序列的图像对比度。物体的运动在最终图像中的影响，可以通过运动的频率进行具体的分析。理论上，任何一种运动频率都会在图像中引入两个对称的伪影，该对称的伪影在图像中的位置由这个运动在k空间(数据空间)的频率(有效运动频率)决定。一般来说，有效运动频率越小，则由于运动而引入的伪影在图像中分散就不明显，反之，当有效运动频率越大，则由于运动而引入的伪影在图像中的分散就会更加明显。临床应用中的运动往往包含多种运动频率，而且每种频率的运动相互之间也可能增强或抵消，最终会形成非常复杂的伪影形状。实际运动频率虽然复杂，但是可以通过改变k空间的填充方式，改变运动在k空间的的有效运动频率，进而修改运动伪影在图像中的显示方式。1985年，bailes等人在“respiratoryorderedphaseencoding(rope)--amethodforreducingrespiratorymotionartifactsinmrimaging.(j.comput.assist.tomogr.9:835-838)”提出了一种“低频重排”方法。在这个方法中，运动位置相近的信号被尽量安排在相近的k空间中，这样，k空间中的有效运动频率非常低。通过这种方法，运动产生的伪影被限制在距离运动发生部位很小的位置；1986年，haaceke等人在“reducingmotionartifactsintwo-dimensionalfouriertransformimaging.(magnresonimaging1986；4:359–376)”提出了一种“中心重排”方法。这种方法尽量在呼吸运动的平稳期平稳期采集对图像对比度影响最大的k空间中心部分，而在运动较多的时期采集对比度影响较小的部分，同时尽量安排运动位置相近的信号在相近的k空间中。这种方法能够提供更为稳定的运动抑制效果；1998年，jhooti等人在“hybridorderedphaseencoding(hope):animprovedapproachforrespiratoryartifactreduction.(j.magn.reson.imaging8:968-980.)”提出了一种“高频重排”方法。在这个方法中，运动位置相差较大的信号被尽量安排在相近的k空间中，这样，k空间中的有效运动频率就非常高。通过这种方法，运动产生的伪影被扩散到了距离运动发生部位很远的位置。如果采用更大的成像视野fov，则这些伪影会被扩散到物体的外部。上述基于k空间重排的方法，能够一定程度上减轻伪影的程度，但是当运动模式不是理想情况，例如运动幅度、频率发生突然变化等情况下，不能起到预期的效果。在“simultaneousacquisitionofspatialharmonics(smash):ultra-fastimagingwithradiofrequencycoilarrays.(magnresonmed1997；38:591–603)”、“sensitivityencodingforfastmri.(magnresonmed1999；42:952–962)”、“generalizedsmashimaging.(magnresonmed2002；47:160–170)”、“generalizedautocalibratingpartiallyparallelacquisitions(grappa).(magnresonmed2002；47:1202–1210)”等中公开由于运动对图像造成的影响和运动持续的时间有直接的关系，因此，磁共振并行成像技术被逐渐应用于消除运动伪影。“artifactreductioninmoving-tableacquisitionsusingparallelimagingandmultipleaverages.(magnresonmed2007；57:226–232)”指出采用并行成像技术结合多次平均的方法，能够在一定程度上消除两次扫描中不相干的运动伪影。多通道相控阵采集情况下，考虑到并行成像技术提供的编码能力，正常采样的数据具有较大的冗余特性，基于这个特点，多个基于数据一致性和数据剔除的方案被应用于运动伪影的抑制和消除。2002年，m.bydder等人在“detectionandeliminationofmotionartifactsbyregenerationofk-space.(magnresonmed2002apr；47(4):677-86)”提出通过k空间数据再生技术检测和消除运动伪影的方法。这种方法将全采集的k空间数据分成奇数行和偶数行两个部分，分别对它们进行并行成像处理。分别检测特定位置采集数据和基于并行成像生成的数据的差异，根据这种差异的大小决定该信号是否符合总体一致性的要求，决定是剔除还是保留。方法“smashnavigatorsdavidatkinson,markbydder,josephv.hajnal,derekhill,davidlarkmanismrm2002p2395”采用k空间的并行成像技术，进一步拓展了信号的数据一致性。方法“care:coil-basedartifactreductionismrm2004abstractnumber96”使用最优化目标函数来描述信号的一致性。方法“dataconvolutionandcombinationoperation(cocoa)formotionghostartifactsreduction.(magnresonmed2010；64:157-166)”采用基于k空间的grappa并行成像技术，在数据一致性判定环节，基于信号差异的统计学指标，做出一致性判断。上述这些方法充分利用了并行成像技术来确保采集数据的一致性，但是，实际上并行成像的输入数据本身包含了运动伪影，导致数据一致性较差；同时，数据一致性的判定需要较为稳定的参考对象，实际中难以精确的得到。基于数据一致性消除运动伪影的方法在针对临床实际复杂运动情况下仍面临很多挑战。技术实现要素：本发明旨在提供一种基于数据一致性和图像-伪影分解的方法，一方面，该方法不依赖于特定的扫描序列，可以和多种其他伪影消除技术结合，适用性广泛；同时，采用基于并行成像的数据一致性方法，能够在消除运动伪影的同时较好的保持图像信噪比；最后，本发明对图像-伪影的分解不依赖于特定的运动模型，能够稳定、有效的消除临床实际中多种运动伪影。为达到上述目的，本发明采用的技术方案如下：本发明公开的基于数据一致性和图像伪影分解技术的运动伪影消除方法，包括以下步骤：a)采集多通道数据a；b)通过多通道数据获得自校准数据，生成nk组基于数据一致性的卷积核；c)根据该nk组卷积核，对采集的多通道数据a进行卷积，获得nk组满足数据一致性的多通道数据b；d)根据采集的多通道数据a，和新生成的nk组满足数据一致性的多通道数据b，计算出真实的图像及伪影；所述a)步骤中，通过多通道接收线圈，采集多通道数据a，数据大小表示为：nx*ny*nc，其中nx代表采集数据的行数，ny表示数据的列数，nc代表接收通道的数目；所述b)步骤中，自校准数据大小表示为：na*ny*nc，其中na为自校准数据的宽度，默认相位编码方向沿着行的方向，通过以下方程求解基于数据一致性的卷积核每一组卷积核的大小为四维数组：wx*wy*nc*nc，代表nc个wx*wy*nc三维卷积核，每一个通道对应一个三维卷积核，这样总的来看就是一个四维的卷积核，其中，j代表目标卷积核中包含的信号的数目，l代表接收通道数目，sl(kr)代表通道l中位置为kr的信号强度，sl′(kr-δkj)代表通道l′中位置为kr-δkj的信号强度，δkj为信号sl′相对于信号sl在该接收通道内的偏移量，卷积核权重是针对于通道l′的加权因子，方程(1)中sl(kr)和sl′(kr-δkj)为已知量；所述c)步骤中，将采集的多通道数据a带入方程(1)进行卷积，和sl′(kr-δkj)为已知量，通过卷积计算求出满足数据一致性的多通道数据b；所述d)步骤中，通过以下方程计算：i′0＝i0+g0(2)is＝i0+gs,s＝1,2,...,k(3)其中，i0表示不包含伪影的复数图像，g0代表由于运动导致的伪影，i′0代表直接采集后包含运动伪影的图像，gs是通过数据一致性处理后的伪影，is是通过数据一致性处理后包含伪影的数据，是运动伪影对基于数据一致性的卷积核的响应函数，as为幅度，θs为相位，在已知的情况下，求解上述(2)-(4)组成的方程组，就能够得到不包含伪影的复数图像i0，求解结果表示为(5)。优选的，所述d)步骤中，在求解方程(5)之前，需要确定卷积核对应的运动伪影的响应函数，即根据任一个输入的三维卷积核，即wx*wy*nc，其对应的响应函数计算方法由如下方程计算：其中，w大小为wx*wy*nc,ifft代表傅立叶逆变换算子。优选的，所述d)步骤中，在求解方程(5)的过程中，当该方程组的条件数大于50时，增加基于信噪比最优和伪影最小化的图像伪影消除步骤。进一步的，所述基于信噪比最优和伪影最小化的图像伪影消除：限定解ic的求解由以下方程计算：其中，ic是基于信噪比最优和伪影最小化的图像，ii是第i个基于数据一致性生成的多通道数据，wi是第i个数据的权重，是图像中像素位置周围范围内的标准差。进一步的，加入迭代过程，即将新生成的多通道数据认为是采集数据，重新执行b)、c)、d)过程，对伪影的抑制逐渐收敛。根据输入的原始数据，生成多组满足数据一致性的数据的目的是多方面的。一方面，基于数据一致性而新生成的数据中，包含的运动调制伪影是不同的，这些不同来自于运动伪影的特殊性，无法事先预知；第二，多组数据本身由于施加了数据一致性，因此，相对于原始采集数据，相对减少了运动伪影的影响(因为运动造成的伪影往往不满足数据一致性，而在重新生成的过程中受到抑制)。因此，在这里，满足数据一致性的数据生成方法都是可以包含进来的，例如sense,grappa,spirit等。方程(2)-(4)是一组超定方程，可以采用多种方法求解，例如广义逆矩阵法，共轭迭代法等。同时考虑到上述中提到的新生成满足数据一致性数据的特点，可以更好的对最终的求解做一定的限制。限制解的求出方法就是基于区域方差的多组数据合成。采用区域方差作为确定权重的原因在于，物理上伪影或者噪声较大的数据，其方差必然较大，因此，应该赋予较小的权重；而伪影或者噪声较小的数据，其方差必然较小，应该赋予较大的权重。类似的，可以采用其他基于图像特征的数据量来决定权重的计算，例如可以采用图像的全变分或者熵等。本发明具有以下有益效果：1.本发明不依赖于特定的扫描序列和重建方法，因此，具有广泛的适用性，能够与现有临床的多种扫描序列及特定重建方法相结合；2.本发明基于数据的一致性生成最终图像，不仅能够较好的抑制运动伪影，而且较好的保证了最终图像的信噪比；3.本发明对图像-伪影的分解不依赖于特定的运动模型，能够适用于刚性和非刚性运动，以及高阶运动模式，因此，能够稳定、有效的消除临床多变的扫描条件下的运动伪影。4.通过迭代过程，能够对伪影的抑制逐渐收敛。附图说明图1为基于k空间非等间隔的数据采集，减少数据采集时间的示意图；图2为基于数据一致性和鬼影消除的运动伪影消除方法流程图；图3为基于信噪比最优和伪影最小化的图像合成；图1中m0(k)代表k空间采集线所对应的物体偏移位置，横轴代表k空间的采集位置，纵轴代表采集该k空间数据的时候对应的物体偏移位置；m1(k)代表经过卷积处理后，k空间采集线所对应的物体偏移位置；图3中，区域的大小为kr*kr，为了求出区域内黑色点位置的权重wi，选择上述区域内的所有点，求出这些点对应的标准差代入方程(7)(8)即可。具体实施方式为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。如图2所示，本发明空开的基于数据一致性和图像伪影分解技术的运动伪影消除方法，包括以下步骤：a)采集多通道数据。采用多通道接收线圈，接收数据。该采样数据可以来自于多种扫描序列。该数据大小可以表示为：nx*ny*nc。其中，nx代表采集数据的行数，ny表示数据的列数，nc代表接收通道的数目。如图1m0(k)中所示，由于物体的运动，导致k空间数据采集的过程中物体的偏移位置发生了变化，这种变化等价于在k空间对数据做了运动的调制，这也是图像中最终出现运动伪影的根本原因。b)基于上述多个通道数据中的自校准数据，生成nk组基于数据一致性的卷积核。自校准数据来自于上述步骤中采集到的多通道数据的中心位置，数据大小可以表示为：na*ny*nc。其中na为自校准数据的宽度，默认相位编码方向沿着行的方向。基于数据一致性的卷积核可通过以下方程求解：其中，j代表目标卷积核中包含的信号的数目，l代表接收通道数目，sl(kr)代表通道l中位置为kr的信号强度，sl′(kr-δkj)代表通道l′中位置为kr-δkj的信号强度，δkj为信号sl′相对于信号sl在该接收通道内的偏移量，卷积核权重是针对于通道l′的加权因子。每一组卷积核的大小为四维数组：wx*wy*nc*nc。在该步骤中，方程(1)中sl(kr)和sl′(kr-δkj)被认为是已知量，卷积核为未知量，通过求解上述线性方程能够计算出该卷积核。c)根据该nk组卷积核，对采集的多通道数据进行卷积，获得nk组满足数据一致性的多通道数据。生成满足数据一致性的方法仍然满足上述方程(1)。在该步骤中，和sl′(kr-δkj)均为已知量，通过卷积计算就能求出满足数据一致性的多通道数据。d)根据采集的多通道数据，和新生成的nk组多通道数据，计算出真实的图像及伪影，可用如下方程计算：i′0＝i0+g0(2)is＝i0+gs,s＝1,2,...,k(3)其中，i0表示不包含伪影的复数图像，g0代表由于运动导致的伪影，i′0代表直接采集后包含运动伪影的图像，gs是通过数据一致性处理后的伪影，is是通过数据一致性处理后包含伪影的数据，是运动伪影对基于数据一致性的卷积核的响应函数，as为幅度，θs为相位。在已知的情况下，求解上述(2)-(4)组成的方程组，就能够得到不包含伪影的复数图像i0。求解结果表示如下：e)求解上述方程之前，需要确定卷积核对应的运动伪影的响应函数，即根据任一个输入的三维卷积核，其对应的响应函数计算方法可用如下方程计算：其中，w是大小为wx*wy*nc,ifft代表傅立叶逆变换算子。f)基于信噪比最优和伪影最小化的图像伪影消除。求解方程组(5)的过程中，如果该方程组的条件数大于50时，则该矩阵接近奇异矩阵，求解的结果有较大的偏差。因此，需要对求出的结果进行限定。限定解ic的求解可以用以下方程计算：其中，ic是基于信噪比最优和伪影最小化的图像，ii是第i个基于数据一致性生成的多通道数据，wi是第i个数据的权重，是图像中像素位置周围范围内的标准差。g)上述求解多组图像合成的过程，输入为s组大小为nx*ny*nc的图像，输出是一组大小为nx*ny*nc的图像。其中，合成后的每个像素都有来自于s组结果的成分，成分的大小来自于合成公式(7)中的权重大小。如图3所示，每一个输入图像像素的权重，依赖于该像素周围区域内图像的方差：方差越大则对应的权重越小，方差越小则对应的权重越大。该合成方法的物理意义在于，当输入图像区域含有较多伪影或噪声时，则该图像区域的方差较大，此时应该具有较小的权重；反之，当输入图像区域含有的伪影或噪声较少，则该图像区域的方差较小，此时应该具有较大的权重。这种求解方法中，平等的对待运动造成的伪影和各种信号噪声，能够对两者同时进行优化。上述步骤f)-g)中，基于区域支持的多组图像合成方法中，kr*kr构成的区域形成了对该区域中心点的支持区域，该区域的图像的方差，决定了合成过程中所需要的加权权重。h)加入迭代过程，即将新生成的多通道数据认为是采集数据，重新执行b)-g)过程，对伪影的抑制逐渐收敛。当然，本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。当前第1页12当前第1页12