0 引言

利用核磁共振(magnetic resonance sounding, MRS)技术探测地下水是20世纪80年代国际上发展起来的一种地球物理新方法^[1]。过去几年，该技术在正演模型^[2]和反演解释^[3-5]方面取得了显著的进展，在水文地质勘查方面的应用也进一步拓展。这种在地面直接探测地下水中氢核丰度的技术，其基本原理是利用人工产生的磁场激发地下水中的氢核形成宏观磁矩；这一宏观磁矩在地磁场中进行旋进运动，当激发场停止后，氢核自旋产生弛豫现象；再用地面线圈记录宏观磁矩进动产生的电磁信号，即可判断地下水是否存在^[6-10]。该电磁信号E(t, q)即为核磁共振信号，一般可用下式来诠释

其中：t为采样时间；q为激发脉冲矩；E₀为与地下水含量成正比的MRS信号初始振幅；T₂^*为含水层介质平均孔隙度的平均弛豫时间；f_L为由测量地点地磁场决定的拉莫尔频率；φ₀为与含水层导电性相关的初始相位。E₀、T₂^*、f_L、φ₀即为MRS信号的四个关键参数。MRS信号关键参数能否被准确提取依赖于测量地点的电磁噪声水平^[11]。当信噪比较低时，MRS信号完全淹没在电磁噪声中，难以获得四个关键参数。尽管核磁共振地下水探测方法具有多种优势，但由于MRS信号极其微弱，对周围环境中的电磁干扰敏感成为其应用的主要瓶颈。

影响信号质量最严重的两类电磁噪声是工频谐波噪声和尖峰脉冲噪声。工频谐波噪声来源于电力线、发电机、变压器、工业及生活用电等，频率固定为50或60 Hz的整数倍，或某一频率的单频干扰，其特点是频率固定，且持续时间长。尖峰脉冲噪声主要由太阳磁暴、雷暴或任何物体突然放电产生，一般持续时间很短，干扰幅度较大，高出信号幅度数倍，干扰频率范围较宽。

核磁共振信号的噪声压制近几年成为非常重要的研究热点。早先的核磁共振地下水探测仪采用单一线圈模式^{[1, 12-15]}，发射线圈和接收线圈分时复用同一线圈。普遍采用的方法是利用采集数据的多次叠加^[11]来提高信噪比。即当噪声强度大时需要叠加的次数多，这将增加测量的时间，导致仪器的工作效率降低。为此，Legchenko等^[16]提出采用正弦对消法、区块对消法和陷波滤波法滤除工频谐波干扰；然而，当测量地点的拉莫尔频率与工频谐波的干扰频率较接近时，MRS信号在陷波后会产生严重失真，降低信号的可靠性。在探测方法方面，通过合理铺设“8”字形线圈^[17-18]可有效抵消噪声干扰。这种方法在噪声均匀分布的环境中应用效果较好，但其激发场的分布会随着线圈形状的改变而发生变化，影响激发效果。此外，铺设“8”字形线圈探测时探测深度将减小为原来的一半，需要增加线圈尺寸才能达到传统方形或圆形线圈的探测水平。由于最新的仪器可以通过铺设多个接收线圈实现信号和噪声的同步采集，因此通过参考消噪方法可明显提高信号的信噪比^[19-20]。参考消噪方法需要用一个探测线圈来测量含有噪声的MRS信号，另外一个或者多个参考线圈在远端同时测量不包含信号的噪声数据，其中最关键的一步是获得使参考线圈数据与探测线圈数据相关的传递函数。2006年，Radic^[21]首次描述了频域上的滤波算法，提出使用远端垂直铺设的参考线圈采集空间磁场的水平方向分量，以此对探测线圈中的相关噪声进行抵消计算；但由于参考线圈感应空间磁场水平分量与探测线圈磁场垂直分量的相关性较低，给出的应用结果中仍有较高水平的残余噪声。2008年，Walsh^[19]研制多通道地面核磁共振找水系统，从多个通道上提取主通道探测信号与其他通道噪声信号构成新的数据序列，基于频域自适应消噪算法实现噪声滤除。2010年，Muller-Petke等^[20]提出了一种基于单参考通道的时域噪声预测滤波器，该滤波器可以自动计算参考通道数据与主通道数据的振幅与相位差，得到含有校正系数的传递函数，实验结果显示该方法可以完整地保留MRS信号，也能对以工频谐波为主的电磁干扰进行压制；然而，通过该方法滤波后的信噪比还有待提升。2012年，Dalgaard等^[22]采用时域自适应滤波器实现参考消噪方法，同时发现类似于尖峰脉冲的噪声干扰影响了通道间的数据相关性，需要先剔除尖峰脉冲干扰。2014年，Muller-Petke等^[23]开展了时域参考消噪方法和频域参考消噪方法的系统研究，认为频域参考消噪方法的滤波效果优于时域参考消噪方法。另外，法国推出的NUMIS^Poly系统也是针对存在强电磁干扰的探测环境，采集多个远端的参考数据，抵消信号中的环境噪声；德国推出的NMR-MIDI还加入了远端垂直方向的参考线圈，进一步消除环境噪声影响；我国研究人员采用宽频带滤波器，通过自适应参考消噪技术，在相同的电磁噪声环境下，利用较少的叠加次数获取了较高信噪比的地下水MRS信号^[24]。

