一种核磁共振t1噪声的压制方法与流程

【技术领域】

本发明涉及核磁共振数据处理技术领域，尤其涉及一种核磁共振t1噪声的压制方法。

背景技术：

核磁共振在生物大分子、分子动力学、药物研发等领域取得了重要应用。蛋白质等复杂结构的大分子使核磁共振谱图重叠严重，分辨率降低，对谱图的辨认造成困难。核磁共振多维谱是解决谱峰重叠的重要手段，但由于仪器的不稳定性使间接维的噪声水平增大，形成t1噪声，对谱图的辨认形成干扰。为提高核磁共振多维谱图的数据质量，如何压制t1噪声是核磁共振数据处理的重要课题。

现有的核磁共振t1噪声压制方法主要分为三种类型：1.在脉冲序列执行期间使用线性梯度场进行相干路径选择(hornetj,morrisga.p-typegradient-enhancedcosyexperimentsshowlowert1noisethann-type[j].magneticresonanceinchemistry,1997,35(10):680-686.)；2.通过抑制对角峰，以减小较强的对角峰形成的t1噪声(denkw,wagnerg,rancem,etal.combinedsuppressionofdiagonalpeaksandt1ridgesintwo-dimensionalnuclearoverhauserenhancementspectra[j].journalofmagneticresonance,1985,62(2):350-355.)；3.在数据处理阶段采用参考去卷积(gibbsa,morrisga,swansonag,etal.suppressionoft1noisein2dnmrspectroscopybyreferencedeconvolution[j].journalofmagneticresonance,seriesa,1993,101(3):351-356)、数据平滑(glasersj,kalbitzerhr.improvementoftwo-dimensionalnmrspectrabyweightedmeant1-ridgesubtractionandantidiagonalreduction[j].journalofmagneticresonance,1986,68(2):350-354.)、数据相关性去噪(simonpoulding,adrianj.charlton,jamesdonarski.removaloft1noisefrommetabolomic2d1h–13chsqcnmrspectrabycorrelatedtracedenoising[j].journalofmagneticresonance,189(2):190-199.)等。其中，基于数据处理的t1噪声压制方法应用较为广泛，但是，目前的方法普遍执行效率较低，且对谱峰的正负、宽窄等不同特征的普适性不强。

鉴于此，实有必要提供一种新型的核磁共振t1噪声的压制方法以克服上述缺陷。

技术实现要素：

本发明的目的是提供一种核磁共振t1噪声的压制方法，可以对噪声区域进行去噪，且压制方法的整个过程速度快，对各种谱峰类型均具有较好的适应性。

为了实现上述目的，本发明提供一种核磁共振t1噪声的压制方法，包括如下步骤：

s1：载入核磁共振多维谱数据s(m,n)，并计算直接维投影p1；其中m表示间接维的点数，即多维谱数据的矩阵s的行数；n表示直接维的点数，即多维谱数据的矩阵s的列数；

s2：将s1中得到的直接维投影p1采用连续小波变换求得一阶导数谱p2，并对一阶导数谱p2采用滑动窗口法得到谱峰区域，并将该区域作为t1噪声区域；

s3：取s2中非t1噪声区域的连续d列数据，分别求取各列的噪声水平，并将所有噪声水平的平均值σ1作为t1噪声压制的目标噪声水平；

s4：取s2中得到的t1噪声区域对应的矩阵s中的某列，采用连续小波变换得到一阶导数p3，并用滑动窗口算法得到该列的噪声区域p4；

s5：对s4中的t1噪声区域p4采用whittaker算法滤波进行平滑，使该区域的噪声水平被压制到目标值σ1；

s6：取s2得到的t1噪声区域对应的矩阵s的下一列数据，回到s4，直到t1噪声区域的平滑全部完成。

在一个优选实施方式中，s1具体包括，

步骤1.1：将多维谱数据s的负信号部分置零后的新矩阵，记为s⁺，然后将m行直接维数据逐条相加，得到正投影p⁺,

s⁺＝s

s⁺(s⁺＜0)＝0

