基于傅里叶变换磁共振成像叠加式峰形的磁共振成像方法与流程

本发明涉及数学变换、信号处理、磁共振成像领域，特别是一种基于傅里叶变换磁共振成像叠加式峰形的磁共振成像方法。

背景技术：

数学变换被广泛地应用于信号处理，光谱分析，数字图像。所有这些技术都一直致力于如何提高信号的分辨率，降低噪音，快速实现数据获取与转换。现有计算机技术的飞速发展，对数据获取，存贮，运算和显示等基本不再是数学变换需考虑的制约因素。

磁共振成像是一项现代广泛应用的医学成像技术，原理基于当把会自旋的原子核放置在磁场中它们能吸收外加射频脉冲能量产生共振射频信号，给这个均匀主磁场内再施加三维空间梯度磁场，通过设定梯度与空间位置的数学关系，就能获取共振射频信号在空间的分布图。水分和脂肪是人体及生物器官组织中含有的主要成分，磁共振成像主要用于测量自旋氢原子核里的质子磁共振响应信号，即质子密度图像。

自旋原子核有一个非常重要的特征常数–磁旋比，磁场强度对应于原子核的共振进动频率。氢核磁旋比＝42.58兆赫兹/特斯拉，在1.5特斯拉磁场强度，氢核自旋进动频率为63.87兆赫；当磁场强度增大到3特斯拉，则氢核自旋进动频率增加到127.74兆赫。

磁共振成像步骤先在主磁场方向上设置身体部位截面磁场梯度，施加射频脉冲，选择需要扫描成像的截面，然后在与之正交的平面上做对应的竖坐标相位扫描和横坐标频率扫描，获取平面梯度编码的共振信号，专业上名称叫k-空间,这些编码按预定轨迹方向填充梯度磁场响应信号形成的一个原始数据矩阵。纵坐标梯度相位是有规律地随时间变化的，实际上也是频率扫描，因此采取对k-空间原始数据进行二维傅里叶变换转换为磁共振图像。

磁共振成像的空间分辨率主要取决于采样时间和信号噪音比。当主磁场强度增加一倍，简单地根据磁旋比可推出氢核自旋共振频率增加一倍，在同等测量时间条件下，磁共振频率测量误差由海森伯格测不准原理决定，因此相对频率误差就缩小一倍；按波尔兹曼能量分布理论，磁共振信号基本增强一倍，所以提高磁共振成像仪的磁场强度不但加快核自旋进动频率而且提升核共振信号强度，磁共振成像空间分辨率就被提高了。现有磁共振成像技术都是采用增大磁场强度来成倍提高成像的空间分辨率。

核磁共振成像是当代非常先进的医学检查手段，为了有效地研究病理，设计了各种各样的实验序列测定生物活性组织里的氢核共振效应，即测量水分子或脂肪的质子密度，还执行例如纵向弛豫成像和横向弛豫成像等等检测方法。磁共振成像做为一项高科技，需要调控很多技术参数来获取高质量的生物组织成像，但最终收集在k-空间里的原始数据都得通过傅里叶变换技术转换成图像。

1807年法国数学物理学家约瑟夫傅里叶首次提出并证明任何随时间周期变化的信号只要满足收敛条件能被展开为余弦和正弦三角频率函数的级数组合。当信号的周期趋于无限大时，则演变为著名的傅里叶变换。对于有限的n个这样的三角函数组合，积分形式的傅里叶变换可以用离散傅里叶变换来表达。

磁共振成像使用双通道正交检测器，磁场梯度共振响应是由实数通道和虚数通道组成的复数信号。一个频率为ω0磁共振梯度信号随时间t变化的函数可以表达成cos(ω0t)+isin(ω0t),其中它的傅里叶变换为

应用欧拉公式e^-iωt＝cos(ωt)+isin(ωt)，傅里叶变换又可表达为

实际信号都是在有限时间域进行分析的，所以傅里叶变换目前都是使用表一中列举的三种基本峰形进行的，即吸收峰形，发散峰形和幅度峰形，这三种傅里叶变换峰形都直接或间接影响磁共振像素之间的分辨。

表一

氢核在梯度磁场中的共振反应以傅里叶变换吸收峰形计算，可按如下公式计算单个像素面积(笛卡尔坐标xy平面的面积)，包括δx以及δy：

式中，γ为原子核的磁旋比，gx和gy为x和y方向的磁场梯度,而tx和ty则是x和y方向的全采样时间。

