基于外插磁化转移信号的CEST数据拟合方法、装置及介质

本发明属于磁共振成像技术领域，具体涉及一种基于外插磁化转移信号的cest数据拟合方法、装置及介质。

背景技术：

化学交换饱和转移(chemicalexchangesaturationtransfer,cest)成像是一种磁共振分子成像技术，可以检测活体里游离的蛋白质和多肽。cest成像已经被应用于多个临床领域，比如脑肿瘤分级、监测放化疗的早期治疗反应、预测肿瘤标志物、检测缺血性中风等。但是cest信号容易受多种效应的影响，国内外研究者提出了多种计算cest信号的方法以去除这些干扰效应，这些方法可以大致分成两类，一是使用参考信号与标记信号作差得到最终的信号，二是模型拟合得到最终的信号。非对称磁化转移(magnetizationtransferrationasymmetry,mtrasym)假设(1)参考信号中不包含cest效应；(2)标记信号包含cest效应；(3)其他效应在参考信号和标记信号中相等。然后用参考信号减去实验信号，就消除了其他效应只保留cest效应。但是其他干扰效应在参考信号和标记信号中并不都相等，比如大分子磁化转移(magnetizationtransfer,mt)效应，所以mtraysm方法得到的信号包含了其他效应。多池洛伦兹拟合假设每个效应的曲线是一条洛伦兹曲线，这些洛伦兹曲线叠加起来形成了实验采集的曲线，通过拟合实验采集到的曲线，可以得出各个效应的值，但是其他效应的曲线并不都是洛伦兹曲线，比如mt，它的曲线是超洛伦兹。外插磁化转移信号(extrapolatedsemisolidmagnetizationtransferreference,emr)拟合出的曲线代表除了cest效应之外的其他效应，将emr拟合的曲线与实验得到的曲线作差，就可以得到cest效应。emr假设实验使用的饱和脉冲是方波而且系统已经到达稳态，但是这两个条件并不总是满足的。

技术实现要素：

现有技术中的模型拟合方法对实验条件作了近似，但是这种近似在许多实验条件下并不成立。因此，本发明的目的在于解决现有技术中存在的问题，并在去除cest成像的干扰效应方面，提供一种鲁棒的、适合于多种实验条件的基于外插磁化转移信号的cest数值拟合(numericalfitofextrapolatedsemisolidmagnetizationtransferreferencesignal,简称nemr)方法。

本发明所采用的具体技术方案如下：

第一方面，本发明提供了一种基于外插磁化转移信号的cest数据拟合方法，该方法针对cest图像中的每一个体素，依次执行s1～s3，得到每一个体素的消除背景干扰信号的cest效应值；其中s1～s3为：

s1、对获得的原始z谱进行主磁场频率b0偏移校正，得到校正后z谱；

s2、取校正后z谱中仅含有mt效应的部分数据点作为拟合样本数据，对预先设置有待定参数初值和上下限的两池模型进行拟合，得到待定参数的拟合值；

所述两池模型表示为

式中：exp表示以自然常数e为底的指数函数；δt表示时刻t的增量；分别为t+δt时刻、t时刻的磁化矢量其中矩阵m＝[mxa,mya,mza,mzb]^t，mxa、mya、mza、mzb分别为水集团的x方向磁化强度分量、水集团的y方向磁化强度分量、水集团的z方向磁化强度分量、mt集团的z方向磁化强度分量；系数矩阵其中矩阵向量c＝[0,0,r1am0a,r1bm0b]^t，δωa表示外界施加的射频脉冲的频率和水的频率差，ω1表示所施加的射频脉冲强度，kab表示从水集团交换到mt集团的速率，kba表示从mt集团交换到水集团的速率，rrfb为mt集团对射频脉冲的吸收速率，δωb表示外界施加的射频脉冲的频率和mt的频率差，r1a＝1/t1a和r2a＝1/t2a分别为水集团的纵向弛豫速率和横向弛豫速率，r1b＝1/t1b和r2b＝1/t2b分别为mt集团的纵向弛豫速率和横向弛豫速率，m0a和m0b分别为水集团和mt集团的稳态磁化强度；

所述两池模型中的待定参数为t1a,t2a,t2b,m0b,kba；