p⁺(k)＝sum(s⁺(:,k))

其中，s⁺(:,k)表示矩阵s⁺的第k列，p⁺(k)表示正投影p⁺的第k个点；

步骤1.2：将多维谱数据s的正信号部分置零的新矩阵，记为s^-，然后将m行直接维数据逐条相加，得到负投影p^-,

s^-＝s

s^-(s^-＞0)＝0

p^-(k)＝sum(s^-(:,k))

其中，s^-(:,k)表示矩阵s^-的第k列，p^-(k)表示负投影p^-的第k个点；

步骤1.3：将正投影p⁺和负投影p^-合并，得到直接维投影p1，

p1＝p⁺+abs(p^-)

其中，abs(p^-)表示负投影的绝对值。

在一个优选实施方式中，s2具体包括，

步骤2.1：设置cwt参数r，

r＝a*n

其中，a为系数常量，通常取0.01～0.02，n为矩阵s的列数；

步骤2.2：使用cwt参数r，对间接维投影p1采用基于harr特征的连续小波变换求得一阶导数p2；

步骤2.3：将p2平均分为b段，分别求每段数据的标准差，取所有标准差的最小值作为p2的噪声水平ε；

步骤2.4：在p2上取总长度的1/c的数据段，计算该数据段最大值与最小值之差η，逐点滑动该数据段，当η＞f*ε时认为是峰区域，否则认为是t1噪声区域。

在一个优选实施方式中，s5具体包括，

步骤5.1：设置正则化参数λ，对噪声区域p4采用whittaker算法平滑：

minf＝min((y-z)²+λ(d²z)²)