此外，通过分析噪声特性，采用建模法去除特定的电磁干扰开辟了磁共振噪声压制方法的新领域。2013年，为解决多源噪声干扰下自适应噪声对消方法去噪能力有限的难题，Larsen等^[25]提出建模法去除工频谐波噪声，该方法与正弦对消法相似，通过建模构建通带内的工频谐波干扰以提高信噪比。2015年，Larsen^[26]使用两个级联的二阶带通滤波器对尖峰脉冲噪声建模，以剔除来源于电网的尖峰脉冲干扰。

吉林大学地下水磁共振探测技术课题组总结国内外研究现状，深入分析噪声的产生机理，针对影响MRS信号质量最严重的两类电磁干扰工频谐波噪声和尖峰脉冲噪声，取得了一系列的噪声压制研究进展。

1 工频谐波噪声压制方法研究进展 1.1 实时参考消噪技术

使用吉林大学自主研制的核磁共振多通道宽频带接收系统^[27]，结合后续数据处理可压制工频谐波噪声；然而，通带内噪声干扰严重容易淹没微弱的MRS信号，导致放大器饱和。针对多通道宽频带接收系统受强电磁干扰放大器易饱和的问题，吉林大学地下水磁共振技术课题组正在开发基于现场可编程门阵列(field-programmable gate array, FPGA)的实时参考消噪技术，将后续的自适应消噪处理环节移植到前端的FPGA中，通过铺设参考线圈采集探测点附近的环境噪声，在信号和噪声未知的情况下，采用变步长最小均方(least mean square, LMS)算法自适应跟踪学习去除相关噪声，在系统不饱和的前提下尽可能提高信噪比，同时可现场实时检测MRS信号。

由于定步长LMS算法的收敛速度与稳态误差相互矛盾，并且对输入信号信噪比要求严格；因此，课题组决定采用变步长LMS算法^[28-34]实现实时参考消噪，其本质是利用误差因子的e指数形式作为步长参数，进行自适应滤波。变步长表达式为

式中：m是离散的采样时间；μ是控制算法稳定性和收敛性的步长因子；e是误差信号；α是控制函数形状的常数；β是控制函数取值范围的常数。

理想的步长参数是：当e(m)较大时，选择较大的步长因子加快收敛速度；当e(m)较小时，选择较小的步长因子减小稳态误差。该e指数变步长表达式初始时刻e(m)较大，步长大小主要取决于参数β，而达到稳态时e(m)很小，步长趋于0；符合理想步长参数要求。通过适当调整参数β，可使算法快速稳定收敛。含有指数因子的变步长LMS算法不仅可以快速稳定地收敛，而且具有较好的鲁棒性；同时对输入信号的信噪比要求降低，是实现实时参考消噪算法较好的选择。

在核磁共振宽频阵列式接收系统^[27]中嵌入变步长LMS滤波模块，以实现实时参考消噪，其原理如图 1所示。前端处理模块由低噪声前放和宽频带带通滤波器构成，实现对信号的前置放大和初步滤波；使用24位超低噪声模数转换器(analog to digital converter, A/D)采集模块对信号进行同步采集并送至FPGA，在FPGA中利用变步长LMS滤波器对含有工频谐波噪声的核磁信号进行处理；最后通过以太网接口模块对处理完的数据进行打包送至上位机，并在Labview中实时显示。

