一种基于信噪比加权的多b值扩散磁共振成像优化方法与流程

本发明涉及一种MRI成像技术应用中非线性最小二乘算法的优化改进方法，特别涉及一种基于信噪比加权的多b值扩散磁共振成像优化方法。
背景技术：
：基于多加权b值的扩散磁共振成像技术越来越成为扩散磁共振技术的发展趋势，包括扩散峰度成像、扩散谱成像等。同时，模型的阶数也在不断提高，例如从二阶的扩散张量成像到四阶的扩散峰度成像模型。因而模型的对于数据信噪比的敏感性也越来越敏感，因而在实际应用中，由于模型拟合的误差出现奇异值、拟合不稳定、无意义解的情况也逐渐增多，影响技术的稳定发挥和应用。本技术发明旨在通过一定的算法改进优化，提高多b值模型的寻优求解稳定性，进一步促进MRI成像技术的发展和应用。技术实现要素：针对现有技术存在的技术问题，本发明提供一种基于信噪比加权的多b值扩散磁共振成像优化方法，该方法是针对多b值扩散加权构建拟合模型，通过最小二乘算法的优化提高模型寻优求解的稳定性，进而使得MRI成像中参数图像的估计更加稳定，降低奇异值或无意义值的出现，进一步促进MRI成像技术的发展和应用。为了解决现有技术中存在技术问题，本发明采用如下步骤予以实施：一种基于信噪比加权的多b值扩散磁共振成像优化方法，包括下列步骤：步骤一，通过多b值加权数据建立拟合模型；步骤二，通过多b值加权数据估算信噪比；步骤三，根据信噪比计算残差修正权值αi，步骤四，建立拟合模型初始值；步骤五，根据残差修正权值αi计算拟合模型残差修正值；步骤六，如果拟合模型残差修正值符合收敛条件，则通过拟合模型的最优解计算提高其稳定性；否则，返回步骤四继续寻求拟合模型的最优解。所述步骤一中的拟合模型建立Si＝f(bi，β)，，其中β为待求解的未知参数，Si>0，i＝0,1,2,3…；b0＝0；bi>0，i＝1,2,3…。所述步骤二中信噪比的估算是包括全局信噪比和局部信噪比，即SNRi和SNR0。所述步骤三中残差修正权值αi采用如下公式：这里，b0＝0；bi>0,i＝1,2,3,….,bi常取值500,1000,1500,2000,2500,…。所述步骤五中拟合模型残差修正值计算采用如下公式：E(β^)=Σi=0nαi×(Si-f(bi,β^))2.]]>所述步骤六中拟合模型最优解计算采用如下公式：min.E(β^)=Σi=0nαi×(Si-f(bi,β^))2.]]>本发明有益效果：第一，本发明通过估计不同加权b值数据的信噪比加权，削减或者消除不同加权b值信噪比差异，进而减弱了噪声和拟合误差对模型拟合稳定性的影响。第二，如图2所示，本发明通过在模型优化求解过程中引入对残差平方和进行修正的权值向量，使得曲线拟合的参数结果更依赖于高信噪比的加权b值，即有限保证更高信噪比(较低)的b值点数据，从而减小高b值中噪声对曲线拟合的影响，使得参数的求解更加稳定。第三，如图3所示，奇异值的情况多出现在灰质和灰质区域，本发明将本来较高的灰质白质峰度值估计成较低值甚至零值。这里分别给出本发明方法对灰质(图4)、白质(图5)奇异值以及正常(图6)组织的效果。图4、图5和图6中的E和Ec，分别是修正前的求解函数曲面(式5)和修正后的求解函数曲面(式6)，图中的白色“+”是两种方法得到的解的位置；相应地，图中下方给出了不同方法得到的拟合曲线，其中蓝色均为使用本发明方法得到的曲线拟合结果。此三图结果明显展示了本发明方法的效果：能够很好解决拟合过程中出现的奇异值现象，同时不影响原本正确稳定的拟合结果。第四，如图7所示，给出了现有方法和本发明方法得到的扩散峰度成像参数图的结果，其中上方两图为现有方法得到的参数图，下方的两图为使用本发明方法得到的参数图。附图说明图1是本发明一种基于信噪比加权的多b值扩散磁共振成像优化方法流程图。图2是本发明一种基于信噪比加权的多b值扩散磁共振成像优化方法对残差以及拟合结果的影响图。图3是本发明一种基于信噪比加权的多b值扩散磁共振成像优化方法的模型拟合结果极易出现的奇异值示意(黑点)图。图4是本发明一种基于信噪比加权的多b值扩散磁共振成像优化方法的灰质或脑脊液的奇异值图。图5本发明一种基于信噪比加权的多b值扩散磁共振成像优化方法的白质部位奇异值图。图6是本发明一种基于信噪比加权的多b值扩散磁共振成像优化方法的优化方案对正常部位的图。图7是本发明一种基于信噪比加权的多b值扩散磁共振成像优化方法的优化方案的实际效果展示图(左：平均峰度；右：径向峰度)。具体实施方式下面结合附图对本发明做出详细地说明：如图1所示，本发明提供一种基于信噪比加权的多b值扩散磁共振成像优化方法，包括下列步骤：步骤一101，通过多b值加权数据建立拟合模型；所述步骤一中的拟合模型建立Si＝f(bi，β)，，其中β为待求解的未知参数，Si>0，bi>0，i＝0,1,2,3…。在实际中，多b值加权成像在目前的MRI技术应用中已经较为常见，其数据包括一个b＝0的无扩散加权的基准信号(图像)，和若干非零b值下的扩散加权成像信号(图像)，且这些非零b值扩散加权信号(图像)来自同一个对象。这里设基准信号为强度为S0＝S(b0)，而加权b值的信号为Si＝S(bi)(i＝1,2,3，…)。对于某个对象(仿体实物、动植物等所有非金属物)进行多b值扩散加权成像时，信号S的强度随着b值变大而不断减小，构成的曲线称之为衰减曲线，实际应用中用不同的模型进行拟合并定量描述对象的结构特性。这里将所有模型一般抽象地用Si＝f(bi，β)表示，其中β为待求解的未知参数，Si>0，bi>0，i＝0,1,2,3…。步骤二102，通过多b值加权数据估算信噪比；所述步骤二中信噪比的估算是包括全局信噪比和局部信噪比，即SNRi和SNR0。这里需要说明，在相同的机器环境中，不同b值数据(信号或者图像)中噪声的水平是相同，而由于不同b值加权得到的信号强度不同，使得不同b值加权数据的信噪比不同。这里需要首先估计不同b值加权图像的信噪比，包括全局信噪比和局部信噪比。其中，全局信噪比SNRi＝MeanSignal/MeanNoise。这里MeanSignal为第i个b值图像中信号(检测对象或目标区域的信号)的均值，MeanNoise为第i个b值图像中空气噪声信号(非检测对象或非目标区域的信号)的均值。局部信噪比SNRi通过局部空间(或者图像的局部区域/邻域)信号的均值除以信号的标准差来估计。步骤三103，根据信噪比计算残差修正权值αi，所述步骤三中残差修正权值αi采用如下公式：这里，b0＝0；bi＞0.i＝1，2，3，…，bi常取值500，1000，1500，2000，2500，··…（1）在数据曲线的拟合精度一定程度上依赖于采样点(已知点)数据的可靠性。因为不同加权b值下得到的信号其信噪比不同，因而其中有用信号或者数据的可信程度也不一样，这会对曲线或这模型的拟合造成一定的偏差。具体来说，采样点不等同的信噪比使得噪声对曲线拟合精度的影响更显著了，要想控制这一影响，首先需要根据信噪比估计不同b值下信号的可靠程度。步骤四104，建立拟合模型初始值；是对拟合模型中β的某个具体值(初始值)，采用如下公式计算：Si=f(bi,β^)+μ---(2)]]>其残差平方和为：E(β^)=Σi=0n(Si-f(bi,β^))2---(3)]]>步骤五105，根据残差修正权值αi计算拟合模型残差修正值；所述步骤五中拟合模型残差修正值计算采用如下公式：E(β^)=Σi=0nαi×(Si-f(bi,β^))2.---(4)]]>步骤六(106，107)，如果拟合模型残差修正值符合收敛条件，则通过拟合模型的最优解计算提高其稳定性；否则，返回步骤四继续寻求拟合模型的最优解。所述步骤六中拟合模型最优解计算采用如下公式：minE(β^)=Σi=0nαi×(Si-f(bi,β^))2.---(5)]]>上式子取极小值的一阶条件：dE(β^)dβ^=-2Σi=1nαi×(Si-f(bi,β^))×(-df(bi,β^)dβ^)=0---(6)]]>引入权值前后对于寻找最优解的过程和方法没有任何影响，只是残差平方和发生了改变，形象来说是最优解的位置发生了变化。改进的残差平方和通过对每个点等权值的累加，而改进后的残差平方和为不同点估计误差的加权平方和，因此在此新的残差和函数上寻求最优解。特别地，该残差和函数的维度与β的位置参数的个数一致。若β有两个未知参数，残差和函数的纬度为2，为一个二次曲面；而残差函数的极小值一阶条件也将变成两个偏一阶导联立等于零。本发明中拟合模型寻求最优解返回步骤四的过程是通过迭代寻优方法，其中以β＝(β0，β1)为例，取极小值的一阶条件为：dE(β0,β1)dβ0=-2Σi=1nαi×(Si-f(bi,β0,β1))×(-df(bi,β0,β1)dβ1)=0---(7)]]>dE(β0,β1)dβ1=-2Σi=1nαi×(Si-f(bi,β0,β1))×(-df(bi,β0,β1)dβ1)=0---(8)]]>扩散磁共振技术中多b值的扩散衰减模型多为非线性的，如二次曲线或者双指数曲线等，因而基本采用非线性最小二乘的求解思路。在具体实践中，也会加入一些对未知参数的约束条件，而采用一些不同的求解过程，大都也是非线性最小二乘的变式。非线性最小二乘法的思路是，通过泰勒级数将均值函数展开为线性模型。即，只包括一阶展开式，高阶展开式都归入误差项。然后再进行OLS回归，将得到的估计量作为新的展开点，在对线性部分进行估计。如此往复，直至收敛。这里给出高斯-牛顿(Gauss-Newton)迭代法求解非线性最小二乘问题的数值求解过程。对原始模型展开泰勒级数，取一阶近似值：f(bi,β^)≈f(bi,β^(0))+df(bi,β^)dβ^|β^(0)(β^-β^(0))---(9)]]>zi(β^)=df(bi,β^)dβ^---(10)]]>E(β^)=Σi=0nαi×(Si-f(bi,β^)-zi(β^(0))(β^-β^(0)))2=Σi=1nαi×(Si-f(bi,β^)+zi(β^(0))β^(0)-zi(β^(0))β^)2=Σi=1nαi×(St~(β^(0))-zi(β^(0))β^)2---(11)]]>构造并估计线性伪模型：St~(β^(0))=zi(β^(0))β^+ϵi---(12)]]>Et~(β^(1))=Σi=0nαi×(St~(β^(0))-zi(β^(0))β^(1))2---(13)]]>估计得到第一次迭代值并继续迭代。完成迭代过程如下：第一步：给出参数估计值的初始值将在处展开泰勒级数，取一阶近似值；第二步：计算和的样本观测值；第三步：采用普通最小二乘法估计模型得到β的估计值第四步：用代替第一步中的重复这些过程，直至收敛。收敛条件的设定有较大的自由度，例如设置为100次连续迭代残差等于0或者不变可认为收敛。本发明方法是力图解决扩散磁共振成像多b值模型拟合问题中，不同b值信噪比不同尤其是高b值信噪比很低对模型拟合结果精度和稳定性造成破坏的问题。以实际应用中扩散峰度成像(diffusionkurtosisimaging，DKI)模型为例，其拟合函数为：f(bi，β)＝-bi×D+(1/6)×(bi×D)2×K＝-bi×D+(1/6)×bi2×D×K这里β＝[DK]。因此该求解问题是一个2参数的函数(式5、6)，其问题为在一个二维曲面上求取最小值，如图3、4、5中的E和Ec。由于K是一个四阶统计量，对模型中的噪声和误差异常敏感，因而极易出现如图2所示的奇异值点。图4、5中，可以看出求解函数的解所在范围在K的维度上跨度更大，更体现了拟合误差在K值估计中的不稳定性。而本发明方法在DKI应用中的目的就是解决DKI对峰度估计的不稳定性问题。DKI使用应用中至少需要不少于2个的非零b值(b值量级为103，单位为s/mm2)，而且其中必须有不低于b＝2000的高b值，这里，高b值必将引入较低信噪比的拟合点，影响峰度的估计。采用本发明方法能够明显改善提高DKI模型的峰度估计，如图6所示。上述实例仅用于说明本发明，其中各部件的结构、材料、连接方式都是可以有所变化的，凡是在本发明技术基础上进行的等同变换和改进，均不应该排除在本发明的保护范围之外。当前第1页1 2 3