z＝(1+λd'd)^-1y

其中，y表示噪声区域的原始数据，z表示平滑之后的数据；

步骤5.2：按照步骤2.3计算步骤5.1平滑后数据z的噪声水平σ2，若σ2>σ1,则增大平滑系数λ；

λ＝α*λ；

其中，α表示λ增大的倍数，通常的取值范围为：1<α<4；若σ2<＝σ1，则进入s6。

在一个优选实施方式中，所述连续d列数据的范围为4<＝d<＝20，b的范围为16<＝b<＝32，c的范围为20<＝c<＝100，f的范围为4<＝f<＝10，λ的范围为200<＝λ<＝500。

与现有技术相比，本发明提供的一种核磁共振t1噪声的压制方法，有益效果在于，对多维谱数据计算各维投影，并对各维投影采用连续小波变换和滑动窗口法得到谱峰区域，将该谱峰区域作为压制t1噪声的目标区域，在目标区域中各条数据再次采用连续小波变换和滑动窗口法得到各条数据的噪声区域，最后对噪声区域采用whittaker算法对噪声进行压制，压制方法的整个过程速度快，各种谱峰类型均具有较好的适应性。

【附图说明】

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中有关的附图作简要介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明提供的核磁共振t1噪声压制方法的流程图。

图2为本发明提供的核磁共振t1噪声压制方法的多维谱数据。

图3为本发明提供的核磁共振t1噪声压制方法的多维谱直接维投影和投影的一维导数谱。

图4为本发明提供的核磁共振t1噪声压制方法的多维谱的t1噪声区域。

图5为本发明提供的核磁共振t1噪声压制方法的多维谱(只包含正峰)t1噪声压制后的谱图。

图6为本发明提供的核磁共振t1噪声压制方法的多维谱(同时包含正负峰)t1噪声压制前的谱图。

图7为本发明提供的核磁共振t1噪声压制方法的多维谱(同时包含正负峰)t1噪声压制后的谱图。

【具体实施方式】

为了使本发明的目的、技术方案和有益技术效果更加清晰明白，以下结合附图和具体实施方式，对本发明进行进一步详细说明。应当理解的是，本说明书中描述的具体实施方式仅仅是为了解释本发明，并不是为了限定本发明。

请参阅图1，本发明提供一种核磁共振t1噪声的压制方法，包括如下步骤：

s1：载入核磁共振多维谱数据s(m,n)，并计算直接维投影p1；其中m表示间接维的点数，即多维谱数据的矩阵s的行数；n表示直接维的点数，即多维谱数据的矩阵s的列数；

s2：将s1中得到的直接维投影p1采用连续小波变换(continouswavelettransform，cwt)求得一阶导数谱p2，并对一阶导数谱p2采用滑动窗口法得到谱峰区域，并将该区域作为t1噪声区域；

s3：取s2中非t1噪声区域的连续d列数据(4<＝d<＝20)，分别求取各列的噪声水平，并将所有噪声水平的平均值σ1作为t1噪声压制的目标噪声水平；

s4：取s2中得到的t1噪声区域对应的矩阵s中的某列，采用连续小波变换得到一阶导数p3，并用滑动窗口算法得到该列的噪声区域p4；

s5：对s4中的t1噪声区域p4采用whittaker算法进行平滑，使该区域的噪声水平被压制到目标值σ1；

s6：取s2得到的t1噪声区域对应的s矩阵的下一列数据，回到s4，直到t1噪声区域的平滑全部完成。

具体的，s1：载入核磁共振多维谱数据s(m,n)，并计算直接维投影p1；其中m表示间接维的点数，即多维谱数据的矩阵s的行数；n表示直接维的点数，即多维谱数据的矩阵s的列数，具体包括如下步骤，

步骤1.1：将多维谱数据s的负信号部分置零后的新矩阵，记为s⁺，然后将m行直接维数据逐条相加，得到正投影p⁺,

s⁺＝s

s⁺(s⁺＜0)＝0

p⁺(k)＝sum(s⁺(:,k))

其中，s⁺(:,k)表示矩阵s⁺的第k列，p⁺(k)表示正投影p⁺的第k个点；

步骤1.2：将多维谱数据s的正信号部分置零的新矩阵，记为s^-，然后将m行直接维数据逐条相加，得到负投影p^-,

s^-＝s

s^-(s^-＞0)＝0

p^-(k)＝sum(s^-(:,k))

其中，s^-(:,k)表示矩阵s^-的第k列，p^-(k)表示负投影p^-的第k个点；

步骤1.3：将正投影p⁺和负投影p^-合并，得到直接维投影p1，

p1＝p⁺+abs(p^-)

其中，abs(p^-)表示负投影的绝对值。

具体的，s2：将s1中得到的直接维投影p1采用连续小波变换(continouswavelettransform，cwt)求得一阶导数谱p2，并对一阶导数谱p2采用滑动窗口法得到谱峰区域，并将该区域作为t1噪声区域；具体步骤如下，

步骤2.1：设置cwt(连续小波变换)参数r，

r＝a*n

其中，a为系数常量，通常取0.01～0.02，n为矩阵s的列数；

步骤2.2：使用cwt参数r，对间接维投影p1采用基于harr特征的连续小波变换求得一阶导数p2；

步骤2.3：将p2平均分为b(16<＝b<＝32)段，分别求每段数据的标准差，取所有标准差的最小值作为p2的噪声水平ε；

步骤2.4：在p2上取总长度的1/c(20<＝c<＝100)的数据段，计算该数据段最大值与最小值之差η。逐点滑动该数据段，当η＞f*ε时认为是峰区域(4<＝f<＝10)，否则认为是t1噪声区域；

具体的，s5：对步骤4中的噪声区域p4采用whittaker算法进行平滑，使该区域的噪声水平被压制到目标值σ1，具体步骤如下，

步骤5.1：设置正则化参数λ(200<＝λ<＝500)，对噪声区域p4采用whittaker算法平滑：

minf＝min((y-z)²+λ(d²z)²)

z＝(1+λd'd)^-1y

其中，y表示噪声区域的原始数据，z表示平滑之后的数据；

步骤5.2：按照步骤2.3计算步骤5.1平滑后数据z的噪声水平σ2，若σ2>σ1,则增大平滑系数λ；

λ＝α*λ

其中，α表示λ增大的倍数，通常的取值范围为：1<α<4。

若σ2<＝σ1,则进入s6。

实施例1：

一种核磁共振t1噪声压制方法，如图1所示，该方法包括以下步骤：

s1：载入核磁共振多维谱数据s(m,n)，m表示间接维的点数，即矩阵s的行数，n表示直接维的点数，即矩阵s的列数。并计算直接维投影p1。(在本实施例中：s为hsqc二维相敏谱，间接维点数m＝512，直接维点数n＝2048，多维谱数据s如图2所示)

步骤1.1、将多维谱数据s的负信号部分置零后的新矩阵，记为s⁺，然后将m行直接维数据逐条相加，得到正投影p⁺,

s⁺＝s

s⁺(s⁺＜0)＝0

p⁺(k)＝sum(s⁺(:,k))

其中，s⁺(:,k)表示矩阵s⁺的第k列，p⁺(k)表示正投影p⁺的第k个点；

步骤1.2、将多维谱数据s的正信号部分置零的新矩阵，记为s^-，然后将m行直接维数据逐条相加，得到负投影p^-,

s^-＝s

s^-(s^-＞0)＝0

p^-(k)＝sum(s^-(:,k))

其中，s^-(:,k)表示矩阵s^-的第k列，p^-(k)表示负投影p^-的第k个点；

步骤1.3、将正投影p⁺和负投影p^-合并，得到直接维投影p1。

p1＝p⁺+abs(p^-)

其中，abs(p^-)表示负投影的绝对值；

s2：将s1中得到的直接维投影p1采用连续小波变换(continouswavelettransform，cwt)求得一阶导数谱p2，并对p2采用滑动窗口法得到谱峰区域，并将该区域作为t1噪声的目标区域，(在本实施例中：p2及谱峰、噪声区域如图3所示)具体步骤如下，

步骤2.1、设置cwt参数r，

r＝a*n(在本实施例中：a＝0.01)

步骤2.2、使用cwt参数r，对间接维投影p1采用基于harr的连续小波变换求得一阶导数p2；

步骤2.3、将p2平均分为b段，分别求每段数据的标准差，取所有标准差的最小值作为p2的噪声水平ε(在本实施例中：b＝16，ε＝4.42*10⁵)；

步骤2.4、在p2上取长度1/c(在本实施例中：c＝100)的数据段，计算该数据段最大值与最小值之差η，逐点滑动该数据段，当η＞f*ε时认为是峰区域，否则认为是噪声区域；(在本实施例中：f＝6；直接维投影p1、一阶导数p2如图3所示，t1噪声区域如图4所示)

s3：取s2中非t1噪声区域的连续d列数据，按照步骤2.3分别求取各列的噪声水平，并将所有噪声水平的平均值σ1作为t1噪声压制的目标噪声水平。(本实施例中：d＝8，σ1＝1.945*10⁴)

s4：取步骤2的t1噪声区域对应的s矩阵的某列，按照步骤2中的连续小波变换得到一阶导数p3，并用滑动窗口法得到该列的噪声区域p4；

s5：对步骤4中的噪声区域p4采用whittaker算法进行平滑，使该区域的噪声水平被压制到目标值σ1，具体步骤如下：

步骤5.1、设置正则化参数λ(本实施例中：λ＝200)，对噪声区域p4采用whittaker算法滤波：

minf＝min((y-z)²+λ(d²z)²)

z＝(1+λd'd)^-1y

其中，y表示噪声区域的原始数据，z表示平滑之后的数据；

步骤5.2、按照步骤2.3计算步骤5.1平滑后数据z的噪声水平σ2，若σ2>σ1,则增大平滑系数λ

λ＝α*λ

其中，α表示λ增大的倍数(本实施例中：α＝2)。

若σ2<＝σ1,则进入s6；

s6：取s2得到的t1噪声区域对应的矩阵s的下一列数据，回到s4，直到t1噪声区域的平滑全部完成。(在本实施例中：t1噪声压制后的谱图如图5所示，本方法整个执行过程耗时2.8秒，测试电脑cpu为inteli57200u)

实施例2：

s1：载入核磁共振多维谱数据s(m,n)，m表示间接维的点数，即矩阵s的行数，n表示直接维的点数，即矩阵s的列数，并计算直接维投影p1。(在本实施例中：s为带谱编辑的hsqc二维相敏谱，谱峰同时包含正峰和负峰，间接维点数m＝512，直接维点数n＝2048，多维谱数据s如图6所示)。

步骤1.1、将多维谱数据s的负信号部分置零后的新矩阵，记为s⁺，然后将m行直接维数据逐条相加，得到正投影p⁺，

s⁺＝s

s⁺(s⁺＜0)＝0

p⁺(k)＝sum(s⁺(:,k))

其中，s⁺(:,k)表示矩阵s⁺的第k列，p⁺(k)表示正投影p⁺的第k个点；

步骤1.2、将多维谱数据s的正信号部分置零的新矩阵，记为s^-，然后将m行直接维数据逐条相加，得到负投影p^-，

s^-＝s

s^-(s^-＞0)＝0

p^-(k)＝sum(s^-(:,k))

其中，s^-(:,k)表示矩阵s^-的第k列，p^-(k)表示负投影p^-的第k个点；

步骤1.3、将正投影p⁺和负投影p^-合并，得到直接维投影p1，

p1＝p⁺+abs(p^-)

其中，abs(p^-)表示负投影的绝对值；

s2：将s1得到的直接维投影p1采用连续小波变换(continouswavelettransform，cwt)求得一阶导数谱p2，并对一阶导数谱p2采用滑动窗口法得到谱峰区域，并将该区域作为t1噪声的目标区域。

步骤2.1、设置cwt参数r，

r＝a*n(在本实施例中：a＝0.015)

步骤2.2、使用cwt参数r，对间接维投影p1采用基于harr的连续小波变换求得一阶导数p2；

步骤2.3、将p2平均分为b段，分别求每段数据的标准差，取所有标准差的最小值作为p2的噪声水平ε(在本实施例中：b＝20，ε＝6.75*10⁶)；

步骤2.4、在p2上取长度1/c(在本实施例中：c＝100)的数据段，计算该数据段最大值与最小值之差η，逐点滑动该数据段，当η＞f*ε时认为是峰区域，否则认为是噪声区域；(在本实施例中：f＝5)

s3：取步骤2中噪声区域的连续d列数据，按照步骤2.3分别求取各列的噪声水平，并将所有噪声水平的平均值σ1作为t1噪声压制的目标噪声水平；(本实施例中：d＝12，σ1＝2.465*10⁵)

s4：取步骤2得到的t1噪声区域对应的矩阵s的某列，按照步骤2中的连续小波变换得到一阶导数p3，并用滑动窗口法得到该列的噪声区域p4；

s5：对步骤4中的噪声区域p4采用whittaker算法进行平滑，使该区域的噪声水平被压制到目标值σ1，具体步骤如下：

步骤5.1、设置正则化参数λ(本实施例中：λ＝200)，对噪声区域p4采用whittaker算法滤波：

minf＝min((y-z)²+λ(d²z)²)

z＝(1+λd'd)^-1y

其中，y表示噪声区域的原始数据，z表示平滑之后的数据；

步骤5.2、按照步骤2.3计算步骤5.1平滑后数据z的噪声水平σ2，若σ2>σ1,则增大平滑系数λ

λ＝α*λ

其中，α表示λ增大的倍数(本实施例中：α＝2.5)。

若σ2<＝σ1,则进入s6；

s6：取步骤2得到的t1噪声区域对应的矩阵s的下一列数据，回到s4，直到t1噪声区域的平滑全部完成。(在本实施例中：t1噪声压制后的谱图如图7所示，本方法整个执行过程耗时2.5秒，测试电脑cpu为inteli57200u。)

本发明并不仅仅限于说明书和实施方式中所描述，因此对于熟悉领域的人员而言可容易地实现另外的优点和修改，故在不背离权利要求及等同范围所限定的一般概念的精神和范围的情况下，本发明并不限于特定的细节、代表性的设备和这里示出与描述的图示示例。