图 2和图 3分别为30阶变步长LMS滤波器理论模型仿真和硬件实时测试结果。测试结果表明，主通道中被工频噪声淹没的核磁信号自适应滤波后获得了自由衰减的MRS信号。由于FPGA的电路效应，滤波器输出会产生短暂延时，输出不能立即随输入变化，出现竞争冒险，产生高频毛刺信号；可通过在FPGA的输出端并联小电容去除该毛刺信号。此外，硬件实时测试结果表明，当参考通道信号和主通道信号在0.20 s后为0时，消噪后的信号经电路延时后也为0；再次证明变步长LMS滤波算法在FPGA中是实时实现的。对基于理论模型的软件仿真结果和基于FPGA的硬件实测结果进行频谱分析，硬件实测结果虽然没有理论仿真效果好；但总体而言，在FPGA中运行变步长LMS算法可以很大程度地实时抑制磁共振信号中的工频谐波干扰，有效防止系统饱和，最大限度地从环境噪声中提取有用的核磁共振信号，且不发生畸变和衰减。

1.2 基于独立分量分析的工频谐波噪声压制技术

陷波法^[16]和参考消噪方法^[19-20]均是基于信号和噪声在时频域上具有不同特征加以分离；然而，噪声的无规律性、不稳定性以及混合机理的不确定性限制了算法的有效应用。独立分量分析(independent component analysis, ICA)方法^[35-40]作为一种有效的盲源分离手段，其基本思路是，在源信号、噪声以及混合机理均未知的情况下，假设一组相互独立的源信号经过线性混合系统后得到一组混合信号，通过选择某种独立性判据，构造一个多通道混合信号的目标函数，再利用合适的优化算法调节解混矩阵以得到最优解，最终得到各分量之间相互独立的输出信号，以分离源信号。主要包括四个步骤。

1)首先对观测信号进行去均值处理，使数据中心化：

式中：x是观测信号；E{·}是求观测信号的均值。

2)对零均值信号进行白化处理，得到数据z(t)。该数据具有单位方差，数据之间的相关性较低。白化处理的表达式为

式中，D和E分别是观测信号的协方差矩阵特征值和特征向量。

3)利用基于负熵最大化原理和牛顿迭代求解解混矩阵w，得到独立成分y：

4)数据重构，保留MRS信号的独立分量，屏蔽噪声分量，恢复信号幅度：

吉林大学地下水磁共振探测技术课题组利用自主研制的磁共振地下水探测仪^{[9, 41]}在长春市烧锅镇开展了基于独立分量分析技术抑制工频谐波噪声的试验研究。采用边长为100 m的方形单匝线圈，激发频率为2 326 Hz，得到相同激发脉冲矩下三个含噪声磁共振信号的时域波形图(图 4ac)；采用基于独立分量分析技术抑制工频谐波噪声后的信号如图 4df所示。为了验证算法的有效性，将独立分量分析方法与陷波法进行比较。图 4gi是陷波法去除工频谐波噪声后信号的时域波形，图 4jl是消噪前后信号的频域对比结果。经对比分析得到结论：独立分量分析技术对工频谐波噪声的抑制效果优于陷波法。

2 尖峰脉冲噪声压制方法研究进展 2.1 基于能量运算的磁共振信号尖峰噪声抑制技术

通常的磁共振信号尖峰噪声处理方法是基于统计原理找到数据中的“粗大误差”，然后进行删除或修正，正如广泛应用的统计叠加方法^[42]。然而，针对干扰幅度较小的尖峰噪声，统计叠加法很难剔除。因此，利用基于能量运算的尖峰噪声抑制技术，通过计算信号的能量，可有效检测尖峰噪声并突出不易识别的小幅度尖峰噪声。主要包括三个步骤。

1)计算单次采集信号的能量，将信号由时域转化到能量域；再将计算后信号的能量通过低通滤波器。设x_n为采样信号，其能量E_n可表示为

2)经过上述能量运算后，尖峰噪声被突出，此时需要选取适当的阈值进行尖峰识别并剔除。基于中位数的绝对偏差(median absolute deviation, MAD)法^[43]是一种非常有效的统计学方法，定义为

式中：x_i为该时刻某次采集的数据点，x_j为某一时刻所有单次采集的数据列；median_j{·}是求该数据列的中位数。计算MAD后，阈值为median{x}+k·MAD(K是常数)。根据尖峰噪声和信号在能量域中的幅值关系，阈值以能有效标记出被突出信号的能量为标准进行选取。如图 5所示，尖峰噪声在能量域中被突出，使用基于中位数的绝对偏差法确定的阈值既跟踪了信号的整体变化趋势，又有效地将幅度较小的尖峰噪声标记出。

