技术领域
本发明属于磁共振并行成像领域，具体涉及一种多支撑向量机模型的磁共振并行成像方 法。
背景技术
磁共振成像(Magnetic resonance imaging，MRI)由于其无核辐射，分辨率高，能多方位 和多参数成像等优点，已成为临床医学影像检查的重要手段之一。然而受到傅里叶编码方式 和奈奎斯特采样定理限制，磁共振成像速度较慢，这不但给患者带来一定的不适，而且容易 产生运动伪影。同时，长的扫描时间限制了MRI对运动物体的成像，如婴儿，血流，心脏 等。经过近几十年的发展，依靠提高硬件性能来加速采集的方式以达到人体承受极限。
并行成像技术使用多个线圈同时采集信号，利用各个线圈的空间敏感度来代替部分傅里 叶编码，从而达到减少扫描时间。磁共振并行成像方法有很多种，其中灵敏度编码和广义自 校准并行采集(Generalized Auto-calibrating Partially Parallel Acquisitions，GRAPPA) 是临床上较为常用的两种。GRAPPA算法是假设线圈敏感度的线性联合能产生空间谐波，达 到傅里叶编码的效果。首先，通过中间区域全采样的数据，作为校准数据，通过构造线性方 程组，求解线圈联合权重函数。再对其它区域未采集到的点进行联合权重填充。对于第l个 线圈坐标位置在(kx，ky+mΔky)的点，可以通过公式：
s l ( k x , k y + mΔ k y ) = Σ j = 1 L Σ h = - Hb Ha Σ b = - Nb Na w l , m ( j , h , b ) s j ( k x + h * Δ k x , k y + b * R * Δ k y ) ]]>进行确定。
其中kx，ky为K空间频率编码和相位编码方向上的坐标值，m为相位编码方向上的偏移 量，Δky为相位编码方向上单位间隔，s代表K空间信号值，Nb and Na为联合的相邻的相 位编码方向上所用的联合权重行数(分别上和下)，Ha和Hb为频率编码方向所用的联合 联合权重列数(分别为左和右)，j代表线圈编号，h为重建的频率编码方向所用联合权重点 序号，b为相位编码方向所用的联合权重点序号，w为联合权重函数，R为采样加速因子。 为了得到权重函数，GRAPPA算法在中间区域全采样，然后假设一些采集到的点为需要拟合 的点，组成一个知道输入和输出的线性系统方程：
AW＝B
通过最小二乘法求解得：W＝argmin||B-AW||2
实际上，由于联合权重函数的维数未知，只能通过人为的选择，因此GRAPPA算法具有 一定的随意性。如何选择最佳的一个相邻子集以及这个集合的大小一直是GRAPPA算法的难 点，一些方法通过迭代或者交叉验证的方法，去选择误差最小的相邻子集，但计算时间太 长，不能满足实时性要求。同时，通过最小二乘法求解，只能使校准区域内的匹配误差最 小，而没考虑联合权重函数的复杂度，容易出现过匹配，造成对未采集点的预测误差增大。
因此，如何提供一种具有最小全局误差，使得重建伪影较少的重建算法，同时计算速度 较快的重建算法已成为业界急需解决的问题。
发明内容
本发明的目的在于提供一种多支撑向量机模型的磁共振并行成像方法，其具有最小全局 误差，使得重建伪影较少、计算速度快。
为了实现上述目的，本发明的技术方案为：一种多支撑向量机模型的磁共振并行成像方 法，其特征在于，包括如下步骤：
(1)用多通道线圈对K空间中间区域进行全采样后，将其划分为训练集和检验集，其 它区域加速采样后作为预测集，并对各集合内的数据进行归一化处理；
(2)将训练集划分为多组训练子集，利用支撑向量机，选择不同的核函数和拟合参数 (C，v)对各训练子集进行训练，得到不同的联合权重函数模型；
(3)利用检验集，对不同的联合权重函数进行检验，选择最佳的几个模型，
(4)使用最佳的几个联合权重函数对预测集进行预测，取其平均值作为未采集点的 值，反归一化处理后，将K空间数据转换为图像。
进一步地，所述训练子集为Si＝{(x，y)i，(x，y)∈T，i＝1，...，Nt}，其是由训练集T＝{(x，y)} 等分或者随机的划分为Nt个子集，其中y为全采样中一些要拟合的值，x为所有线圈在坐标 上相应的邻域内采集到的值，将所有的y和对应的x列在一起组成训练集，
所述检验集为A＝{(x，y)，(x，y)∈T}，其是从全采样的数据中随机选取一部分，
所述预测集为P＝{(x，y0)}，其中y0为未采集到的点，设为0。
进一步地，在步骤(1)中，对数据进行归一化处理包括：将K空间数据映射到[-1，1] 上，其采取的公式为： x ′ = 2 * ( x - min ( x ) ) max ( x ) - min ( x ) - 1 , ]]> y ′ = 2 * ( y - min ( y ) ) max ( y ) - min ( y ) - 1 , ]]>
其中x为参考点的原始值，y为目标点的原始值，x′为归一化后参考点的值，y′为目标 点归一化后的值。
进一步地，步骤(2)中，选择一种改进的支撑向量机(v-SVR)，对不同的训练子集选择 不同的核函数和拟合参数(C，v)进行训练，解其约束优化问题：
min imizeτ ( w , ξ ( * ) , ϵ ) = 1 2 | | w | | 2 + C · ( vϵ + 1 l Σ i = 1 l ( ξ i + ξ i * ) ) , ]]>
subjectto((w·xi)+b)-yi≤ε+ξi
y i - ( ( w · x i ) + b ) ≤ ϵ + ξ i * ]]>
ξ i * ≥ 0 , ]]>ε≥0；
其中，设求解的线性函数为y＝wx+b，x为训练集中参考点的原始值，w为各参考点的 联合权重值，b为线性函数在Y轴上的截距，y为训练集中目标点的原始值，C是平衡因 子，调节目标函数的平滑度和误差，ε为误差精度控制项，可以被另一个变量v自动的调 节，v∈(0，1)，代表误差点所占比例的上限，l为训练样本总数，ξ和ξ*为松驰变量，通过 引入核函数并解此优化问题的对偶形式：
min α , α i 1 2 Σ i , j ( α i - α i * ) ( α j - α j * ) k ( x i , x j ) + Σ i y i ( α i - α i * ) ]]>
subject to Σ i ( α i - α i * ) = 0 ]]>
Σ i ( α i + α i * ) ≤ Cgv ]]>
αi， α i * ∈ [ 0 , C / l ] , ]]>i＝1，...l
其中，αi，为待求解的拉格朗日变量，k(xi，xj)为核函数。
进一步地，对于不同的训练子集，核函数分别选取线性核、径向基函数核和正则化的 傅里叶核，线性核的形式为：k(xi，x)＝xiHx，
径向基函数核形式为： k ( x i , x j ) = exp ( - | | x i - x j | | 2 / γ 2 ) , ]]>
傅里叶核形式为： k ( x i , x j ) = 1 - q 2 2 ( 1 - 2 q cos ( x i - x j ) + q 2 ) ; ]]>0＜q＜1，
其中γ和q分别为高斯核函数和傅里叶核函数的宽度参数，控制函数的径向作用范围。
进一步地，根据求解的拉格朗日变量αi，得到联合权重函数为 f ( x ) = Σ i = 1 l ( ai - a i * ) k ( x i , x ) + b . ]]>
进一步地，对于不同的训练子集，选取不同的平衡因子C和变量v，平衡因子C的范围 从0.005到500，其变化步长为10，变量v从0.1到0.9，其变化步长为0.1。
进一步地，在步骤(3)中将步骤(2)中训练得到的各个子集的权重函数应用到检验 集，对于每一个权重函数模型，其在检验集上的均方根误差定义为： RMSE = 1 n Σ i = 1 n ( y ^ - y ) 2 , ]]>
其中，n为检验集的大小，为通过权重函数计算得到的估计值，y是检验集里本身采 集到的值，最后选取均方根误差最小的Nm个子模型。
进一步地，在步骤(4)中利用步骤(3)选出的最佳模型分别对预测集内的点进行预 测，最后通过线性联合平均其中为每一个子模型预测结果，然后对数据进 行反归一化处理；并将K空间数据经过二维的快速离散傅里叶变换转换为图像，最后，通过 平方和的方法将各个线圈的图像联合成一幅最终图像。
进一步地，步骤(4)中的数据反归一化公式为： y = ( 1 + y ^ ) g ( max ( y ) - min ( y ) ) 2 + min ( y ) , ]]>且将数据从[-1，1]映射到原始空间。
本发明与现有技术相比具有如下优点：选择一种改进的支撑向量机(v-SVR)，因引入变量 v，能自动调节不敏感损失函数中参数ε的大小，使用起来比支撑向量回归算法得方便和准 确；利用支撑向量机求解得到的线圈联合权重函数，具有良好的泛化能力，对K空间噪声具 有良好的抑制能力，同时，可以选择非线性的训练核函数进行非线性的重建，使得当采样加 速因子过大时能较好去地除去混叠伪影，并用多个模型同时拟合，减少运算量，选择最佳的 几个模型，使重建误差和伪影更少。
附图说明
图1为并行成像图像重建算法的流程图；
图2为MRI加速扫描数据填充方式和GRAPPA重建的示意图。
具体实施方式
下面结合具体的实施例及附图对本发明进行详细描述
如图1所示，一种多支撑向量机模型的磁共振并行成像方法，包括以下步骤：
对于自校准类的并行成像算法，勿需另外单独扫描得到线圈敏感度，只需在扫描时进行 混合采样，对图2所示，其在K空间中，一部分相位编码按奈奎斯特采样速率进行采样，其 它部分则进行加速采样，对于R倍的加速，则每采集一条相位编码后，隔R-1条相位编码步 后再采集一条相位编码线。
对混合采样后的K空间数据划分为不同的集合，把全采样的数据作为训练集，其中一些 点作为要拟合的值即y，而其所有线圈在坐标上相对应的邻域内采集到的值作为x(如图 2)，将所有的y和对应的x列在一起组成训练集T＝{(x，y)}，把训练集等分或者随机的划分 为Nt个子集Si＝{(x，y)i，(x，y)∈T，i＝1，...，Nt}，同时从全采样的数据中随机选取一部分作为 检验集A＝{(x，y)，(x，y)∈T}，而把全部未采集到的数据作为预测集。
为了加快训练速度，需对训练集，检验集和预测集分别进行归一化处理，其归一化采取 公式：
x ′ = 2 * ( x - min ( x ) ) max ( x ) - min ( x ) - 1 ]]>
y ′ = 2 * ( y - min ( y ) ) max ( y ) - min ( y ) - 1 . ]]>
利用v-SVR，对不同的训练子集选择不同的核和参数进行训练，求解凸二次优化问题：
min α , α i 1 2 Σ i , j ( α i - α i * ) ( α j - α j * ) k ( x i , x j ) + Σ i y i ( α i - α i * ) ]]>
s . t . Σ i ( α i - α i * ) = 0 , Σ i ( α i + α i * ) ≤ Cv ; ]]>
αi， α i * ∈ [ 0 , C / l ] , ]]>i＝1，...，l
其中，αi，为拉格朗日变量，C是平衡因子，调节目标函数的平滑度和误差，ε为误差 精度控制项，可以被另一个变量v自动的调节，v∈(0，1)，代表误差点所占比例的上限，l为 训练样本总数。
核函数k(xi，xj)的选取不仅可以选择线性的核，还可以选取正则化的傅里叶核，以及径 向基函数核。在本实例中，对于不同的训练子集，分别选取线性核，径向基函数核和正则化 的傅里叶核，其形式为：
线性核：k(xi，x)＝xiHx；
径向基函数核：k(xi，x)＝exp(-||x-xi||2/γ2)；
傅里叶核： k ( x i , x j ) = 1 - q 2 2 ( 1 - 2 q cos ( x i - x j ) + q 2 ) ; ]]>0＜q＜1。
其中γ和q分别为高斯核函数和傅里叶核函数的宽度参数，控制函数的径向作用范围。 在优选的实施方式中，γ和q分别设置为0.05和0.3。
在优选的实施方式中，对于不同的训练子集，还有选取不同的参数C和v，C的范围从 0.005到500，其变化步长可以10，v从0.1到0.9。其变化步长可为0.1。
利用求解得到此的拉格朗日变量αi，可以得到联合权重函数为：
f ( x ) = Σ i = 1 l ( a i - a i * ) k ( x i , x ) + b . ]]>
通过上述核函数拟合得到的联合权重函数为非线性的，还可以选择如多项式，B样条等 其它非线性核。
对于不同的子模型可以得到不同的f(x)，然后将这些函数模型对检验集上的数据进行 预测，得到
将上述训练得到的各个联合权重函数模型应用到检验集，对于每一个权重函数模型，其 在检验集上的均方根误差定义为：n为检验集的大小，为通过权 重函数计算得到的估计值，y是检验集里本身采集到的值。最后选取均方根误差最小的Nm个子模型，经验显示，一般选择最佳子模型数达到10到15个时，其重建误差减小的趋势随 着子模型数的增加将变得缓慢。
利用选出的Nm个子模型分别对预测集内的点进行预测，最后通过线性联合平均 然后对数据进行反归一化处理，将数据从[-1，1]映射到原始的范围，其反归 一化公式为： y = ( 1 + y ^ ) g ( max ( y ) - min ( y ) ) 2 + min ( y ) . ]]>
并将K空间数据经过二维的快速离散傅里叶变换转换为图像。最后，通过平方和的方法将各 个线圈的图像联合成一幅最终图像。