引言

目前，核磁共振技术广泛应用于生物医学、无损检测、石油勘探等领域．根据主磁场强度的不同，核磁共振技术可分为高场核磁共振技术和低场核磁共振技术．高场核磁共振技术一般应用于大分子结构分析，例如蛋白质的结构分析^[1]；在核磁共振成像方面，主磁场强度越高成像质量越好，因而核磁共振成像通常也是高场核磁共振技术的主要应用之一．与高场相比，低场核磁共振技术主要应用在能源勘探、地下水源寻找、灾害防治、食品化工的在线无损检测等方面．同时低场核磁共振在使用弛豫信息进行分析和检测的应用领域中具有无可替代的地位．首先，为了准确得到被测样品的弛豫谱，必须给予足够多的数据，而高场强情况下单次扫描的脉冲个数受到限制^[2]只能提供有限的数据．如果使用低场核磁共振检测技术，可以很容易的获取多个回波数据．其次，当样品中存在顺磁性物质或铁磁性物质时，高场强下诱发的强内部梯度场会严重影响到主磁场的均匀性，因此只能使用低场设备进行检测．最后，低场核磁设备易于小型化，更容易进行磁屏蔽和电磁屏蔽^[3]，而且价格低廉，这些特点推动了低场核磁共振检测技术在各行各业中的广泛应用．

CPMG(Carr Purcell Meiboom Gill,CPMG)序列是低场检测技术中最常用的序列之一，通过对CPMG回波数据进行反演可以得到与被测样品的结构、性质等相关的T₂谱^[4]．但是由于主磁场强度低，接收到的信号微弱^[5]，真实的回波信号容易淹没在背景噪声中，严重影响到后续反演等操作的精度和准确性，这一问题在高温、高压等极端环境下进行的磁共振测井实验中尤为明显．因此，有效抑制噪声以提高信噪比，成为反演工作最重要的预处理环节．在众多算法中，Buades等人^[6]于2005年提出的非局部均值(non-local means,NLM)滤波方法取得了良好的滤波效果，已成为图像去噪领域的一个新的研究热点．该算法利用自然图像自身大量的冗余信息，用非局部自相似性来抑制图像噪声．自相似性是信号的一个基本性质^[7]，一维信号也具有自相似性，因此将NLM算法应用于一维信号去噪具有可行性．文献^{[8, 9, 10]}将其应用于心电信号的去噪并取得了很好的效果．文献^[11]利用NL M 对Blocks等标准一维测试信号进行去噪，在信噪比和视觉效果上都高于包括小波去噪在内的传统算法．

NLM算法中搜索窗宽度、相似窗宽度和衰减参数等3个参数对去噪效果起到决定性作用^[12]．衰减参数是一个控制指数衰减程度的参数，通过影响权重的大小而决定算法的滤波程度．如果取值太小，则滤除噪声不彻底；如果取值太大，则会导致滤波信号过于平滑，即产生过滤波现象．理论上来说，搜索窗宽度应该尽可能的大，但是增大搜索窗宽度的同时计算量也会增加．相似窗宽度决定了在多大的尺度上表示信号的结构，如果选的太小就不足以体现每个点的结构特征，对噪声的鲁棒性也不好；如果选的过大，就不容易找到相似的邻域结构，限制了算法的去噪能力．实际应用中，通常需要依据经验或者反复尝试来选择合适的参数，有时浪费了大量的时间也很难获得较好的去噪效果．针对这一问题，本文提出了利用Stein无偏风险估计(Stein’s unbiased risk estimated,SURE)的自适应参数选取算法．

另外，CPMG回波信号具有特殊性，信号前端变化较快，因而相似窗宽度不宜太大，否则无法体现回波信号间的相似性，降低去噪效果；而尾端平坦部分几乎全是噪声，应采用较大的相似窗宽度来去除噪声以达到平滑的作用．由于CPMG回波信号的特性，采用固定的相似窗宽度并不能达到最优的去噪效果，因此本文在自适应参数选取算法的基础上提出了自适应相似窗宽度的NLM算法．

1 CPMG回波信号分析

