CN102854481A - 一种计算1h-mrs中物质浓度的方法 - Google Patents

Abstract

本发明公开一种准确计算¹H-MRS中各物质浓度的算法，适用于所有¹H-MRS的扫描序列。本方法包括小波阈值法去噪的环节，提高了信号的信噪比；采用提取波峰数据点并逐个对波峰进行曲线拟合的方法计算波峰的峰下面积；还包括采用内标准的方法实现峰下面积与物质浓度的转换。本方法可根据临床采集的实际数据对计算过程进行编辑处理，并可获得计算过程中的峰位、峰半高宽、峰下面积等参量，具有可编辑、准确高、重复性好等特点。通过本方法可对大脑内所有的代谢物的浓度变化进行准确计算及在体监测，并进行动态跟踪。这对研究大脑某些疾病的发病机制，如成瘾、癫痫、老年痴呆、肿瘤等，具有重要的临床价值与意义。

Description

一种计算1H-MRS中物质浓度的方法

技术领域

本发明涉及一种准确计算质子磁共振波谱（Proton Magnetic ResonanceSpectroscopy，¹H-MRS）中各物质浓度的方法。

背景技术

质子磁共振波谱（Proton Magnetic Resonance Spectroscopy，¹H-MRS）成像技术是一种在体、无创探测大脑代谢物的有效方法，并可定量测量这些代谢物的浓度。观察大脑代谢物浓度的变化，可为了解大脑的生理病理活动提供重要的信息，因此，¹H-MRS在临床中已有广泛的应用。健康大脑中，通常能被在体观察的代谢产物包括N-乙酰天门冬氨酸（NAA）、胆碱（Cho）、肌酸（tCr）、肌醇（myo-inosito，Ins）、脂类（Lipid）、谷氨酸复合物（Glu\Gln——Glx）、γ-氨基丁酸（GABA）等。

临床常用PRESS或STEAM序列在磁共振扫描仪上采集¹H-MRS数据。所扫描区域内的代谢物在不同的峰位形成一系列波峰,每个波峰代表一种代谢物,其峰下面积则表征该种物质的浓度。通常由于机器硬件等因素的影响,得到的是一系列含有噪声和有用信号的混合波峰信号。如何降低噪声、并准确地计算出¹H-MRS中各物质的浓度，是目前国际上¹H-MRS研究的重点之一。

LC-model(linear combination of model)是目前临床上最新的且应用最广泛的用来探测¹H-MRS中代谢物浓度的软件。该软件是由Stephen Provencher博士等主持编写的¹H-MRS后处理软件,具有高信噪比和窄的体外线宽，只要输入¹H-MRS数据，就能自动拟合出与原始数据接近的平滑的谱线形状和基线，计算出大脑中各物质的浓度。

但LC-model也存在一些问题:

1)LC-model是一个黑箱式的软件,其核心代码是封闭的,不是开源的,用户无法对其计算过程进行编辑处理。通常LC-model的处理方式都是由用户输入一组数据,确认运行后软件经过黑箱处理即给出最终结果。对于处理过程中一些重要的参量如峰位、峰半高宽、峰下面积等，软件都没有给出,而这些参量对于¹H-MRS的定量研究都是非常重要的。

2)对于¹H-MRS中一些相互重叠的波峰,如在3.00ppm处重叠的高浓度代谢物Cr与低浓度代谢物GABA等,LC-model则无法将重叠峰GABA分离提取出来,也就不能准确计算出该物质的浓度。

3)由于LC-model算法的限制,在LC-model标准库外的一些物质所代表的波峰,则无法对其进行计算处理。大脑中代谢物种类繁多，随着磁共振扫描仪软硬件的提高，在体1H-MRS扫描中一些新发现的物质所代表的波峰不断涌现，LC-model只能较好的处理一些常见物质所代表的波峰,但是对于一些新的波峰,LC-model则不能进行计算。

4)LC-model主要用来处理PRESS或STEAM序列所采集的数据，对于其它的脉冲序列,LC-model无法识别,需要开发配套的版本。随着1H-MRS研究的不断发展,目前开发了一些专用的探测某些特定代谢物浓度的脉冲序列,如用于检测GABA的mega-PRSS序列,LC-model无法对这些序列采集的数据进行处理。因此，需要开发一些普适的算法或软件，可对所有¹H-MRS的序列进行定量计算。

发明内容

本发明所要解决的技术问题是提供一种普适的计算¹H-MRS中物质浓度的方法。