按现有理论傅里叶变换磁共振成像仪的空间分辨率与主磁场强度成正比。增大磁场强度一倍就能提高空间分辨率一倍，但成像仪的制造成本和技术要求都需要添加不止一倍的投入，而且仪器设备也将变得庞大笨重。考虑到电磁辐射对人体健康的影响，国际医疗机构规定临床使用的磁场强度不能超过3特斯拉。当然根据方程式1，磁共振像素也可以通过延长取样时间(或者加大梯度)缩小面积，但是像素变小到超过可容忍程度会使信号噪声比严重下降，空间分辨率反而降低。

技术实现要素：

本发明的目的在于提供一种基于傅里叶变换磁共振成像叠加式峰形的磁共振成像方法，以克服现有技术中存在的缺陷。

为实现上述目的，本发明的技术方案是：一种基于傅里叶变换磁共振成像叠加式峰形的磁共振成像方法，按照如下步骤实现：

步骤s1：通过一磁共振成像装置获取磁共振信号；

步骤s2：对所述步骤s1获取的磁共振信号进行傅里叶变换处理，获取傅里叶变换吸收峰形、傅里叶变换发散峰形以及傅里叶变换幅度峰形，并分别经离散化处理，获取离散化的吸收峰形、离散化的发散峰形以及离散化的幅度峰形；

步骤s3：通过一叠加函数，对所述步骤s2中所获取的峰形进行峰形叠加处理，分别获取傅里叶变换叠加吸收峰形傅里叶变换叠加发散峰形以及傅里叶变换叠加幅度峰形，并分别经离散化处理，获取离散化的叠加吸收峰形、离散化的叠加发散峰形以及离散化的叠加幅度峰形；

步骤s4：根据经峰形叠加处理后的信号，进行磁共振成像。

在本发明一实施例中，在所述步骤s3中，所述叠加函数为：

在本发明一实施例中，在所述步骤s1中，记磁共振梯度频率为ω0，强度为k，磁共振产生的k空间信号为：

f(t)＝2πk[cos(ω0t)+isin(ω0t)]0≤t≤t；

在所述步骤s2中，所述傅里叶变换吸收峰形为：

当有n个频率组合的k空间信号时，角频率级数ω＝2mπ/t和ω0＝2nπ/t；其中，m和n＝0，1，2，……，n-1，且所述离散化的吸收峰形为：

在所述步骤s3中，令x＝ω-ω0，通过所述叠加函数进行峰形叠加，得所述傅里叶变换叠加吸收峰形为：

所述离散化的叠加吸收峰形为：

在本发明一实施例中，，在所述步骤s1中，记磁共振梯度频率为ω0，强度为k，磁共振产生的k空间信号为：

f(t)＝2πk[cos(ω0t)+isin(ω0t)]0≤t≤t；

在所述步骤s2中，所述傅里叶变换发散峰形为：

且所述离散化的发散峰形为：

在所述步骤s3中，令x＝ω-ω0，通过所述叠加函数进行峰形叠加，得傅里叶变换叠加发散峰形为：

所述离散化的叠加发散峰形为：

在本发明一实施例中，，在所述步骤s1中，记磁共振梯度频率为ω0，强度为k，磁共振产生的k空间信号为：

f(t)＝2πk[cos(ω0t)+isin(ω0t)]0≤t≤t；

在所述步骤s2中，所述傅里叶变换幅度峰形为：

且所述离散化的幅度峰形为：

在所述步骤s3中，令x＝ω-ω0,通过所述叠加函数进行峰形叠加，得傅里叶变换叠加幅度峰形为：

所述离散化的叠加幅度峰形为：

在本发明一实施例中，在所述步骤s3中，所述峰形叠加处理还包括：根据所述磁共振信号的采样频率组分ω0进行频分组，对每个频分组内的所述傅里叶变换吸收峰形、傅里叶变换发散峰形以及傅里叶变换幅度峰形，采用所述叠加函数进行叠加运算。

在本发明一实施例中，在所述步骤s3中，所述峰形叠加处理还包括：对谐波组成的磁共振信号f(t)作离散化取样，记取样数为n，截取一组分立信号点f(0)，f(1)，f(2)，···，f(k)，···，f(n-1)；通过离散化傅里叶变换获取n个数据：f(0)，f(1)，f(2)，···f(k)，···，f(n-1)，通过一傅里叶变换矩阵表示如下：

其中，n×n傅里叶变换矩阵的w＝exp(-i2π/n)；

在所述傅里叶变换矩阵中增加一对角线叠加矩阵，获取经叠加运算后的傅里叶变换矩阵：

在本发明一实施例中，以采用逐列或在预置分辨本领条件δn区域内，决定所述对角线叠加矩阵元素是2，还是0,或者是接近0的小数。