s3、基于代入所述拟合值的两池模型得到包含mt效应的拟合曲线，将所述拟合曲线与所述校正后z谱上感兴趣频率处的差值，作为消去背景干扰信号的cest效应值。

作为第一方面的优选，所述s2中，两池模型中待定参数的拟合具有两轮，第一轮拟合时对5个待定参数在其预先设定的上下限范围内均利用所述拟合样本数据进行参数优化，获得第一轮拟合值；第二轮拟合时，以第一轮拟合值为中心，选择相对于第一轮拟合上下限范围更窄的上下限，重新使用和第一轮拟合中相同的拟合样本数据对所述两池模型进行第二轮精度更高的拟合，得到待定参数的最终拟合值。

再进一步的优选是，在进行第二轮拟合时，在所述两池模型的5个待定参数中仅选择部分参数调整其上下限，剩余参数的上下限与第一轮拟合保持一致；

更进一步的优选是，在进行第二轮拟合时，需改变上下限的参数组合为t1a,t2b,m0b或t1a,t2b或t2b,m0b。

作为第一方面中两轮拟合方案的优选，所述s2中，仅含有mt效应的数据点为频率位于80～6ppm范围内的数据点，优选为频率位于80～20ppm范围内的数据点。

作为第一方面中两轮拟合方案的优选，所述s2中，进行第一轮拟合时，待定参数t1a,t2a,t2b,m0b,kba预设的上限是[1.5，0.25，12e^-6，0.16，40]，预设的下限是[0.6，0.065，8e^-6，0.001，20]。

作为第一方面中两轮拟合方案的优选，所述s2中，在进行第二轮拟合时，对于需改变上下限的参数，其上限值改为该参数的第一轮拟合值的120％～150％，下限值改为该参数的第一轮拟合值的80％～50％；进一步的，上限值优选改为该参数的第一轮拟合值的120％，下限值改为该参数的第一轮拟合值的80％。

作为第一方面中两轮拟合方案的优选，所述s2中，mt集团对射频脉冲的吸收速率其中g(·)是一个线形函数，包括洛伦兹型、高斯型或超洛伦兹型。

进一步的优选是，所述s3中，所述感兴趣频率为3.5ppm或-3.5ppm，3.5ppm的感兴趣频率对应的cest效应值为apt^#值，-3.5ppm的感兴趣频率对应的cest效应值为noe^#。

第二方面，本发明提供了一种计算机可读存储介质，所述存储介质上存储有计算机程序，当所述计算机程序被处理器执行时，实现如第一方面中任一项所述的基于外插磁化转移信号的cest数据拟合方法。

第三方面，本发明提供了一种基于外插磁化转移信号数值拟合的cest数据拟合装置，其包括处理器和存储介质，所述存储介质上存储有计算机程序，当所述计算机程序被处理器执行时，实现如第一方面中任一项所述的基于外插磁化转移信号的cest数据拟合方法。

第四方面，本发明提供了一种去除cest成像干扰效应的磁共振成像装置，其特征在于，包括磁共振扫描器以及控制单元，所述磁共振扫描器用于通过磁共振cest成像获取cest图像；所述控制单元能获取所述cest图像且控制单元中存储有计算机程序，当所述计算机程序被执行时，用于实现如第一方面中任一项所述的基于外插磁化转移信号的cest数据拟合方法，输出每一个体素消除背景干扰信号的cest效应值。

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

本发明提出的拟合方法，本质上是对bm方程的分步求解，该方法有两大优点，首先，由于nemr在时间维度上进行迭代，可以将实际使用的射频脉冲的波形纳入到求解过程中，而不是假定波形为矩形波。其次，nemr方法采用bm方程仿真来描述cest效应，而不是经过简化的解析表达式，使得在系统非稳态时也适用，bm方程是描述cest效应最基本的方程，在此基础上进行仿真和迭代可以准确地还原生理参数值。

附图说明

图1为本发明的拟合方法流程图；

图2为实施例中实验采集的一条曲线和它的拟合曲线图；

图3为实施例中一位缺血性中风患者的t2加权图、弥散加权图、apt^#图和noe^#图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步阐述和说明。