本发明的思路是：首先采用小波阈值法对原始¹H-MRS数据降噪，然后采用数学拟合的方法准确地提取出某种代谢物的波谱形状；在此基础上，再准确计算波峰的峰下面积，最后将峰下面积转换成在体内的浓度。

本发明解决上述技术问题所采用的技术方案，包括以下步骤：

1）对临床采集的¹H-MRS数据采用小波阈值法去噪，提高信号的信噪比；

2）提取各物质的波峰数据，其中被大信号掩盖的小信号波峰和相互重叠的小信号波峰可采用本研究小组自主研发的新小波阈值法提取，亦可采用其他数据后处理的方法提取；

3）提取波峰的数据点，逐个对波峰进行曲线拟合，计算波峰的峰下面积；

4）采用内标准的方法实现峰下面积与物质浓度的转换。

步骤3）中，计算波峰的峰下面积时，对于单峰，逐个提取构成该峰的数据点，按照尽可能多的原始数据点落在曲线上的原则，选择Lorentz线型或Guass线型进行曲线拟合，再通过优化参数得到准确的峰下面积；

对于多重峰，当多重峰中每个单峰的数据点大于等于5时，将多重峰拆开为单峰分别进行拟合，然后叠加其峰下面积得到多重峰的峰下面积；当多重峰中的每个单峰的数据点小于5时，将多重峰拟合成一个大峰，拟合的原则是：中间最高峰的峰位保持不变，让尽可能多的左右小峰的数据点落在拟合后的曲线上。

下面进一步详细介绍本发明方法的计算步骤及原理如下：

1.准确提取临床¹H-MRS中的波峰

1)去噪：要获得准确的定量结果，首先必须对临床¹H-MRS数据进行去噪处理。不同于LC-model中算法，本文采用小波分析中的小波阈值去噪方法。

小波阈值去噪方法的基本原理如下:小波变换，特别是正交小波变换，具有很强的去数据相关性，能使信号的能量在小波域内集中在一些较大的小波系数中；而噪声能量却分布于整个小波域内。因此，经小波分解后，幅值较大的小波系数一般以信号为主，而幅值较小的系数在很大程度上是噪声信号。采用阈值的办法可将信号系数保留，而使大部分噪声系数减少至零。在此基础上，再重构信号,即可得到去噪后的“纯”信号。其计算公式见（1）式：

$c(\mathrm{scale},\mathrm{position})={\∫}_{-}^{+}f\left(t\right)\mathrm{\ψ}(\mathrm{scale},\mathrm{position},t)\mathrm{dt}---\left(1\right)$

其中f(t)是待处理的信号，即临床采集到的原始数据;ψ表示小波基函数；小波系数C是缩放因子（scale）和平移（position）的函数。通常噪声的小波系数较小。

在处理过程中可对具体的参数进行调节优化以得到最好的结果，如小波基函数db1~db8；sym1~sym8等可供选择,分解层数也可进行调节。

2)提取重叠峰：建议采用一种新的小波阈值法提取被大信号掩盖的小信号及重叠信号，如3.00ppm处被高浓度代谢物Cr掩盖的低浓度代谢物GABA的波峰；在2.30ppm处相互重叠的Glu\Gln的波峰；等。该方法见本研究小组申请的专利（申请号：201210080499.2）。也可采用其他的数据后处理方法提取小信号。

2.准确计算各波峰的峰下面积

提取出各物质的波峰后，下一步即为准确地计算出该波峰的峰下面积。以往MRS研究中对于提取MRS的重叠波峰研究较多，但对于如何准确地计算峰下面积往往重视不够，但这一过程对于定量计算物质浓度有极大的影响。本文提出了一整套可准确计算代谢物波峰的峰下面积的算法。

1)临床常见高浓度代谢物NAA、Cr、Cho等峰下面积的计算

①将表征某波峰的数据点截取下来。MRS信号由一系列波峰叠加而成,且临床数据受到采样时间的限制，是由一系列离散点组成的，并不是一条连续的曲线。将处理后的波谱图(频域)中某待处理波峰截取出来进行单独计算。这一过程可以有效地避免一张频谱图中含有多个波峰时存在的基线、相位、相互叠加等问题，而这些问题往往在很大程度上影响计算的准确程度。

②选取一种适当的函数来拟合这些数据点。¹H-MRS中物质波峰的信号理论模型一般为Lorentz信号（见公式（2）），即：

$y=y_0+\frac{2A}{\mathrm{\π}}\frac{w}{4{(x-x_c)}^2+w^2}---\left(2\right)$