在本发明一实施例中，基于磁共振成像以计算机二进制为基准设置像素数目，所述对角线叠加矩阵按如下矩阵，以2，0，2，0…排序方式，执行所述峰形叠加处理中左叠加运算；

所述对角线叠加矩阵按如下矩阵，以0，2，0，2…排序方式，执行所述峰形叠加处理中右叠加运算。

在本发明一实施例中，在所述步骤s3中，所述峰形叠加处理还包括：通过同步对相邻谐波信号进行前后左叠加运算或/和右叠加运算。

相较于现有技术，本发明具有以下有益效果：

(1)不同于使用复杂且成本极高的仪器构造来提高信号的辨识度，而是通过采用叠加处理达到百分之百提高信号的分辨率和强度。

(2)利用经典变换峰形(吸收峰形，发散峰形，幅度峰形)对称性采用叠加技术，实现缩窄一半现有傅里叶变换的峰宽度，增强信号强度一倍和提高分辨本领一倍。相当于用同等采样时间获取磁场增大一倍的效果，或者在同台仪器上，可缩短一半的采样时间，获取同等成像效果。

(3)对重迭相邻峰同步操作左叠加和右叠加，最佳可以提高4倍的分辨率。

(4)叠加函数可以优化应用，可采取全部频率组分都进行叠加。

(5)同样适用于基于k空间的数字成像技术，声波成像技术(超声波成像和声呐成像)，雷达成像等。

(6)本发明所提出的方法系国际首创，可以广泛地应用通讯，光谱学，数字成像等领域，具有适用性强，成本低，效率高和操作方便等优点。

附图说明

图1为本发明一实施例傅里叶变换吸收峰形示意图。

图2为本发明一实施例中傅里叶变换右叠加吸收峰示意图。

图3为本发明一实施例中傅里叶变换左叠加吸收峰示意图。

图4为本发明一实施例中傅里叶变换发散峰形示意图。

图5为本发明一实施例中傅里叶变换右叠加发散峰型示意图。

图6为本发明一实施例中傅里叶变换幅度峰形示意图。

图7为本发明一实施例中傅里叶变换右叠加幅度峰形示意图。

图8为本发明一实施例中磁共振成像仪中身体对应的区域示意图。

图9为本发明一实施例中傅里叶变换峰宽对图8中信号的影响。

图10为本发明一实施例中均等分布的10个磁共振频率做一组，采用现有傅里叶幅度变换与本实施例中提出的傅里叶变换做左叠加后的结果示意图。

图11为本发明一实施例中图10中信号经一傅里叶变换左叠加矩阵处理后的变换结果示意图。

图12为本发明一实施例中采用的人造膜磁共振原始k-空间图。

图13为本发明一实施例中采用的人造膜磁共振成像原图。

图14为本发明一实施例中对图12通过对角线叠加矩阵灰度增强后的对比图。

图15为本发明一实施例中同步对相邻谐波信号实行前后左(右)叠加和右(左)叠加处理后的示意图。

图16为本发明一实施例中对图12的原图经扩展一倍采样像素后用傅里叶叠加变换获取的清晰一倍的图象。

具体实施方式

下面结合附图，对本发明的技术方案进行具体说明。

本发明一种基于傅里叶变换磁共振成像叠加式峰形的磁共振成像方法，按照如下步骤实现：

步骤s1：记核磁共振成像装置的磁共振梯度频率为ω0，强度为k，磁共振产生的k空间信号为：

f(t)＝2πk[cos(ω0t)+isin(ω0t)]0≤t≤t；

步骤s2：对步骤s1获取的信号进行傅里叶变换处理，获取傅里叶变换吸收峰形、傅里叶变换发散峰形以及傅里叶变换幅度峰形，并分别经离散化处理，获取离散化的吸收峰形、离散化的发散峰形以及离散化的幅度峰形；

步骤s3：通过一叠加函数，对步骤s2中所获取的峰形进行峰形叠加处理，分别获取傅里叶变换叠加吸收峰形、离散化的叠加吸收峰形、傅里叶变换叠加发散峰形、离散化的叠加发散峰形、傅里叶变换叠加幅度峰形以及离散化的叠加幅度峰形；

步骤s4：根据傅里叶变换磁共振成像叠加式峰形处理后的信号进行磁共振成像。

进一步的，在本实施例中，对步骤s1获取的时间域信号按现有傅里叶变换处理会产生三种经典的峰形：

(1)傅里叶变换吸收峰形，如图1所示，其数学表达方程式就是常见的sinc函数：