本发明提供了一种去除cest成像干扰效应的方法，称为外插磁化转移信号的数值拟合(numericalfitofextrapolatedsemisolidmagnetizationtransferreferencesignal,nemr)。该基于外插磁化转移信号的cest数据拟合方法将mt效应的线型融入到bloch-mcconnell(bm)方程中，通过对该方程代表的模型的拟合，可以得出mt效应，将拟合得出的曲线与实验采集到的曲线作差，就可以去除干扰效应。

本发明的核心是通过对bm方程的合理近似从而将mt的线型融入到bm方程，形成一个用于进行后续拟合的两池模型。下面先对该两池模型的原理进行详细描述。

磁共振的成像过程可由如下的bloch方程来描述，

dmx/dt＝-mxr2-(ωw-ω)my

dmy/dt＝(ωw-ω)mx-myr2-ω1mz

dmz/dt＝ω1my-mzr1+m0r1

其中mx、my和mz分别为x、y和z三个方向的磁化强度分量，ωw和ω分别为水的磁共振频率和外界施加的射频脉冲的频率，ω1为所施加的射频脉冲的强度，r1＝1/t1和r2＝1/t2分别为纵向和横向弛豫速率，m0表示稳态磁化强度。该方程只描述了水的磁化强度在射频脉冲影响下随时间改变的情况。对于具有两个交换集团(水集团a；mt集团b)的系统，可由以下方程来描述(需注意的是下列公式中下标a代表水集团的相关参数，而下标b代表mt集团的相关参数)：

其中δωa＝ω-ωa表示外界施加的射频脉冲的频率和水的频率差；kij(i＝a或b；j＝a或b)则表示了i和j两个集团之间的交换速率(例如kab表示从水集团a交换到mt集团b的速率)。

由于mt集团的t2在微秒量级，mxb,myb会快速达到稳态，也即

0＝-δωbmyb-r2bmxb

0＝δωbmxb-r2bmyb-ω1mzb

可以解得如果定义为mt集团对射频脉冲的吸收速率，myb可以进一步表示为ω1myb＝-rrfbmzb，那么关于mzb的微分方程可以写为

上述mt集团对射频脉冲的吸收速率rrfb还可以用一个线型函数来表示，线形函数g(·)可以使用洛伦兹型、高斯型、超洛伦兹型等多种形式。其中：

洛伦兹型为

高斯型为：

超洛伦兹型为：

经过上述简化之后，原有的六个方程变成四个：

这四个微分方程可以表示为：

其中，m＝[mxa,mya,mza,mzb]^t，

c＝[0,0,r1am0a,r1bm0b]^t

进一步可以将其转化为齐次微分方程：

其中，

上述齐次微分方程的解可以表示为

其中为t+δt时刻的磁化矢量。

上述齐次微分方程的解即为本发明的两池模型，该模型有5个未知数，分别是t1a,t2a,t2b,m0b,kba，其余参数都是已知量。通过该模型的拟合，可以得出包含mt效应的拟合曲线，将此拟合曲线和经过主磁场频率b0偏移校正后的z谱曲线作差，就可以得到更为干净的cest效应。

需说明的是，本发明可以用于获取不同的cest效应值，但由于不同的cest效应值对应于曲线上的不同频率，因此为了便于表示，本发明中将z谱上能够反映某种cest效应的频率称为感兴趣频率。例如，假如需要计算的cest效应值为apt^#(amideprotontransfersignalfromnemranalysis，使用外插磁化转移信号分析得到的酰胺质子转移信号)值，那么其对应的感兴趣频率为3.5ppm；假如需要计算的cest效应值为noe^#(numericalfitofextrapolatedsemisolidmagnetizationtransferreferencesignal，使用外插磁化转移信号分析得到的核奥氏效应信号)，那么其对应的感兴趣频率为-3.5ppm。当然，假如z谱上具有其他频率也可以反映某种cest效应，那么也可以作为感兴趣频率，通过计算拟合曲线与校正后z谱上感兴趣频率处的差值，即可得到该感兴趣频率对应的消去背景干扰信号的cest效应值。本发明在实现过程中，可以单次计算一种cest效应值，亦可同时计算多种cest效应值，对此不作限定。