CPMG是低场NMR分析技术中最常用的序列，该序列由一个90°脉冲和多个180°脉冲组成，如图 1所示．通过施加90°脉冲，使磁化矢量反转到XOY平面，经过一定时间磁性核开始散相，在未达到完全散相时，施加180°脉冲使磁性核的运动方向相反，散相过程变为聚相过程，会出现一个峰值，得到一个回波．后续的180°脉冲都会产生相同的效果，但是随着能量的损失，回波波峰的幅度会逐渐减小，最终达到0（如图 2所示）．由各个回波波峰组成的指数衰减曲线就是T₂衰减曲线（图 2中的虚线），可以利用此衰减曲线进行反演得到T₂谱，但是采集到的CPMG回波信号都是含有噪声的，对后续的反演工作会造成很大的影响．同时不同回波时刻的数据对反演结果的影响程度不一，前端数据保存了绝大部分信息对反演谱较重要，因此要避免前端信号出现损失，后面信号量几乎衰减为0的部分对反演结果影响较小．由于CPMG信号的这种特性，不宜套用普通的算法进行去噪，而应根据每一部分的特性进行滤波算法设计．

CPMG回波波峰幅度变化相对平滑，连续的信号段具有相近的幅度，尾端信号变化趋势平缓，相似性更加明显，而NLM算法正是利用信号间的相似性进行滤波．同时CPMG信号的尾端几乎全是噪声，与前段信号的相似性很小，所以以信号幅度为相似性判断依据的NLM去噪中对前端的影响也较小．综上所述，NLM算法适合用于CPMG信号的去噪．

2 NLM滤波算法

NLM处理的是包含噪声的信号，假设该模型为v = u + n，其中n为噪声信号，u为真实的信号，v为实际采集到的信号数据．对于任意位置s处的信号值，其去噪后的值$\hat u(s)$是在搜索范围内其它点处幅值的加权平均和．公式如下：

其中$z(s) = \sum\limits_{t \in N(s)} {\omega (s,t)}$，是归一化因子，权重$\omega (s,t) = \exp \left\{ { - \frac{{\sum\limits_{\delta \in B} {{{[v(s + \delta ) - v(t + \delta )]}^2}} }}{{2L{\lambda ^2}}}} \right\}$，N(s)是以s点为中心的搜索区域，B是以s点为中心的相似邻域，L是邻域B中所包含点的个数，例如B = [-10,10]，则L为21．

NLM算法的核心思想是在一定的搜索范围内搜寻与被滤波点相似或匹配的其它点参与到滤波过程中，以实现更好的滤波效果．该算法的基本操作为在以被滤波点s为中心的搜索区域N(s)内，计算该区域内每一个点t与点s的相似性来决定点t的权重大小，最后将搜索区域内每个点的幅值与该点相应权重乘积的和作为点s去噪后的幅值．搜索区域内的两个点s和t的相似性，则通过中心分别位于这两个点、称为相似窗的两个线段N_s和N_t中所有点对应位置的幅值差的平方和来计算．该方法有3个关键参数：相似窗半径R_sim、搜索窗半径R_sear、衰减参数l．这3个参数的选择直接影响到NLM的滤波效果．示意图如下：

3 基于SURE的参数选取策略

在已知真实信号u的情况下，均方误差(Mean Squared Error,MSE)是评价去噪效果最直接和有效的标准，MSE的计算公式如下：

(2)式中，I表示信号中的所有点的集合．对于去噪后的信号，MSE的值越小就表示信号的去噪效果越好，因此可以通过MSE判断NLM算法中选取怎样的参数可以达到最优的去噪效果．然而，在实际应用中采集到的都是含噪信号，所以均方误差并不能真正用于参数的选择．

SURE算法是对MSE的无偏风险估计，可以在真实信号未知的情况下对MSE做出准确的估计．SURE算法由下面的解析表达式表示^[13]：

其中${\text{di}}{{\text{v}}_v}\{ \hat u\}$是NLM算法的微分，${\text{di}}{{\text{v}}_v}\{ \hat u\} = \sum\limits_{s \in I} {\frac{{\partial {{\hat u}_s}}}{{\partial {v_s}}}}$，SURE算法假定噪声为高斯白噪声，并假定已知噪声方差${\sigma ^2}$．在实际当中，很容易从测得的数据中估算出噪声方差．微分表达式在整个SURE计算公式中具有重要的作用．对于NLM算法，能够从(1)式中得到这个微分表达式，计算公式如下（详细推导过程见参考文献^[12]）： 其中${W_s} = \sum\limits_{t \in N(s)} {\omega (t,s)}$，是对含噪信号v的平方进行NLM去噪的结果，公式为：$u_s^2 = \frac{1}{{Z(s)}}\sum\limits_{t \in N(s)} {\omega (s,t)v_t^2}$，式中的权重因子与(1)式中的权重大小一致．将(4)式代入(3)式得到SURE的最终表达式：

