结合k空间和图像空间重建的高鲁棒快速成像方法和装置

1.本发明属于磁共振成像领域，可以针对需要多帧图像数据的磁共振参数测量 实现高鲁棒快速成像。


背景技术：

2.磁共振成像(mri)可用来测量各种生理参数，已被证明其可以用来定量评 估各种病理信息，例如纵向弛豫时间(t1)，横向弛豫时间(t2)，表观扩散系数 (adc)，酰胺质子转移(apt)等。定量测量生理参数需要采集具有脉冲序列 参数调制的多部数据集，例如回波时间(te)、翻转角度(fa)、反转时间(ti) 或饱和频率。然后将多部分数据集中的信号演化按照参数模型进行拟合以实现逐 像素的参数估计。由于多次数据采集而导致mri参数测量扫描时间长，这限制 了其在临床上的应用。较长的采集时间不仅会降低患者吞吐量，还可能由于患者 的运动而影响参数定量的效果。
3.对其扫描时间较长这一局限，已有很多研究者开发出相关的方法来加速多帧 成像。以t1、t2定量为例，有压缩感知(cs)，repcom，llr-spirit，k-t slr， k-t pca，磁共振指纹(mrf)，aloha等方法。压缩感知利用了图像在特定变 换域中的稀疏性，例如小波变换等，一些技术已成功应用于mr参数映射。而由 于时间依赖性可以建模为时空数据的低秩性，因此低秩性通常作为优化问题的强 制约束。k-t pca利用矩阵分解将训练数据分解为一组主成分，从中估计参数。 mrf使用基于仿真的字典和模型匹配算法，以同时估计多个参数映射。
4.并行成像是在临床上应用最广的加速方法，其中又以grappa和sense最 为广泛。grappa将k空间中的缺失点建模所有通道的相邻的获得数据点的线性 组合；sense将重建问题视为图像空间的线性求逆。理论上来说，如果获得的 灵敏度图是准确的，sense方法是最优的，可获得最小二乘意义上的精确解。 但是灵敏度图通常是不准确的。最近，一种应用在多帧成像领域的变加速敏感度 编码方法(vsense)被提出以加速参数定量。通过假设在多帧成像中不同帧的 灵敏度图是一致的，vsense可以实现较高倍数的加速。作为对vsense方法的 改进，kipi方法实现了自动校准，不需要采集额外的灵敏度图，也无需全采样 帧，且允许更为灵活的变提速策略。但是kipi方法获取的灵敏度图并不总是准 确的，特别在图像的信号较低区域，例如反转恢复(ir)弛豫时间映射序列中的头 皮区域。这些地方的错误将会导致信号在图像中混叠，造成伪影。
5.因此，在多帧成像领域，提出一种能够自动优化灵敏度图的方法具有重要意 义。


技术实现要素：