基于上述原理，本发明提供了两种基于外插磁化转移信号的cest数据拟合方法。如图1中的(a)和(b)所示，两种做法的区别在于对于两池模型的参数拟合方法不同，第一种方式为单轮拟合，第二种方式为两轮拟合，其余的做法均相同。下面对两种做法分别通过实施例进行展示说明。

实施例1

在本实施例中，如图1中(a)所示展示了采用上述第一种拟合方式的nemr拟合方法实现流程。该方法需要先获取cest图像，然后针对cest图像中的每一个体素，依次执行s1～s3，得到每一个体素的apt^#值和noe^#值，从而组成apt^#图和noe^#图。下面对于s1～s3的具体实现过程进行展开描述。

s1、对获得的原始z谱进行主磁场频率b0偏移校正，得到校正后z谱。主磁场频率b0偏移校正属于现有技术，对此不再赘述。

s2、取校正后z谱中仅含有mt效应的部分数据点作为拟合样本数据，对预先设置有待定参数初值和上下限的两池模型进行拟合，得到待定参数的拟合值。

如前所述，本发明中的两池模型表示为

式中：exp表示以自然常数e为底的指数函数；δt表示时刻t的增量；分别为t+δt时刻、t时刻的磁化矢量其中矩阵m＝[mxa,mya,mza,mzb]^t，mxa、mya、mza、mzb分别为水集团的x方向磁化强度分量、水集团的y方向磁化强度分量、水集团的z方向磁化强度分量、mt集团的z方向磁化强度分量；系数矩阵其中矩阵向量c＝[0,0,r1am0a,r1bm0b]^t，δωa表示外界施加的射频脉冲的频率和水的频率差，ω1表示所施加的射频脉冲强度，kab表示从水集团交换到mt集团的速率，kba表示从mt集团交换到水集团的速率，参数其中g(·)为超洛伦兹型，δωb表示外界施加的射频脉冲的频率和mt的频率差，r1a＝1/t1a和r2a＝1/t2a分别为水集团的纵向弛豫速率和横向弛豫速率，r1b＝1/t1b和r2b＝1/t2b分别为mt集团的纵向弛豫速率和横向弛豫速率，m0a和m0b分别为水集团和mt集团的稳态磁化强度。

上述两池模型中，需要对5个待定参数t1a,t2a,t2b,m0b,kba进行后续的参数拟合。仅含有mt效应的部分数据点可以根据实际的原始z谱进行分析确定，本实施例中优选为频率位于80～20ppm范围内的数据点。

对于第一种方式而言，可以直接利用拟合样本数据对5个待定参数t1a,t2a,t2b,m0b,kba进行拟合，在5个待定参数各自的上下限范围内，获取使整个两池模型方程拟合程度最高的5个待定参数的拟合值。具体的拟合方法可以利用matlab、spss等软件或者其他现有技术实现，拟合程度可以用最小均方误差来表示。

在进行拟合时，待定参数t1a,t2a,t2b,m0b,kba预设的上下限以及初值可以根据实际情况调整，只要符合各参数的正常取值范围即可。在本实施例中，数t1a,t2a,t2b,m0b,kba预设的上限设置为[1.5，0.25，12e^-6，0.16，40]，预设的下限设置为[0.6，0.065，8e^-6，0.001，20]，初值均设为各自的下限。

s3、基于代入5个拟合值的两池模型，即可得到包含mt效应的拟合曲线，双池模型中的频率ω即为曲线的横坐标，双池模型中的mza即为曲线的纵坐标。用该拟合曲线减去主磁场b0频率校正后的z谱，即可得到一条差值曲线，在这条差值曲线上取横坐标为3.5ppm和-3.5ppm对应的纵坐标值，分别记为apt^#值和noe^#值。当然，从简化角度来看，差值曲线不属于必须要获取的，只要能够计算拟合曲线上3.5ppm对应的纵坐标值与主磁场b0频率校正后的z谱上3.5ppm对应的纵坐标值间的差值，即可作为当前体素的apt^#值，同时计算拟合曲线上-3.5ppm对应的纵坐标值与主磁场b0频率校正后的z谱上-3.5ppm对应的纵坐标值间的差值，即可作为当前体素的noe^#值。