微分公式的计算可以和NLM算法同时进行，此外运用(5)式需要两个额外的向量来存储$u_s^2$和W_s，其计算复杂度仅仅为$O(L \cdot N)$，和NLM算法的复杂度$O(L \cdot N \cdot S)$相比是可以忽略的（S为搜索区间点的个数）．

下面用实验验证SURE算法对MSE的估计效果．首先使用中心位于T₂ = 10 ms处的高斯峰进行正演，然后向其中添加一定方差的高斯白噪声，最终得到含噪声的模拟信号．在计算SURE公式的时候使用了两个噪声方差，即真实的噪声方差和估计的噪声方差．CPMG信号末端信号量理论值为0，其末端数据可以认为是噪声，可用于估计噪声方差．图 4显示了在对模拟信号添加不同方差的噪声后，MSE和两组SURE值（分别由真实噪声方差和估计噪声方差得到的值）随l的变化曲线．图 5显示了在噪声方差固定，只改变相似窗宽度的情况下，MSE和两组SURE值随λ的变化曲线．

从图 4中可以看出NLM算法的去噪效果依赖于λ，MSE的变化趋势可以用SURE来替代．最优的λ与噪声标准差有关，也就是最优λ/σ的值随着噪声标准差的改变而变化，同时从图中可以粗略的看出最优λ在噪声标准差σ值的附近，因此可以利用噪声标准差粗略的估计出λ的范围，以减少迭代次数．在实际当中可以不用进行这么多次的迭代，可以用较少的迭代然后根据变化趋势利用黄金分割搜索的方法来确定变化曲线的最低点．

从图 5和表 1中可以看出，当相似窗宽度变化时，最优λ几乎不变．但随着相似窗宽度值的增加，信噪比和计算时间也增加，所以选择参数时要有所折中．综上所述，SURE公式可以代替MSE用于参数的选择．

4 回波信号去噪优化

传统NLM算法在进行滤波时总是选用固定的相似窗宽度，但这种方法并不能对所有的信号都取得良好的去噪效果．问题在于如果相似窗宽度取的较大，那么将会使信号中的一些细节部分被滤除掉，还会增加两个点的相似性权重以致某些特殊的点（不含相似性区域的点，如极大值或极小值点）被平滑；如果相似窗宽度较小，那么在信号变化平坦的区域就不足以体现每个点的结构特征，对噪声的鲁棒性也不好．也就是说在信号缓慢变化的区域（方差较小）应该给予较大的相似窗宽度，而在包含细节的区域（方差较大）需要较小的相似窗宽度才能取得良好的去噪效果．对于CPMG信号而言，前端数据保存了绝大部分信息，因此去噪时要避免前端信号出现损失，所以相似窗半径不宜过大；尾端部分信号量趋于平坦，应采用较大的相似窗半径来去噪以达到平滑的作用．

由于CPMG信号的上述特性，采用固定的窗宽半径，去噪效果并不是十分理想，为了达到更好的去噪效果，应该根据区域的特点自适应的改变相似窗宽的大小．

对于相同的窗宽区域，CPMG信号前端数据的数据方差较大，后端信号数据方差较小，因此可以根据数据方差的不同来自适应的改变不同去噪点的相似窗宽度．首先对于信号中的每个点选择一个固定的窗宽大小并计算窗宽内的信号幅值的方差，作为该点的方差．将这些点对应的方差存储在一个与该信号同样大小的向量当中．最后将方差最大的点赋予最小的窗宽半径，方差最小的点取最大的窗宽半径，其余的点依据每个点的方差大小来决定去噪时的相似窗口半径大小，计算公式如下：

其中${f_{\max }}和{f_{\min }}$分别为相似窗口半径的最大值和最小值，其值一般分别取60和40．$\sigma _{\max }^2$和$\sigma _{\min }^2$分别为所有点的方差最大值和最小值．下图显示了它们的对应关系：

由于CPMG信号前后两端特性的不同，去噪过程中应该合理的选取搜索窗宽度．如果取值较小，则无法找到足够的相似点以致影响去噪效果；取值过大则会增加计算量，同时也会将尾端的噪声部分引入到前段信号的去噪过程，对滤波结果造成干扰．通过对比大量实验结果，我们发现当搜索窗宽度取从首点幅度到其值的30％之间点的个数时降噪结果较优．