其中x表示离散数据点横坐标的值，y表示离散数据点纵坐标的值；x_c,y_c为波峰信号最高点的坐标,w为半高宽,y₀为信号最低点的纵坐标,A为波峰的峰下面积（见图1）。令x=x_c,y=y_c,可得:

$A=\frac{\mathrm{\π}(\mathrm{yc}-y0)w}{2}---\left(3\right)$

A即为拟合后的Lorentz信号的峰下面积。可见，对于一个标准的Lorentz信号峰，只要知道yc、y0、w等参数，即可计算出曲线的峰下面积。

临床采集的数据受到各种硬件条件的限制,往往并不是标准的Lorentz信号。这就需要对待处理的波峰的数据点按照Lorentz函数进行拟合。拟合的原则是：使得拟合出来的Lorentz曲线有尽可能多的原始数据点落在上面。也有部分研究表明某些物质波峰可能更符合Guass信号模型，本文对部分数据同时进行了Guass拟合和Lorentz拟合，发现某些情况下Guass拟合优于Lorentz拟合。因此，在进行曲线拟合时，也可以对这两种函数按照拟合的原则进行优化选择。为了与一般软件所采用的Lorentz模型对照，下面按照Lorentz信号模型进行说明。

③通过优化参数得到准确的峰下面积

在拟合过程中一系列参数（主要为拟合后的Lorentz曲线的x_c,y_c,w）均为可调的。当参数选择不同时，拟合出来的曲线的峰下面积也有所差别，会影响最终计算出来的代谢物的浓度。经过对大量临床数据的计算分析发现，对于NAA、Cr、Cho等临床常见高浓度代谢物，若采用Lorentz拟合，则采取固定半高宽w和峰位x_c的方式得到的结果一般较稳定可靠。在采用这种方式处理某些数据得到的结果不够理想时（即不能满足“有尽可能多的原始数据点落在曲线上”的原则），再针对这些特异性数据选用适合的拟合参数。

2)多重峰的峰下面积计算：

1H-MRS还有一类波峰并非如NAA、Cr、Cho那样的单峰，由于J耦合作用而分裂成多重峰。这类波峰的峰下面积该如何准确计算，以往文献没有提到过。本文针对这种情况提出了特殊的处理方法。本文以3.00ppm处的GABA三峰举例说明（见图2）。

①根据¹H-MRS原理，应该将此GABA三峰拆开为3个单独的波峰分别进行计算，均参照1)的步骤拟合计算出各自的面积,然后叠加求得总的峰下面积（见公式（4））：

$A={\mathrm{\Σ}}_{i=1}^N{A}_i---\left(4\right)$

其中A_i为拆开后的某个单独峰的峰下面积,N为拆开后单峰的个数,A为叠加后的总面积。

但是临床数据采集由于磁场强度、噪声等因素的限制，当在¹H-MRS中某个峰的数据信号点过少(一般小于等于5)时，将无法进行准确的拟合，从而导致无法准确地计算出峰下面积。经过临床计算发现：对于一般临床采集到的¹H-MRS数据提取出的GABA信号峰，当将其拆开成3个单峰时，其左右2个单峰的数据信号点的个数往往小于等于5。

②参照某些扫描序列探测GABA信号峰只能得到一个大峰的方法(因分辨率不够。如mega-PRESS，在3.00ppm处只能采集到一个大峰)，当无法将获取到的GABA三峰分别拟合成3个单独的波峰时，可以将其拟合成一个大峰，再计算这个大峰的峰下面积。拟合的原则是：中间最高峰的峰位保持不变，让尽可能多的左右2个小峰的数据点落在拟合后的曲线上。

③参照1)中所列方法对拟合出的单个大峰进行峰下面积计算。

3)峰下面积与物质浓度的转换

计算出波峰的峰下面积后，还需将其转换为实际采集中的浓度。本文采用内标准的方法计算浓度，即：将某一¹H-MRS中含量较稳定的物质作为参考，将待测代谢物与这一参考值相比较，从而求出待测物质的浓度。内标准的数值是已知的，并且理论上应在人体各种生理和病理的情况下保持恒定。水在人体各种生理病理的情况下数值变化相对较小(约15%)，共振峰在4.7ppm。由于水信号的测量准确且可重复性强，可避免代谢物数值计算的系统误差，因此将水作为内标准。以水为内标准根据公式(5)即可将峰下面积转换为浓度。即：