6.本发明的目的在于解决多帧成像技术扫描时间较长，现有加速技术的线圈灵 敏图不准确等问题，并提供一种结合k空间和图像空间重建的高鲁棒快速成像方 法(英文名字为joint k-space and image-space parallel imaging with highrobustness，后续可简称为robustkipi方法)和装置，可实现自动优化灵敏度 图，从高度欠采样数据中精确
重建图像。
7.本发明所采用的具体技术方案如下：
8.第一方面，本发明提供了一种结合k空间和图像空间重建的高鲁棒快速成像 方法，用于多帧成像中欠采样图像帧进行重建，所述欠采样图像帧包含带有自动 校准信号的第一欠采样图像帧和加速因子不低于2的第二欠采样图像帧，且第一 欠采样图像帧的加速因子不高于第二欠采样图像帧的加速因子；
9.所述成像方法步骤如下：
10.s1：针对第一欠采样图像帧，利用k空间自动校准的重建方法重建得到完整 图像帧作为校准帧；对校准帧的k空间数据进行傅里叶变换后得到每个通道的通 道图像，将所有通道图像进行合并得到通道合并图像；
11.s2：对第一欠采样图像帧的自动校准信号使用espirit方法，得到espirit 线圈灵敏度图和espirit特征值图，基于espirit特征值图识别得到完全覆盖 目标成像区域的第一支撑区域；将每个通道图像除以所述通道合并图像，得到每 个通道的原始线圈灵敏度图；基于所述通道合并图像识别得到所述第一支撑区域 内去除噪声和空洞区域后的第二支撑区域；所述espirit线圈灵敏度图和原始 线圈灵敏度图均由对应的幅值图和相位图组成；
12.s3：将处于第一支撑区域但不在第二支撑区域内的espirit线圈灵敏度幅 值图与第二支撑区域内的原始线圈灵敏度幅值图相加，并通过平滑滤波后得到优 化的灵敏度幅值图；基于第二支撑区域内的原始线圈灵敏度相位图对所述噪声和 空洞区域进行插值补全，得到优化的灵敏度相位图；将优化的灵敏度幅值图和优 化的灵敏度相位图相结合，最终得到优化的灵敏度图；
13.s4：利用优化的灵敏度图，在所述校准帧上进行加速因子为1的sense重 建，生成参考图像；同时，对所述校准帧进行回顾性欠采样且回顾性欠采样的加 速因子与所述第二欠采样图像帧的加速因子相同，再对回顾性欠采样后的数据进 行sense重建，生成具有潜在混叠伪影的回顾性重建图像；基于所述参考图像 和所述回顾性重建图像得到校正因子图，且校正因子图中每个位置的像素值为参 考图像和回顾性重建图像中对应位置像素值的商；
14.s5：针对除第一欠采样图像帧之外的其余第二欠采样图像帧，利用所述优化 的灵敏度图进行sense重建，将sense重建后的图像与所述校正因子图相乘， 得到抑制伪影的完整图像帧。
15.作为第一方面的优选，所述第一欠采样图像帧和第二欠采样图像帧均为二维 图像；所述第一欠采样图像帧的加速因子优选为2，所述第二欠采样图像帧的加 速因子优选为2到4。
16.作为第一方面的优选，所述第一欠采样图像帧和第二欠采样图像帧均为三维 图像；所述第一欠采样图像帧的加速因子优选为4，所述第二欠采样图像帧的加 速因子优选为4到12。
17.作为第一方面的优选，所述s1中，所述k空间自动校准的重建方法优选为 grappa方法，通过计算出线圈权重并应用于欠采样区域，重建得到完整图像帧； 优选的，所述grappa方法为带有tikhonov正则化的grappa。
18.作为第一方面的优选，所述s1中，通过平方根方法或者自适应方法合并所 有的线圈图像。
19.作为第一方面的优选，所述s2中，espirit线圈灵敏度图和espirit特征 值图均为espirit方法特征值分解后的最大特征值对应的一套espirit线圈灵 敏度图和espirit特征值图(每个通道各自有一个espirit线圈灵敏度图，但 所有通道对应于1个espirit特征值图)。
20.作为第一方面的优选，所述s4中，所述校正因子图需经过滤波器滤波，以 去除异常值。
21.作为第一方面的优选，所述s5中，利用所述优化的灵敏度图进行sense 重建时，采用截断奇异值正则化方法。
22.第二方面，本发明提供了一种磁共振成像数据处理装置，其包括存储器和处 理器；
23.所述存储器，用于存储计算机程序；
24.所述处理器，用于当执行所述计算机程序时，实现如第一方面任一项所述的 结合k空间和图像空间重建的高鲁棒快速成像方法。
25.第三方面，本发明提供了一种磁共振成像装置，其包括磁共振扫描器以及控 制单元，所述控制单元中存储有计算机程序，当所述计算机程序被执行时，用于 实现如第一方面任一项所述的结合k空间和图像空间重建的高鲁棒快速成像方 法；所述磁共振扫描器用于获取所述第一欠采样图像帧和第二欠采样图像帧数据。
26.本发明相对于现有技术而言，具有以下有益效果：
27.现有技术需要精确识别噪声区域和信号区域，并且需要大幅外推灵敏度图， 这通常会造成灵敏度图的出现较高的异常值。本发明可通过espirit方法产生 的线圈灵敏度图自动补全空洞或信号周围区域的线圈灵敏度图，实现自动优化。
28.现有技术受限于数据信噪比，当信噪比较差时，线圈灵敏度图也较差。本发 明通过引入espirit方法，将直接相除获得的灵敏度图和通过模型特征值分解 获得的灵敏度图相结合，产生更加鲁棒和精确的线圈灵敏度图。
29.相比于现有技术，本发明在保证了加速效果的同时，提升了图像重建效果， 所重建的图像与真实值一致，且没有明显伪影。
附图说明
30.图1为rubustkipi组合原始线圈灵敏度幅值图和espirit线圈灵敏度幅值 图的流程图。第一支撑区域a为最终要形成的灵敏度范围；第二支撑区域b为 第一支撑区域a中排除了原始灵敏度图的噪声和空洞区域之后的剩余区域，a 区域和b区域的差值部分即为噪声和空洞区域c。组合的线圈灵敏度幅值图是b 区域内的原始灵敏度幅值图和c区域内的espirit线圈灵敏度幅值图相结合。
31.图2为三种不同的线圈灵敏度图(通道7)的幅值和相位。其中，左列:原始 线圈灵敏度图，将每个通道的图像除以均方根图像。中间:用espirit计算的线 圈灵敏度图。右图:将原始和espirit相结合，得到组合的灵敏度图并进行拟合。 白色箭头表示由于噪声或空洞导致的而错误灵敏度图区域。
32.图3为从欠采样ir tse数据中得到的重建的水模源图像和对应的t1参数图 像；其中(a):使用grappa、espirit、kipi和rubustkipi从4倍欠采样数据重 建ir tse图像(ti＝300ms)。kipi分别使用原始(raw)灵敏度图和espirit灵 敏度图。(b):全采样真实值与重建图像的误差，底部为nrmse。(c-d):t1参数图 像(c)和重建误差(d)。(a-b)中箭头指示源图像中的伪影，(c)中白色箭头指 示t1参数图像中的伪影。其中源图像误差图像放大20倍，t1参数图像误差图放 大10倍。nrmse:归一化均方根误差。
33.图4为使用grappa、espirit、kipi和rubustkipi进行重建的ti＝50ms 的ir-tse图像(a)。第二行(b)显示与全采样真实值的差别。kipi分别使用原始 (raw)灵敏度图和espirit灵敏度图。灰色箭头指示grappa和espirit结 果的混叠伪影。白色箭头指示相对较大的误差，源图像误差图像为放大20倍的 结果。nrmse值显示在底部。
34.图5为从图4所示的源图像计算得到的t1参数图像(a)。第二行(b)显示 显示与全采样真实值的差别。kipi分别使用原始(raw)灵敏度图和espirit 灵敏度图。白色箭头指示使用raw和espirit线圈灵敏度图的kipi中相对较大 的伪影，t1参数图像误差图为放大10倍的结果。nrmse值显示在底部。
35.图6为本发明rubustkipi方法的算法流程图，其中完成重建的第一欠采样 图像帧即为校准帧。
具体实施方式
36.下面结合附图和具体实施方式对本发明做进一步阐述和说明。本发明中各个 实施方式的技术特征在没有相互冲突的前提下，均可进行相应组合。
37.在需要多帧成像的磁共振参数测量领域，数据集由从相同体积中获取的一系 列不同脉冲序列参数的图像组成，这意味着不同帧的线圈灵敏度图理论上是相同 的。传统的sense方法显式地使用来自单独参考扫描的灵敏度图，而grappa 方法通过自动校准信号(acs)估计的k空间相关性隐式地利用灵敏度编码。然而， 在sense重建中很难获得精确的灵敏度图，这限制了其在实际中的最大加速度 因子。对于grapaa，所有图像帧都需要重复获取acs数据，导致扫描时间相 对延长。
38.最近提出的kipi方法利用不同cest帧的可变加速因子(af)，并结合了 sense和grappa的优点。kipi以低加速度因子(如af＝2)对cest帧进行采样， 使用grappa填充未采样的k空间，从grappa生成的全k空间获得灵敏度图， 然后使用sense和伪像抑制算法重建其他高度欠采样(af》2)帧。kipi每个通道 的图像除以通道组合图像得到的线圈灵敏度图并不总是准确的，尤其是在某些序 列中会出现空洞区域，例如反转恢复序列中的头皮区域；在信噪比较低的应用中， kipi获得的灵敏度图也相对不准确。
39.espirit方法通过对构造的校准矩阵a进行svd分解，通过截断较小的奇 异值，获得子空间和零空间。espirit模型将重建问题视作子空间补全问题，最 终可对于每个像素位置都有一个小的校准方程，如下：
40.ψqcsmq＝csmqꢀꢀꢀ
1.41.其中下标q代表位置，为该处的相关矩阵，为该像 素位置处的线圈灵敏度向量，nc为通道数量。因此线圈灵敏度向量csmq是相关 矩阵ψq的特征值为1对应的特征向量。对ψq进行特征值分解，可以得到nc个 特征值和对应的nc个特征向量，所
有特征值不会大于1，因此有nc套线圈灵敏 度图，但是特征值最大的那套相对最准确。迭代所有像素位置，可以得到模型分 解出的espirit线圈灵敏度图。其具有光滑，且噪声不敏感等特点，但是espirit 线圈灵敏度图相对比较模糊且相位信息不准确。
42.本发明通过结合kipi灵敏度图和espirit线圈灵敏度图，得到更加准确更 加鲁棒的灵敏度图以供重建使用。本发明提供的结合k空间和图像空间重建的高 鲁棒快速成像方法流程如图6所示，下面对其具体实现过程进行展开描述。
43.该成像方法针对磁共振机器采集到的欠采样图像帧进行重建。其中分为第一 欠采样图像帧，以较低的加速因子欠采样并包含acs；和第二欠采样图像帧， 以较高的加速因子(af》2)欠采样。第二欠采样图像帧的加速因子不低于第一 欠采样图像帧。第一欠采样图像帧只有一帧，其余皆为第二欠采样图像帧。利用 第一欠采样图像帧可以获得准确的灵敏度图和校正因子图，可对较高加速因子欠 采样的第二欠采样图像帧进行精确重建，实现了图像质量和采集速度的平衡。
44.需要说明的是，磁共振机器采集的多帧图像数据可以是2d数据也可以是3d 数据，变加速因子根据实际情况进行改变。若为2d图像，所述第一欠采样图像 帧的加速因子优选为2，所述第二欠采样图像帧的加速因子优选为2到4；若为 3d图像，所述第一欠采样图像帧的加速因子优选为4，所述第二欠采样图像帧 的加速因子优选为4到12。
45.该结合k空间和图像空间重建的高鲁棒快速成像方法，具体步骤如下：
46.s1：针对第一欠采样图像帧，利用k空间自动校准的重建方法重建得到完整 图像帧作为校准帧；对校准帧的k空间数据进行傅里叶变换后得到每个通道的通 道图像，将所有通道图像进行合并得到通道合并图像。
47.本发明中，k空间自动校准的重建方法可以采用grappa方法，其通过计算 出线圈权重并应用于欠采样区域，重建得到完整图像帧。具体而言，grappa 属于现有技术，其通过acs计算出权重将各通道数据线性组合重建完整图像帧。 进一步，若使用grappa方法，本发明推荐采用带有tikhonov正则化的grappa。
48.但需要说明的是，在本步骤中，虽然推荐使用grappa来实现第一欠采样 图像帧的重建，但是根据情况也可以采取sake，spirit，loraks等其他的k 空间中自校准的并行成像方法。
49.另外，在本步骤中，所有的线圈图像合并可以通过平方根(rss)方法实现， 当然也可以采用自适应合并(adaptive combine)方法实现。
50.s2：对第一欠采样图像帧的自动校准信号使用espirit方法，得到espirit 线圈灵敏度图和espirit特征值图，基于espirit特征值图识别得到完全覆盖 目标成像区域的第一支撑区域；将每个通道图像除以所述通道合并图像，得到每 个通道的原始线圈灵敏度图；基于所述通道合并图像识别得到所述第一支撑区域 内去除噪声和空洞区域后的第二支撑区域。
51.需要说明的是，上述espirit线圈灵敏度图和原始线圈灵敏度图均为复数 形式，可以分解为对应的幅值图和相位图。espirit线圈灵敏度图由espirit线 圈灵敏度幅值图和espirit线圈灵敏度相位图组成，原始线圈灵敏度图由原始 线圈灵敏度幅值图和原始线圈灵敏度相位图组成。
52.由于本发明中不考虑成像物体大于fov的情况，espirit线圈灵敏度图和espirit
特征值图均取espirit方法特征值分解后的最大特征值对应的一套线圈 灵敏度图和特征值图。
53.需要说明的是，第一支撑区域相当于最终要形成的目标灵敏度图的范围区域， 其应当完全覆盖目标成像区域。例如对于脑部磁共振成像而言，目标成像区域即 为大脑区域。但是第一支撑区域在覆盖目标成像区域时，边缘应当留有余量避免 边缘处产生成像噪声误差。所谓边缘应当留有余量，即第一支撑区域的边缘相对 于目标成像区域的边缘略微外扩，一般外扩5～20像素即可。且根据两个支撑区 域的识别过程可知，第二支撑区域应当完全被包含于第一支撑区域内。
54.另外，在本步骤中，第一支撑区域和第二支撑区域优选为阈值识别方法。所 谓阈值识别方法是指通过设置一个阈值来对大于阈值的区域进行提取作为相应 的支撑区域。用于提取第一支撑区域和第二支撑区域的具体阈值可根据实际进行 调整，在本发明中均可通过以下阈值识别方法来实现第一支撑区域和第二支撑区 域的识别提取：
55.1)第一支撑区域的识别提取：对通道合并图像使用otsu方法(最大类间方 差法)确定otsu区域，计算此otsu区域在整幅通道合并图像中的占比为a1，将 espirit特征值图中的像素值从大到小进行排序，选择百分位为a1的像素值作 为第一阈值，将espirit特征值图中像素值大于第一阈值的区域识别为第一支 撑区域；
56.2)第二支撑区域的识别提取：从通道合并图像中找到最大像素值q，然后 根据预设的比例值b1，计算第二阈值为q*b1，将通道合并图像中大于第二阈值 的区域识别为第二支撑区域。上述b一般可取值为0.1～0.2，以能够尽可能识别 出灵敏度图内的噪声和空洞区域为准。
57.s3：将处于第一支撑区域但不在第二支撑区域内的espirit线圈灵敏度幅 值图(即espirit线圈灵敏度幅值图中处于第一支撑区域但不在第二支撑区域 内的部分)与第二支撑区域内的原始线圈灵敏度幅值图(即原始线圈灵敏度幅值 图中位于第二支撑区域内部分)相加，并通过平滑滤波后得到优化的灵敏度幅值 图；基于第二支撑区域内的原始线圈灵敏度相位图(即原始线圈灵敏度相位图中 位于第二支撑区域内部分)对所述噪声和空洞区域进行插值补全，得到优化的灵 敏度相位图；将优化的灵敏度幅值图和优化的灵敏度相位图相结合，最终得到优 化的灵敏度图。
58.需要说明的是，本发明中基于第二支撑区域内的原始线圈灵敏度相位图对所 述噪声和空洞区域进行插值补全，可以采用多种不同形式的插值算法，例如局部 加权多项式回归算法，只要能够实现利用有值像素来对空值像素进行插值补全即 可。
59.在本步骤中，将处于第一支撑区域但不在第二支撑区域内的espirit线圈 灵敏度幅值图与第二支撑区域内的原始线圈灵敏度幅值图相加是指在对应像素 位置相加。优化的灵敏度幅值图和优化的灵敏度相位图相结合是指组合为优化的 灵敏度图(复数)，而且可以使用拟合以进一步平滑。
60.s4：利用优化的灵敏度图，在所述校准帧上进行加速因子为1的sense重 建，生成参考图像；同时，对所述校准帧进行回顾性欠采样且回顾性欠采样的加 速因子与所述第二欠采样图像帧的加速因子相同，再对回顾性欠采样后的数据进 行sense重建，生成具有潜在混叠伪影的回顾性重建图像；基于所述参考图像 和所述回顾性重建图像得到校正因子图，且校正因子图中每个位置的像素值为参 考图像和回顾性重建图像中对应位置像素值
的商。
61.需要注意的是，直接相除得到的校正因子图可能存在异常值，因此最好将其 经过滤波器进行滤波去除异常值后，供下一步使用。
62.需要注意的是，本发明中第二欠采样图像帧有多帧，且加速因子可能不同， 因此需要对每个用到的加速因子都产生一个对应的校正因子图。例如，如果一帧 第二欠采样图像帧的加速因子为x，那么在s4中需要也基于加速因子x对校准 帧进行回顾性欠采样，然后通过sense重建回顾性重建图像，并基于参考图像 和回顾性重建图像得到加速因子x对应的校正因子图；当另外一帧第二欠采样 图像帧的加速因子为y时，那么在s4中需要也基于加速因子y对校准帧进行回 顾性欠采样，然后通过sense重建回顾性重建图像，并基于参考图像和回顾性 重建图像得到加速因子y对应的校正因子图。
63.s5：针对除第一欠采样图像帧之外的其余第二欠采样图像帧，利用所述优化 的灵敏度图进行sense重建，将sense重建后的图像与所述对应校正因子图 相乘，得到抑制伪影的完整图像帧。
64.需要说明的是，本步骤中sense重建时，优选选用截断奇异值正则化方法， 丢弃部分较小的奇异值。
65.需要说明的是，本步骤中将sense重建后的图像与所述对应校正因子图相 乘是指对应位置上的两个像素值相乘，作为完整图像帧中对应位置的像素值。
66.需要注意的是，本步骤中将sense重建后的图像与所述对应校正因子图相 乘时，应该使用对应加速因子的校正因子图，例如sense重建后的图像是从加 速因子为x的欠采样数据重建的，就要用加速因子x对应的校正因子图与其相 乘；若sense重建后的图像是从加速因子为y的欠采样数据重建的，就要用加 速因子y对应的校正因子图与其相乘。
67.由此，上述s1～s4即构成了本发明的结合k空间和图像空间重建的高鲁棒 快速成像方法(robustkipi方法)。
68.同样的，基于同一发明构思，本发明的另一较佳实施例中还提供了与上述实 施例提供的结合k空间和图像空间重建的高鲁棒快速成像方法对应的一种磁共 振成像数据处理装置，其包括存储器和处理器；
69.所述存储器，用于存储计算机程序；
70.所述处理器，用于当执行所述计算机程序时，实现如前述s1～s5所示的结 合k空间和图像空间重建的高鲁棒快速成像方法。
71.可以理解的是，上述存储器可以包括随机存取存储器(random accessmemory，ram)，也可以包括非易失性存储器(non-volatile memory，nvm)，例 如至少一个磁盘存储器。同时存储介质还可以是u盘、移动硬盘、磁碟或者光 盘等各种可以存储程序代码的介质。当然，随着云服务器的广泛应用，上述软件 程序也可以搭载于云平台上，提供相应的服务，因此计算机可读存储介质并不限 于本地硬件的形式。
72.可以理解的是，上述的处理器可以是通用处理器，包括中央处理器(centralprocessing unit，cpu)、网络处理器(network processor，np)等；还可以是数字 信号处理器(digital signalprocessing，dsp)、专用集成电路(application specificintegrated circuit，asic)、现场可编程门阵列(field-programmable gate array， fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。
73.另外需要说明的是，所属领域的技术人员可以清楚地了解到，为描述的方便 和简洁，上述描述的装置的具体工作过程，可以参考前述方法实施例中的对应过 程，在此不再赘述。在本技术所提供的各实施例中，所述的装置和方法中对于步 骤或者模块的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方 式，例如多个模块或步骤可以结合或者可以集成到一起，一个模块或者步骤亦可 进行拆分。
74.此外，上述的存储器中的逻辑指令可以通过软件功能单元的形式实现并作为 独立的产品销售或使用时，可以存储在一个计算机可读取存储器中。基于这样的 理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方 案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储器 中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者 网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。
75.同样的，基于同一发明构思，本发明的另一较佳实施例中还提供了与上述实 施例提供的结合k空间和图像空间重建的高鲁棒快速成像方法对应的一种磁共 振成像设备，其包括磁共振扫描器以及控制单元。其中：
76.所述控制单元中存储有计算机程序，当所述计算机程序被执行时，用于实现 如前述s1～s5所示的结合k空间和图像空间重建的高鲁棒快速成像方法；所述 磁共振扫描器用于获取所述第一欠采样图像帧和第二欠采样图像帧数据。
77.上述磁共振扫描器可采用现有技术实现，可采用成熟商用产品，具体型号不 限。
78.需要说明的是，上述控制单元中除了存储有实现robustkipi方法的程序之 外，还应当具有实现多帧成像所必要的成像序列以及其他软件程序。
79.当然，上述控制单元可以是独立的控制单元，也可以是磁共振扫描器自带的 控制单元，即上述结合k空间和图像空间重建的高鲁棒快速成像方法可以以数据 处理程序的形式集成于磁共振成像设备的控制单元中，使得无需额外新增控制单 元即可由磁共振扫描器直接输出重建结果。
80.下面将通过一个实施例来进一步展示本发明前述结合k空间和图像空间重 建的高鲁棒快速成像方法(robustkipi方法)所能够实现的技术效果，以便本领 域技术人员更好地理解本发明的实质。
81.实施例
82.本实施例将上述s1～s5所示的结合k空间和图像空间重建的高鲁棒快速成 像方法(robustkipi方法)应用于具体实例中，下面详细描述robustkipi方法 的具体数据来源以及各步骤的具体实现方式如下：
83.1、mri实验
84.水模和人体实验在一台3t西门子磁共振扫描仪(magnetom prisma)上 进行，使用20通道头部接收线圈。水模使用西门子提供的标准均匀水模。人体 实验得到了当地机构审查委员会的批准。所使用的序列为反转恢复(ir)tse成 像序列，测量的mri参数为纵向弛豫时间(t1)。
85.对于人体和水模实验，序列参数如下：tr＝3000ms，te＝7.5ms,fov＝ 220
×
220mm2,分辨率＝1.7
×
1.7mm2，层厚＝5mm，采集矩阵＝128
×
128，涡 轮系数＝16。采集了7个不同的反转时间的全采样数据，ti＝50,150,300,500, 800,1300,2000ms。
86.2、图像重建与分析
87.所有的处理和分析都是使用在pc电脑(3.2ghz)上编写的matlab (mathworks,natick,ma)软件离线进行的。
88.对于水模，选择第七帧(ti＝2000ms)为第一欠采样图像帧，af＝2。其他六 帧为第二欠采样帧，af＝4。对于第一欠采样图像帧保留了中间24条相位编码线 作为acs，第二欠采样图像帧则没有acs。以上述方式对全采样数据进行回溯 性欠采样，得到变加速因子的数据。对于人体实验，也采取上述回溯性欠采样方 式。
89.本实施例的robustkipi分为以下几步：
90.第一步：对第一欠采样图像帧进行grappa重建，并且将acs填入，得到 校准帧。对校准帧进行傅里叶变换，得到每个通道的通道图像；将所有通道图像 使用平方根(rss)合并方法得到合并后的rss图像。其中grappa内核大小 为2
×
5，即2条采集到的相位编码线和5个频率编码点，每个点使用10
×
20(通 道数目)个点来拟合。
91.第二步，计算espirit、原始灵敏度图；第一、第二支撑区域。对第一欠采 样图像帧的自动校准信号使用espirit方法，得到espirit线圈灵敏度图和 espirit特征值图。其中espirit的内核大小为6
×
6；子空间阈值为0.02，即 大于最大奇异值的2％的奇异值对应的右奇异向量才列为子空间的基；第一支撑 区域通过对espirit特征值图实施阈值识别操作确认，本例中阈值为0.99。 espirit线圈灵敏度图和espirit特征值图均为espirit方法特征值分解后的最 大特征值对应的一套线圈灵敏度图和特征值图，无需多套。
92.将每个通道图像除以rss图像，得到每个通道的原始线圈灵敏度图。对rss 图像进行阈值操作识别出第二支撑区域，本例中阈值为0.1。第二支撑区域含于 第一支撑区域。
93.第三步，计算优化的灵敏度图。将处于第一支撑区域但不在第二支撑区域内 的espirit线圈灵敏度幅值图与第二支撑区域内的原始线圈灵敏度幅值图相加， 得到优化的灵敏度幅值图。由于espirit线圈灵敏度的相位与原始相位有较大 差距，因此单独对第二支撑区域内的原始线圈灵敏度相位图进行外推插值得到优 化的灵敏度相位图；将优化的灵敏度幅值图和优化的灵敏度相位图组合为复数灵 敏度图后，再进行拟合，得到最终的优化的灵敏度图。在本例中，外推插值和拟 合都使用一个带有三次加权核的局部加权多项式回归(lwpr)操作进行。本实施 例采用的多项式为二阶多项式，窗宽为12。
94.第四步，计算用来伪影抑制的校正因子图。利用优化的灵敏度图，在校准帧 上进行加速因子为1的sense重建，生成参考图像ρ0；同时，对所述校准帧进 行回顾性欠采样(回顾性欠采样的加速因子为2或4)，再对欠采样后的数据进 行sense重建，生成具有潜在混叠伪影的回顾性重建图像ρ1；基于参考图像和 回顾性重建图像得到校正因子图，且校正因子图中每个位置的像素值为参考图像 和回顾性重建图像中对应位置像素值的商。此外，本实施例中的校正因子图需要 被一个窗口为3
×
3的中值滤波器进行滤波，以去除异常值。校正因子图c可以 表述如下：
95.c＝ρ0./ρ196.在本例中，由于第二欠采样图像帧存在两个加速因子(af＝2或4)，因此这个过 程要进行两次，以得到af＝2对应的校正因子图和af＝4对应的校正因子图。
97.第五步，重建所有第二欠采样图像帧。利用优化的灵敏度图对第二欠采样图 像帧进行sense重建得到具有潜在伪影的重建图像，并使用截断奇异值正则化， 丢弃小于最大
值2％的奇异值。随后进行伪影抑制，将第二欠采样图像帧与对应 的校正因子图相乘，得到最终的重建好的抑制伪影的图像。
98.上述第一步到第五步过程可参见图6所示的算法流程图。
99.为了进行比较，对回顾性欠采样的数据，进行grappa，espirit，和kipi 重建。其中grappa，espirit方法使用从第一欠采样图像帧获得的acs计算 权重或者espirit灵敏度图，重建所有图像帧。kipi方法使用两种灵敏度图， 一种是直接相除得到的原始线圈灵敏度图，另一种是espirit线圈灵敏度图。 本发明提出的rubustkipi使用上述提到的自动优化的线圈灵敏度图。
100.t1参数图像通过拟合得到，其中y是重建好的像素值，u和v 是模型适应变量。已知7组数据，求解包括t1在内的3个未知数。
101.通过比较grappa、espirit、kipi、rubustkipi与全采样真实值的归一 化均方根误差(nrmse)来评估rubustkipi方法的准确性。
102.3、结果分析
103.图1显示了rubustkipi组合原始线圈灵敏度幅值图和espirit线圈灵敏度 幅值图的流程图。第一支撑区域a为最终要形成的灵敏度范围，其应当完全覆 盖待成像的大脑区域范围且略微有一定的边缘外扩余量；第二支撑区域b为区 域a中排除了原始灵敏度图的噪声和空洞区域之后的区域。对比后可以看到， 第一支撑区域要大于第二支撑区域。而处于第一支撑区域但不在第二支撑区域内 的地方可视为噪声和空洞区域c。rubustkipi可以将准确的原始灵敏度图保留， 将噪声和空洞区域c内灵敏度图用espirit灵敏度图补全。
104.图2显示了通道7的三种不同的线圈灵敏度图的幅值和相位。分别为原始线 圈灵敏度图、espirit线圈灵敏度图、rubustkipi优化的线圈灵敏度图。可以看 出，原始线圈灵敏度图具有嘈杂的边缘和不平滑的跳变(箭头)，而espirit灵敏 度图非常平滑，但有点模糊。rubustkipi优化后的灵敏度图保留了细节，避免了 噪声和空洞的干扰，其相位图与原始的相位图有较高的一致性。
105.图3显示了从欠采样ir tse数据中得到的重建的水模源图像和对应的t1参数图像。再这种情况下，由于没有空洞区域且信噪比较高，水模本身比较均匀 kipi和rubustkipi都获得了高质量的结果，而grappa和espirit方法在源图 像和t1参数图像中都产生了明显的伪影(箭头)。
106.图4为grappa、espirit，kipi和rubustkipi重建的ti＝50ms的ir tse 源图像。kipi使用原始(raw)灵敏度图和espirit灵敏度图。grappa和espirit 的结果有明显的展开伪影，由第二列第三列的箭头表示，这在kipi的结果中是 没有的。使用espirit灵敏度图的kipi也表现出了一定的伪影，而rubustkipi 的误差最小(nrmse＝0.021)。
107.图5显示了从图4的源图像计算出的t1参数图像。grappa和espirit方 法重建的结果中可以看到明显的展开伪影；由于空洞和噪声区域，使用原始(raw) 灵敏度图的kipi方法表现出来轻微的伪影；而由于espirit灵敏度图本身缺失 一些细节的原因，使用espirit灵敏度图的kipi方法也表现出一定的伪影。而 rubustkipi结果中则没有伪影，灰质和白质中的误差非常小，总的误差也是最小 的(nrmse＝0.0686)。rubustkipi方法生成的t1参数图像最符合真实值。
108.以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。 有
关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以 做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案， 均落在本发明的保护范围内。