当有n个频率组合的磁共振信号时，角频率级数ω＝2mπ/t和ω0＝2nπ/t，其中，m和n＝0，1，2，……，n-1，离散化的吸收峰形表达式为：

(2)傅里叶变换发散峰形，如图4所示，

离散化的发散峰形表达式为：

(3)傅里叶变换幅度峰形，也可以称为绝对值峰形，如图6所示，

离散化的幅度峰形表达式为：

进一步的，在本实施例中，傅里叶变换叠加处理采用如下所示的对叠加(superimpose)函数:

上述叠加函数能相当于执行把对称函数和反对称函数的一半叠加到函数本身的另一半，即右叠加或左叠加。

进一步的，在本实施例中，傅里叶变换里常见的符号函数定义为：

叠加函数跟符号函数有如下关系：

simp(x)＝1±sgn(x)(方程式7)

在实数域，其中带加号(+)的叠加函数正好是两倍于傅里叶变换里另一常见的阶跃函数h(x)：

simp1(x)＝2h(x)(方程式8)

阶跃函数定义为：

进一步的，令x＝ω-ω0，采用上述傅里叶变换叠加函数可衍生出三种新的基本峰形。

(1)傅里叶变换叠加吸收峰形，如图2和图3所示：

离散化的叠加吸收峰形表达式为：

(2)傅里叶变换叠加发散峰形：

如图5所示，右叠加发散峰型示意图。

离散化的叠加发散峰形表达式为：

(3)傅里叶变换叠加幅度峰形：

如图7所示，为右叠加幅度峰形示意图。

离散化的叠加幅度峰形表达式为：

进一步的，在本实施例中，本发明所提出的基于傅里叶变换的磁共振成像叠加式峰形处理方法，是对现有傅立叶变换基本理论和技术作了根本性的改革，采用数学手段叠加傅立叶变换的基本峰型，如吸收峰型，发散峰型，和幅度峰型(或功率峰型)等。利用这些峰型的对称性叠加增强信号的辨识度，达到百分之百提高分辨度，信号强度和信噪比。图1是角频率为ω0余弦信号进行傅里叶变换的吸收峰型，它是轴对称的，因此向右翻转叠加后，峰宽减少一半，峰高增大一倍，在翻转轴(原对称轴)左侧的干扰峰(吉布斯现象)降为零，见图2。当然，也可以如图3所示向左翻转叠加吸收峰，得到一样的效果。图4是对同样信号作傅立叶变换得到发散峰型。它是中心对称，可以通过旋转180°进行叠加。因为幅度峰型是吸收峰型和发散峰型各自平方加和后的平方差，所以图6中的幅度峰型也是轴对称的。本发明的叠加原理同样适合于发散峰型和幅度峰型。本发明所提出的傅里叶变换叠加原理改变了通常峰型观念，在保证信息无丢失的情况下追求不对称峰型表达方式，从而提升信号分析的质量。傅里叶变换最大的缺点是主峰会产生旁瓣谐波，如图1、图2以及图3中显示的峰形，也称为吉布斯现象，所以一般都采用截趾函数抑制它们。本发明所提出的方法能减少信号一侧的吉布斯现象，并通过截趾函数经矫正运算对另一侧信号进行处理。

进一步的，在本实施例中，对信号强度为k(k＝任意实数)，频率为ω0的磁共振梯度信号实施了时间为t的测量,表二和表三通过数值计算对比现有傅里叶变换理论与本发明提供的傅里叶变换理论的主要峰形技术参数。计算结果显示三种峰形都能被增高一倍的峰高，缩窄一半的峰宽，所以成倍提高了余弦信号的可辨识度。

表二

表三

进一步的，在本实施例中，方程式1说明磁共振成像的空间分辨率很大程度上依赖于磁场梯度和采样时间。受检身体部位是完全置于磁共振梯度线圈内，只要磁共振信号接受检测器有足够的灵敏度区分梯度的变化，那么分辨率就要取决于采样时间，见表二，因为傅里叶变换后的峰宽由采样时间决定。以图8中的四个紧邻小方块“11”，“12”，“21”和“22”代表磁共振成像仪中身体某个区域。在四个小方块内，氢核梯度场磁共振信号(省略自由感应衰减项)为：

“11”:k11[cos(ω11t)+isin(ω11t)]

“12”:k12[cos(ω12t)+isin(ω12t)]

“21”:k21[cos(ω21t)+isin(ω21t)]

“22”:k22[cos(ω22t)+isin(ω22t)]