其中n_水为水分子中¹H的个数，即n_水=2；A_水为水峰的峰下面积；N_水为水的标准浓度；n为待测物质分子中¹H的个数；A为待测物质的峰下面积；N为待测物质的浓度；p为修正因子，理论值为1，实际计算过程中的值可通过实验测量得到。

总之，本发明方法适用于所有¹H-MRS扫描序列，具有可编辑、准确高、重复性好等特点，可为在体研究代谢物浓度的变化，并进一步在体研究疾病的生理病理作用机制打下基础。本算法的优势在于单独对每个波峰进行提取，从而有效地避开了¹H-MRS后处理技术中复杂的基线问题，提高了计算的准确程度；另外，本发明方法采用开源设计，用户可根据临床采集的实际数据对计算过程进行编辑处理，并可获得计算过程中的峰位、峰半高宽、峰下面积等参量。

附图说明

图1、标准Lorentz信号示例。

图2、GABA理论峰型结果,在3.00ppm处为三峰状态。

图3、临床¹H-MRS数据,未去噪。

图4、图3数据经小波分析去噪后的结果

图5、将单个波峰(Cr)信号的离散数据点提取出来后的图示。

图6、将图5数据优化拟合成Lorentz函数并计算出峰下面积。

图7、本发明方法的流程示意图。

具体实施方式

（采用Lorentz拟合进行说明）

实施例1、临床常见高浓度代谢物，如NAA,Cr,Cho等的浓度计算：

在西门子3T磁共振仪上用PRESS序列采集10组人脑¹H-MRS数据。先采用小波阈值法对数据进行去噪处理，然后采用本文所述单峰峰下面积的计算方法逐个计算Cr的峰下面积，并转换成浓度。

下面将以一例具体的临床¹H-MRS数据举例来说明本算法的具体计算过程。图3为临床采集到的¹H-MRS波谱数据，图4为经过小波分析后得到的¹H-MRS波谱图。

将经过小波分析处理得到的¹H-MRS波谱图中的波峰Cr(3.00ppm)至起始点提取出来，见图5，横坐标表示峰位，纵坐标表示波谱峰值。

再将这些信号点拟合成Lorentz函数，按照“有尽可能多的原始数据点落在曲线上”的原则,调节拟合参数得到优化的结果，见图6。

根据算法中步骤3），将图6中拟合计算得到的峰下面积再转换为浓度。另，利用LC-model软件对相同的数据进行分析处理得到Cr的浓度，并将LC-model软件的计算结果与上述算法的结果进行比较（见表1）。由表1可见，两种方法具有一致性，且相对误差在可接受的范围内。

表1计算Cr浓度时本专利算法和LC-model软件结果对比

实施例2、多重峰，如GABA等，低浓度代谢物的峰下面积的计算：

在西门子3T磁共振仪上用PRESS序列和mega-PRESS序列分别采集7个人大脑同一部位（2种序列均采用相同体素）的数据，再分别计算2种序列得到的GABA峰下面积。PRESS序列采用本文计算方法；mega-PRESS序列是一种临床专门探测GABA的序列，用来与本文计算方法进行对比。表2为两种方法的计算结果。可见，两种计算结果具有高度的一致性。

表2计算GABA浓度时本专利算法与mega-PRESS结果的对比

Claims (3)

1.一种计算¹H-MRS中物质浓度的方法，其特征在于,包括以下步骤：

1）对临床采集的¹H-MRS数据采用小波阈值法去噪，提高信号的信噪比；

2 )提取各物质的波峰数据； 

3）采用提取的波峰数据点，逐个对波峰进行曲线拟合，计算波峰的峰下面积；

4）采用内标准的方法实现峰下面积与物质浓度的转换。

2.根据权利要求1所述的方法，其特征在于，

步骤3）中，计算波峰的峰下面积时，对于单峰，逐个提取构成该峰的数据点，按照尽可能多的原始数据点落在曲线上的原则，选择Lorentz线型或Guass线型进行曲线拟合，再通过优化参数得到准确的峰下面积；

对于多重峰，当多重峰中每个单峰的数据点大于等于5时，将多重峰拆开为单峰分别进行拟合，然后叠加其峰下面积得到多重峰的峰下面积；当多重峰中的每个单峰的数据点小于5时，将多重峰拟合成一个大峰，拟合的原则是：中间最高峰的峰位保持不变，让尽可能多的左右小峰的数据点落在拟合后的曲线上。

3.根据权利要求1或2所述的方法，其特征在于，步骤4）中，采用水作为内标物。

Images