整张cest图像中每一个体素均获得各自对应的apt^#值和noe^#值后，即可形成一张apt^#图和一张noe^#图，apt^#图中的每个位置的值为该位置体素对应的apt^#值，noe^#图中的每个位置的值为该位置体素对应的noe^#值。

为了进一步展示上述方法的技术效果，将上述nemr拟合方法在一例缺血性中风病人的大脑cest图像上进行了测试。在本实施例中，采集图像所使用的cest成像序列包含三个模块：(1)cest饱和模块，该模块包含10个高斯波形饱和脉冲，每个脉冲持续时间约为100ms，幅度为1μt；(2)频谱预饱和反转恢复压脂模块；(3)快速自旋回波采集模块。

如图1(a)所示，首先对采集到的z谱数据(如图2中的实验曲线所示，频率范围80～-6ppm)进行主磁场频率b0偏移校正，然后进行上述nemr拟合。第一轮拟合使用频率为80～20ppm的数据作为拟合样本数据(拟合用到的点如图2中所示)，未知参数t1a,t2a,t2b,m0b,kba上限是[1.5，0.25，12e^-6,0.16,40],下限是[0.6，0.065，8e^-6,0.001,20]，初值等于下限。最终拟合曲线如图2所示。

本实施例的实验结果如附图3所示，可以看出，一步拟合的nemr方法拟合得出的apt^#图和noe^#图除了能够清楚地显示出其他模态上可以看到的病灶(t2加权图像上的病灶，白色实线箭头)，还能够显示出其他模态上不可见的病灶(noe^#图上的病灶，白色虚线箭头)。

实施例2

在本实施例中，如图1中(b)所示展示了采用上述第二种拟合方式的nemr方法实现流程。同样的，该方法需要先获取cest图像，然后针对cest图像中的每一个体素，依次执行s1～s3，得到每一个体素的apt^#值和noe^#值，从而组成apt^#图和noe^#图。

本实施例的s1和s3与实施例1做法相同，区别仅在于步骤s2中，两池模型中待定参数的拟合具有两轮：

第一轮拟合与实施例1相同，先对5个待定参数在其预先设定的上下限范围内均利用前述拟合样本数据进行参数优化，获得第一轮拟合值。5个待定参数的初始值也设置为各自的下限。

但是多参数拟合受变量上下限影响大，有些参数会等于给定的边界值，因此本实施例需要通过设置第二步拟合来减轻这一问题。第二轮拟合时，可以以第一轮拟合值为中心，选择相对于第一轮拟合上下限范围更窄的上下限，重新使用和第一轮拟合中相同的拟合样本数据对上述两池模型进行第二轮精度更高的拟合，得到待定参数的最终拟合值。在具体实现过程中，两池模型的5个待定参数不需要全部进行上下限调整，可以在两池模型的5个待定参数中选择部分参数进一步扩大其上下限值，剩余参数的上下限依然与第一轮拟合保持一致。本实施例在第二轮拟合时，需改变上下限值的参数组合为t1a,t2b,m0b。

本实施例中，通过将部分参数(如t1a,t2b,m0b)限制在第一轮结果的附近，可以使得这些参数在一个较小的范围内变动，从而减轻上下限的不当设置对最终结果的影响。

另外，对于需要改变其上下限值的第一轮拟合值，其上下限值的改变幅度不宜过大但也不宜过小。本实施例中，对于需要改变其上下限值的三个参数t1a,t2b,m0b，第二轮拟合所用的参数上限值优选改为该参数第一轮拟合值的120％，第二轮拟合所用的参数下限值优选改为该参数第一轮拟合值的80％，从而在第一轮拟合值附近对这几个参数进行进一步优化。

为了进一步展示上述方法的技术效果，将上述nemr拟合方法在一例缺血性中风病人的大脑cest图像上进行了测试。在本实施例中，采集图像所使用的cest成像序列与实施例1相同。

如图1(b)所示，首先对采集到的数据进行主磁场频率b0偏移校正，然后进行nemr拟合。第一轮拟合使用频率为80～20ppm的数据作为拟合样本数据，未知的待拟合参数t1a,t2a,t2b,m0b,kba上限是[1.5，0.25，12e-6,0.16,40]，下限是[0.6，0.065，8e-6,0.001,20]，初值等于下限。第二轮拟合时，使用和第一轮拟合中一样的拟合样本数据，将t1a,t2b,m0b的上限设为第一轮得到的t1a,t2b,m0b的120％，下限设为第一轮得到的t1a,t2b,m0b的80％，其余参数的上下限与第一轮拟合保持一致，所有参数的初值设定为该参数值的下限。