现将本文算法步骤归纳如下：

(1) 根据输入CPMG信号，计算出从信号首点开始，到回波幅度下降至首点幅度的30%之间的点的个数，将该值作为搜索窗宽度；

(2) 计算出输入信号中各点的数据方差的大小，并按照(6)式计算出各点进行NLM滤波算法时的相似窗口半径的大小；

(3) 使衰减参数λ从小到大按照固定的步长变化，同时将其值与步骤(1)和(2)中求得的搜索窗宽度的大小和相似窗宽度的大小代入到(1)式和(3)式中对含噪声信号进行NLM滤波和计算对应的SURE值，最后可以得到一条随λ变化的SURE值曲线，将该曲线的最低点所对应的NLM算法的滤波结果作为最优的滤波结果．

算法流程如图所示：

5 实验结果与结论 5.1 仿真数据

本实验是在MATLAB 2012b环境下运行的，首先使用中心T₂ = 10 ms处的高斯峰进行正演，然后向其中添加一定方差的高斯白噪声，最终得到含噪声的仿真数据．对于包含不同方差噪声的仿真数据，滤波过程中{f_{\max }}和{f_{\min }}$均分别取值为60和40，同时原始NLM算法的相似窗口半径为50，反演部分使用了文献^[14]中提及的方法．表 2显示了在添加不同方差的噪声时滤波后信号的信噪比的对比情况．图 8显示了在不同方差噪声下原始NLM算法（图中虚线）、本文改进的NLM算法（图中点线）滤波后反演谱、原始含噪声数据（图中菱形）的反演谱和正演谱（图中实线）对比，其中图 8(a)、(b)分别为噪声标准差为15、20的对比情况，图 8中左上角的图片是图中原始NLM算法、改进的NLM算法反演谱和正演谱的顶端放大图．从表 2和图 8中可以看出，本文改进的NLM算法从信噪比和反演结果两个方面都要优于原始的NLM算法．

5.2 真实数据

为了检验上述参数选择策略和NLM改进方法在实际应用中的效果，采用实际采集到的多组数据进行处理，其中数据来源于上海纽迈公司的NMI-20低场核磁共振分析仪，样品为某浓度的Mn²⁺造影剂，TE为1.009 ms，数据点个数为15 000．由于篇幅限制，本文仅对其中一组数据的实验结果进行显示．此次滤波过程中和均分别取值为60和40，同时原始NLM算法的相似窗口半径为50．图 9显示了原始数据和滤波后的效果图，图 9(a)为NMI-20采集的某样品数据，图 9(b)为本文改进算法滤波后的效果图．图 10显示了原始信号和滤波后信号的反演谱，其中图 10(a)中的3条曲线分别为原始采集的数据反演谱、原始NLM算法滤波后反演谱和改进的NLM算法滤波后的反演谱，图 10(b)为图 10(a)顶端局部放大后的图像．

对于该采集到的CPMG回波数据，其滤波前后反演谱对于时间轴的积分面积是相对不变的，与总信号量一定的事实相符合，同时可以通过反演谱线的半高宽判断反演效果的优劣，即若改进NLM算法比原始NLM算法的反演谱宽度窄且峰值高就可判断改进算法的优越性．为了进一步对比算法的效果，对图 10(a)中反演谱的顶端进行放大，从放大后的图 10(b)中可以明显的看到改进算法比原始算法对应的反演谱的峰值高，且对图 10(a)中反演谱的下半段进行放大可知改进算法比原始算法反演谱的宽度窄，因此在实际应用中改进的NLM算法的滤波效果优于原始NLM算法降噪效果，经过本文算法去噪之后，可以从反演谱中对样品的T₂值进行更加精确的定界．

6 总结与展望

近年来，低场核磁共振在实际检测技术当中发挥着不可替代的作用，但是其接收到的信号微弱和信噪比低的缺陷影响了低场核磁共振的应用，因此其降噪算法的研究必须受到重视．NLM算法在图像降噪领域中应用较为成熟，本文利用该算法的原理将其应用在一维CPMG信号的降噪中．此外，本文利用SURE方法对NLM算法中的参数选取部分提供了依据，同时根据CPMG信号的特点提出了相似窗宽度可变的NLM算法，通过对仿真数据和真实数据的实验证明了本文提出的改进NLM算法的优越性．NLM算法在图像降噪方面能取得良好的效果，但其在一维信号降噪中的研究较少，本文对其在一维信号中的应用提供了一定的参考价值．