其中，k11，k12，k21和k22为各自方块内信号强度；ω11，ω12，ω21和ω22为各自方块内梯度频率。当采样时间使傅里叶变换后的峰宽不足以区分开紧邻方块的梯度变化，见图9所示，用现有傅里叶变换技术就会把ω11和ω12的梯度频率的幅度峰形(细实线)在横向的扫描重合在一起(粗实线)。竖向扫描亦产生同样结果，所以空间分辨不佳时将导致这四个小方块(像素)里的梯度信号混成一个组合。核磁共振响应是带有自由感应衰减的时间信号，现有傅里叶变换方法的峰宽实际比表二列出的基本峰形峰宽还要宽得多，运用本发明的叠加处理方法能有效地缩窄一半各种峰形的峰宽。

在本实施例中，上述三种傅里叶变换基本峰形的叠加可通过如下几种方法实现，但不限于本实施例中所提供的方法。傅里叶变换磁共振成像通常采用幅度峰形，以下实施例中计算实现叠加式的傅里叶变换都以幅度叠加为例。

实施例一

现有磁共振成像是依照事先设置好的磁场梯度频率(或梯度相位)和采样像素做傅里叶变换，只要采样点数足够大，信号频率就不失真。而根据上述提供的傅里叶变换的幅度峰形中的方程，包括：方程式10，方程式11或方程式12，不需要预测额外参数，如图8所示，k-空间收集的原始数据即提供信号强度也有梯度频率数值。为了保证事无巨细全体测量范围内的梯度质子响应峰，可以对所有采样频率组分ω0做叠加运算，相当于有n个组分叠加就要耗时n倍时间来完成。可以进一步优化新型傅里叶变换运算，对频率组分进行适当分组。根据实际情况，采取两个频率一组、四个频率一组或八个频率一组等等，对各频率组根据实际要求达到的指标分别执行傅里叶变换叠加，可以大大有效地缩短运算时间，尤其现代计算机运算速度已大大提升，优化后的叠加运算不会增加太多。如图10所示，为均等分布的10个磁共振频率做一组，用现有傅里叶幅度变换得到10个幅度峰，以灰色细线显示，由于磁场梯度关系，在经过截趾函数处理后，峰与峰之间还有不同程度的重迭；当采用本实施例中提出的傅里叶变换做左叠加后，以粗实线表示，这10个峰就完全分开了。

实施例二

按现有傅里叶变换理论需对谐波组成的磁共振信号f(t)作离散化取样，设取样数为n，截取到一组分立信号点f(0),f(1),f(2),···,f(k),··,f(n-1)。离散化傅里叶变换则得出n个数据f(0),f(1),f(2),··f(k),···,f(n-1)，以矩阵方式表达

式中n×n傅里叶变换矩阵的w＝exp(-i2π/n)。

根据叠加理论在原傅里叶变换矩阵加一个对角线叠加矩阵，得：

进一步的，在本实施例中，可以采用逐列，或在预置分辨本领条件δn区域内，决定对角线矩阵元素是2还是0,或者是接近0的小数。由于磁共振成像以计算机二进制为基准设置像素数目，因此，矩阵对角可按如下方程式14.1，以2，0，2，0…排序，执行左叠加。

矩阵对角也可按方程式14.2，以0，2，0，2…排序，执行右叠加。

如图10所示的10个组分的磁共振梯度信号用上述傅里叶变换左叠加矩阵(方程式14,1)就得到如图11所示的变换结果，10个谱峰也被辨识和完全分开。

进一步的，为磁共振成像利用灰度直方图做成像矫正提供了理论基础。图12和图13所示的人造膜磁共振k-空间原始数据及其成像，获取自siemensverio3ttim核磁共振成像仪，配置3特斯拉超导磁场。主要仪器工作参数是自旋回波脉冲序列，矢状位图，层厚4.0毫米，驻留时间15.6微秒，脉冲重复时间600毫秒，回波时间6毫秒，像素带宽250赫兹，512x512像素。图14是通过对角线叠加矩阵灰度增强后的对比图，明显增强了成像反差。

进一步的，在本实施例中，由于叠加对角线矩阵列扫描跟矩阵行没有关联，故同样适用于快速傅里叶变换(fft)对叠加。但快速傅里叶变换执行的必须是方块矩阵。

实施例三

根据上述提出的两种对称叠加函数，见方程式5.1和方程式5.2，可同步对相邻谐波信号实行前后左(右)叠加和右(左)叠加，如图15所示，虚线所示局部重迭相邻的两个峰能够达到相当于提升四倍傅里叶变换的分辨本领，实线显示叠加变换的效果。

图13中的人造膜磁共振成像在扩展一倍的k-空间数据后用方法一进行傅里叶叠加变换，获得极为清晰如图16所示的1024x1024像素空间分辨率的图像。

以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。