本实施例的实验结果如附图3所示，可以看出，两步拟合的nemr方法拟合得出的apt^#图和noe^#图除了能够清楚地显示出其他模态上可以看到的病灶(t2加权图像上的病灶，白色实线箭头)，还能够显示出其他模态上不可见的病灶(noe^#图上的病灶，白色虚线箭头)，且其显示效果比实施例1更优。

需要注意的是，上述两个实施例仅仅是本发明的两个实例，但并非仅限于此。例如，本发明的上述s2中的拟合样本数据的选择频率范围，可以选用其他范围，如80～20ppm，80～6ppm等。同样的，在进行参数拟合时，对于参数的上下限限制，可以根据实际情况调整，上限值可以改为该参数在第一轮得到的第一轮拟合值的120％～150％，下限值可以改为该参数的第一轮拟合值的80％～50％。再另外，其中第二轮拟合时，需改变上下限值的参数组合可以是多样的，例如除了t1a,t2b,m0b之外，还可以选择t1a,t2b或t2b,m0b等等。在进行拟合时，参数的初值设为上下限范围内的合理数值，也可以通过在上下限范围内随机或者均匀设置几个不同的待选值，所得拟合误差最小的待选值即可选定为最终的初值。

另外，在本发明的另一优选实施例中，上述基于外插磁化转移信号的cest数据拟合方法，可以以计算机程序的形式存储与计算机可读存储介质。当计算机程序被处理器调用并执行时，可以按照上述方法实现基于外插磁化转移信号的cest数据拟合方法。

同样的，在本发明的另一优选实施例中，还提供了一种基于外插磁化转移信号数值拟合的cest数据拟合装置，其包括处理器和存储介质，存储介质上存储有计算机程序，当计算机程序被处理器执行时，实现上述基于外插磁化转移信号的cest数据拟合方法。

计算机可读存储介质一般以存储器硬件形式提供，存储器可以包括随机存取存储器(randomaccessmemory，ram)，也可以包括非易失性存储器(non-volatilememory，nvm)，例如至少一个磁盘存储器。

上述处理程序的处理器可以是通用处理器，包括中央处理器(centralprocessingunit，cpu)、网络处理器(networkprocessor，np)等；还可以是数字信号处理器(digitalsignalprocessing，dsp)、专用集成电路(applicationspecificintegratedcircuit，asic)、现场可编程门阵列(field-programmablegatearray，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。

当然，随着云服务器的广泛应用，上述软件程序也可以搭载于云平台上，提供相应的服务，因此计算机可读存储介质并不限于本地硬件的形式。

在另一优选实施例中，可以将上述基于外插磁化转移信号的cest数据拟合方法，以程序的形式集成于磁共振成像装置的控制单元中。磁共振成像装置应当包括常规的磁共振扫描器以及控制单元，磁共振扫描器可采用现有技术实现，属于成熟商用产品，不再赘述。控制单元中除了存储有上述计算机程序之外，还应当具有实现cest成像所必要的成像序列以及其他软件程序。磁共振扫描器用于通过磁共振cest成像获取cest图像，而控制单元能获取磁共振扫描器得到的cest图像，且控制单元中存储有前述的计算机程序，当计算机程序被执行时，用于实现上述基于外插磁化转移信号的cest数据拟合方法，输出每一个体素消除背景干扰信号的cest效应值。

本领域的技术人员应当知道，本发明中所涉及的各模块、功能可以通过电路、其他硬件或者可执行的程序代码来完成，只要能够实现相应功能即可。若采用代码，则代码可存储于存储装置中，并有计算装置中的相应元件执行。本发明的实现便不限制于任何特定的硬件和软件结合。本发明中的各硬件型号均可采用市售产品，可根据实际用户需求进行选择。当然，磁共振cest成像序列及装置中，也需要配合必要的其他硬件或软件，此处不再赘述。

以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。