3)将能量域中的信号与阈值进行比较，超出阈值的部分视作尖峰噪声，用该时刻无尖峰噪声数据列的中位数进行替代，以达到剔除尖峰噪声的目的。

为了验证算法的实用性，在长春市郊区烧锅镇进行了野外实验。使用核磁共振地下水探测仪^[9]进行信号采集。当地的拉莫尔频率为2 326 Hz，图 6a中的蓝色线为仪器采集的原始信号时域波形图。从图 6a中可以看出：该信号含有较多的尖峰噪声，并且从信号的时域波形上不能明显看出MRS信号指数型的衰减趋势。经计算，原始信号的信噪比为-37.709 0 dB。

依次操作能量运算法的三个步骤，对图 6a所示的实测数据进行尖峰噪声剔除，得到消噪后的时域信号，如图 6a中的红色线所示。可见，原始信号中的尖峰噪声被明显剔除，消噪后的信号呈现指数形式衰减。再从频域上进行比较，结果是消噪后尖峰噪声的频率分量被压制，突出了信号的频率，如图 6b所示。经计算，消噪后信号的信噪比为-29.138 0 dB，剔除尖峰噪声后信噪比提高了8.571 0 dB。

2.2 基于同步压缩小波变换和自寻优非线性阈值的补偿技术

基于统计学原理识别尖峰噪声，进行数据替代时很容易损失MRS信号的有效信息，可采用基于同步压缩小波变换^[44-46]和自寻优非线性阈值补偿技术^[47]抑制尖峰噪声的方法同时重构MRS信号，提高信噪比。

实现同步压缩小波变换，首先要设定所要采用的小波母函数和所要划分的压缩小波域的长度，对所划分不同压缩小波域内的磁共振信号求解压缩小波系数。同步压缩小波变换表达式为

式中：T_x(ω_l, b)为同步压缩小波系数；W_x(a_i, b)为连续小波变换的小波系数；a_i为尺度因子a的离散值；b为平移因子；ω_l为中心频率，在一个连续区间内变化；Δω=ω_l-ω_l-1; n_v为自定义参数，决定了尺度的数量。

对同步压缩小波变换后信号用自寻优算法搜寻尖峰噪声在信号中的位置，并将其剔除。自寻优算法表示为

式中：为去噪后的同步压缩小波系数；ω_{a, b}为压缩小波系数；T为压缩小波系数中的最大值；n为任意正整数，是一个可调的参数。

运用非线性阈值选取算法对尖峰噪声所在位置的MRS信号进行补偿，确定所需非线性阈值函数：_{j, k}=sign(ω_{j, k})×[(ω_{j, k}-k)·(ω_{j, k})²+λ²/(2λ+k)k],

式中：ω_{j, k}为原始数据的压缩小波系数；_{j, k}为去除尖峰噪声后获得的新的压缩小波系数；λ为预先设定的阈值；σ为可调参数。

最后，对经过自寻优算法以及非线性阈值选取算法处理后的MRS信号进行压缩小波重构，得到最终的剔除尖峰噪声后的MRS信号。

为了证明算法的有效性，在长春市净月潭进行试验。使用核磁共振地下水探测仪^[9]进行信号采集。当地的拉莫尔频率2 330 Hz，图 7为应用基于同步压缩小波变换和自寻优非线性阈值补偿技术后重构的MRS信号。

3 磁共振地下水探测噪声压制技术展望

伴随多通道磁共振探测系统的引入，使用宽频带数据采集，为进一步提高信噪比提供了更多的潜能，目前已经在地下水磁共振信号的工频谐波噪声和尖峰脉冲噪声压制方面取得了有效的研究成果和进展。然而，较宽的通频带会引入更多的随机噪声，若采用传统的叠加方法，需要增加叠加次数，且滤波效果有限，严重降低探测效率。因此，亟需研究核磁共振信号随机噪声的压制方法。

由于随机噪声在时间域中会淹没微弱的核磁共振信号，影响核磁共振信号的检测与提取，在频率域中，其功率谱分布在整个通频带范围内，与有效信号的功率谱重叠，难以在时间域或频率域内有效压制随机噪声；因此，如何把含随机噪声的信号投影到某一特定空间，进而有效区分信号和噪声，利用信号和噪声的不同特性压制随机噪声，将是磁共振地下水探测噪声压制方法未来具有挑战性的发展方向。