0 引言

与其他地球物理方法相比，地面磁共振技术(Surface Nuclear Magnetic Resonance，简称SNMR)可以直接探测地下水，具有定性定量的特点，已广泛应用于地下水勘探领域，为缓解旱灾提供了必要的技术支持(林君等，2011; Behroozmand et al., 2015; 蒋川东等，2018).该方法也称为磁共振测深(Magnetic Resonance Sounding，简称MRS)，主要利用磁共振弛豫原理，基于地下水中的氢原子在磁场中会发生进动，因此只要检测到磁共振信号即可证明探测点地下水的存在.通过磁共振响应的关键参数——初始振幅、平均弛豫时间和初始相位，可以分别评估出地下含水层位置、含水量大小、介质孔隙度及地层渗透率等地质水文信息(林婷婷等，2013; 陈斌等，2014).由于SNMR信号极其微弱，只有纳伏级，易受噪声干扰，因此如何从背景噪声中检测出有效信号一直都是该技术领域的研究热点(林君等，2016; Lin et al., 2017).

近些年，国内外相关学者从软硬件两个方面对滤波技术和关键参数提取方法进行了大量研究.地面磁共振探测仪由最初分时复用同一线圈发射和接收信号发展到如今的多通道同步采集系统.相应地，SNMR信号处理技术也从最开始的“8”字线圈法(Chalikakis et al., 2008)、区块对消法(Dalgaard et al., 2012)等发展至今，形成了相对完整的SNMR信号处理流程(Behroozmand et al., 2015)(如图 1所示)：1)消除尖峰噪声；2)消除相干工频谐波噪声；3)消除随机噪声；4)拟合包络曲线提取有效信号.首先，采用非线性能量算子法(Dalgaard et al., 2012；万玲等，2016)识别尖峰噪声并剔除，减小尖峰噪声对通道之间数据相关性的干扰.然后，利用主通道与参考通道之间噪声的相关性，基于自适应参考消噪算法(Walsh, 2008; Jiang et al., 2011; Dalgaard et al., 2012)，压制工频谐波.Larsen等(2014)和田宝凤等(2018)根据工频噪声的数学特征，提出了一种工频噪声建模方法，可以有效减弱SNMR信号中的工频噪声成分.对于随机噪声干扰，Lin等(2018)、林婷婷等(2018)提出了一种时频峰值滤波方法代替传统的统计叠加方法，利用解析信号时频分布沿瞬时频率集中的特性，在大幅度消除随机噪声干扰的同时，提高了信号探测效率.最后，通过包络拟合提取关键参数，Legchenko和Valla(1998)使用全局非线性最小二乘参数拟合方法来提取信号，以减少误差，提高信号准确度.在此基础上，Müller-Petke等(2016)总结了提取信号的同步检测方法，该方法利用希尔伯特变换得到解析信号，与参考信号相乘后应用低通滤波器进行处理.然而，当噪声频率与SNMR信号的拉莫尔频率相同时，也即是同频噪声干扰情况下，现有的数据处理方法无法正确分辨噪声成分与SNMR信号，在消除同频噪声时不可避免地会对MRS信号进行消减，导致信号信息丢失和参数提取错误.Liu等(2018)提出采用分时段建模法，第一阶段利用多通道仪器同步采集噪声，确定同频噪声成分在主通道和参考通道之间的振幅比和相位差；第二阶段在探测工作时，结合第一阶段得到的通道间噪声振幅比和相位差估算此时主通道中的同频谐波并进行消减.但是该算法是基于电力线基频和谐波噪声的幅值相位在短时间内为常量的假设之上提出的.在实际采集和处理数据的过程中，由于工频噪声非线性引发高阶谐波，以及电力线上负载的瞬时变化，主通道与参考通道线圈获取的同频噪声幅度和相位往往是非线性的，而且噪声采集的实际效果还依赖于野外实验过程中参考线圈的摆设方式，其操作过程复杂且耗费时间.因此亟需一种有效压制同频噪声干扰的SNMR信号检测方法，为提取精准的关键参数提供技术支撑.

结合SNMR信号的特殊方程，本文基于相关建模检测技术(Correlation Modeling Detection Technology，简称CMDT)提出一种新的检测方法，可以在含有同频噪声干扰的数据中提取出高精度的关键参数.该方法根据相关原理，将原始数据与同频参考信号相乘，对其结果进行积分，也即低通滤波；再通过数学建模消除滤波结果的干扰项和同频噪声项，获取只含有SNMR信号关键参数的线性方程组；然后，联立方程组提取高精度关键参数.本文其余内容如下：首先，分析了同频噪声干扰对提取SNMR信号关键参数的影响；其次，介绍了本文提出的基于相关建模检测的算法原理和流程；最后，把已知的SNMR信号加入仿真噪声数据和实测野外噪声数据对新方法进行定性定量测试验证，并分析其性能特点及限制条件.

1 同频噪声对磁共振信号影响分析

核磁共振地下水探测方法的原理是基于地下水中的氢原子具有核子顺磁性，在激发电流作用下会产生与地磁场方向垂直的磁化强度，发生能级跃迁.当撤去激发电场时，氢原子在地磁场作用下作旋进运动，并释放能量回落到低能级状态，在接收线圈中感应出具有拉莫尔频率的SNMR信号：

(1)

式中E(t, q)为感应电动势，与数据采集时间t和激发脉冲矩q相关，改变激发脉冲矩q大小可探测不同深度的地下水；E₀(q)是SNMR信号的初始振幅，反映了氢质子的宏观数量，表征含水量大小；T₂^*(q)是氢原子的弛豫时间，与含水层岩石孔隙度有关；θ(q)是初始相位，与地下介质电导率相关.f_L是氢原子进动的拉莫尔频率，由地磁场强度决定，是对氢原子施加交流电场引发磁共振现象的频率条件.

在野外实际探测中，SNMR信号十分微弱，易被环境噪声淹没，虽然SNMR技术已经形成了特有的信号处理策略和流程，但是当存在与信号频率同频的噪声时，由于同频噪声在时域和频域均与信号重叠，因此现有算法难以正确分辨噪声成分与SNMR信号，在消除同频噪声时不可避免地会对SNMR信号进行消减，导致信号信息丢失和参数提取错误，使反演得到的含水量信息产生偏差.例如，在SNMR信号拉莫尔频率处加入一个同频噪声，如图 2所示.图 2a中蓝色曲线表示SNMR信号，其初始振幅E₀=200 nV，弛豫时间T₂²=150 ms，拉莫尔频率f_L=2350 Hz，初始相位θ=30°.其中，同频噪声幅值为50 nV，加入噪声后的时域曲线为图 2b所示，可以看出在噪声干扰下SNMR信号时域曲线产生了一定程度的偏差，尤其晚期信号，完全淹没在噪声中.为了观察信号和噪声的频谱特征，在图 2c中用蓝色曲线和红色曲线分别描绘了加入同频噪声前后SNMR信号的频域曲线，可见在加入同频噪声之后，噪声与SNMR信号重叠，信号频谱幅度增大.图 2d所示是对理想的SNMR信号和同频噪声干扰下的SNMR信号进行拟合得到的曲线，以便于根据拟合结果提取信号特征参数，结果表明同频噪声干扰下提取的信号发生畸变，无法直接提取有效关键参数.

2 基于CMDT的同频噪声抑制算法 2.1 含有同频噪声的SNMR全波信号检测策略

由于SNMR信号极其微弱，需要结合当前的数据处理流程进行去尖峰、工频谐波和随机噪声.在此基础上，再对含有同频噪声干扰的SNMR全波信号进行处理，抑制同频噪声干扰并提取有效信号参数.相关检测技术(Correlation Detection Technology，简称CDT)可以从低信噪比背景中提取有效信号，广泛应用于声学、地学、光学等弱信号检测领域.

假设含有同频噪声的全波磁共振信号可表示为

(2)

其中t为采样时刻，s(t)为SNMR全波信号，n_c(t)表示同频噪声.本文基于CDT理论，并结合SNMR信号特征，提出基于相关建模的磁共振信号同频噪声抑制算法，进而从含有同频噪声的数据中提取出有效信号s(t).

基于相关建模的磁共振信号同频噪声抑制算法流程如图 3所示，其过程可分为三个步骤：首先将包含同频噪声的SNMR全波信号数据做相关检测处理，即原始数据分别与同相参考信号和正交参考信号相乘实现频谱迁移，再进行积分变换或低通滤波获取低频分量；其次，进行相关建模，即对两相通道的输出结果求导以消除同频噪声干扰项，再与指数因子相乘后求导变换，获取只含有SNMR信号参数的方程；最后，通过两个通道输出结果，建立方程组，利用最小二乘指数拟合方法提取SNMR信号的所有关键参数.整个过程既能压制同频噪声干扰，又能提取出信号参数.

2.2 基于相关建模的同频消噪算法

在某一个脉冲矩激励下产生的磁共振信号数学表达式可表示为e(t)=E₀e^(-t/T₂*) cos(ω_Lt+θ)，同频噪声信号n(t)=V_ncos(ω_Lt+θ_n)，两者混合后的信号为x(t)=s(t)+n_c(t).设定相关检测的同相通道参考信号是r(t)=cos(ω_Lt)，正交通道参考信号为r(t)=sin(ω_Lt).其中，噪声、参考信号与磁共振信号为同一频率，即都为拉莫尔角频率ω_L=2πf_L.

首先，将包含同频噪声的SNMR全波信号数据x(t)做相关检测处理，下面以同相通道为例进行说明.x(t)与参考信号乘积后的结果可表示为

(3)

其中，u_low(t)=0.5E₀e^(－t/T₂*)cos(θ)+0.5V_ncos(θ_n)为包含了信号和同频噪声的低频分量，u_high(t)=0.5E₀e^(－t/T₂*)cos(2ω_Lt+θ)+0.5Vcos(2ω_Lt+θ_n)]为两者的高频分量，且频率为拉莫尔频率的2倍.再对u(t)做积分变换或低通滤波，为了便于计算，设定积分变换的传递函数为G(s)=ω_c/(s+ω_c)，其中ω_c为积分变化因子，也即低通滤波器的截止角频率，为了消除u(t)高频分量的影响，设定积分变化因子为极小的值(ω_c < 1 rad).故，忽略高频分量后，积分变换后的结果可表示为

(4)

式中，C是包含同频噪声参数的常量.

其次，对相关检测处理的结果进行建模，通过数学模型变换，包括求导、乘以指数因子和再求导，以消除同频噪声参数，并获取只含有信号参数的方程.为了消除同频噪声参数的影响，对u₀(t)求导，求导后的结果为

(5)

然而，式(5)中指数项依然为干扰项，式(5)乘以质数因子e^(ω_ct)，再对乘积后的结果进一步求导可得

(6)

同理，通过正交通道的相关建模处理可以得到另一个方程

(7)

最后，联立公式6和7，建立方程组，并利用最小二乘指数拟合方法可以分别提取SNMR信号的关键参数，包括初始振幅E₀、弛豫时间T₂^*和初始相位θ.由上述推导可知，运用CMDT算法处理含有同频噪声的SNMR全波信号时，可以有效抑制同频噪声的干扰，并获取原始SNMR信号的关键参数.

3 仿真实验

当采集的数据中含有同频噪声时，经传统SNMR数据处理流程，难以提取出有效的信号参数.为了解决同频噪声干扰问题，本文提出了CMDT算法，通过该算法对同频噪声的抑制情况和提取信号参数的精度验证新算法的性能.并对信噪比、初始振幅提取误差、弛豫时间提取误差和初始相位提取误差进行量化评估，来衡量MRS信号提取结果的准确性.信噪比的单位是dB，定义式如下所示：

(8)

初始振幅提取误差(Initial Amplitude Extracted Error，IAEE)定义为

(9)

式中，为初始振幅的拟合值.

弛豫时间提取误差(Relaxation Time Extracted Error，RTEE)定义为

(10)

式中，为弛豫时间的拟合值.

初始相位提取误差(Initial Phase Extracted Error，IPEE)定义为

(11)

式中，为初始相位的拟合值.

3.1 包含同频噪声的SNMR信号处理仿真实验

通过模拟信号和噪声仿真CMDT算法的性能，仿真模型参数如下：CMDT算法中积分变化因子为ω_c=0.01，SNMR信号的初始振幅E₀=200 nV，弛豫时间T₂²=100 ms，拉莫尔频率f_L=2350 Hz，初始相位θ=30°，加入幅值V_n=100 nV的同频噪声，计算可知信噪比SNR=-14.18 dB.加入同频噪声的SNMR全波信号时频域波形如图 4所示，由于噪声较大，时域上的早期信号幅度较大依然可以看出磁共振信号衰减趋势，而晚期信号完全被淹没，频域也完全重合，传统方法无法直接提取.

根据上一章节介绍的CMDT算法流程对信噪比为-14.18 dB的原始数据进行处理：首先将包含同频噪声的SNMR全波信号数据做相关检测处理；其次，通过求导、指数因子乘积、二次求导进行相关建模变换，获取只含有SNMR信号参数的方程；最后，建立方程组，并利用最小二乘指数拟合方法提取SNMR信号的所有关键参数.提取的关键参数拟合结果为：= 202.42 nV、= 100.03 ms、= 29.98°，与原始信号相比，三个参数的提取误差分别是IAEE=1.21%、RTEE=0.03%、IPEE=0.067%，提取误差较低，提取的信号包络与原始SNMR信号包络曲线如图 5所示，基本吻合，验证了CMDT方法压制同频噪声干扰的有效性.

下面我们进一步探究在不同信噪比情况下CMDT算法的性能.通过改变同频噪声幅度V_n调节SNR，并采用同样的模拟SNMR信号参数进行仿真.仿真结果如表 1所示，噪声幅度由20 nV到1600 nV，信噪比由6.88 dB到-31.17 dB，所有参数的最大提取误差不大于1.22%，新算法都能提取出高精度的信号参数.这一仿真结果表明采用本文提出的CMDT算法能够压制同频噪声干扰的同时，提取出高精度的磁共振信号，而不受同频噪声强度的限制.

3.2 含有随机噪声和同频噪声的SNMR信号检测实验

野外工作时，经过传统SNMR数据处理后依然会有部分残余噪声，主要是随机噪声，因此为了模拟实际情况，进一步验证CMDT算法性能，我们在SNMR信号数据中加入同频噪声和随机噪声.即在随机噪声的干扰下，测试CMDT算法处理含有同频噪声数据的效果.基本仿真模型参数：CMDT算法中为积分变化因子ω_c=0.01，SNMR信号的初始振幅E₀=200 nV，弛豫时间T₂^*=100 ms，拉莫尔频率f_L=2350 Hz，初始相位θ=30°，加入幅值为20 nV的同频噪声和不同信噪比的高斯白噪声.以信噪比SNR=-10.12 dB为例，进行仿真说明，该信噪比的噪声包含了同频噪声和高斯白噪声，其中高斯白噪声与信号的信噪比SNR=-9.96 dB.如图 6所示，我们要处理的原始数据包含信号和两种噪声，从时域波形(图 6a)可以看出SNMR信号幅度小于整体原始数据，完全淹没在噪声中，从频域波形(图 6b)可以看出拉莫尔频点的信号与同频噪声叠加在一起，无法采用传统的检测方法直接提取信号，且随机噪声在每个频点的强度不大，但分布在拉莫尔频率附近的整个频带内.从原始数据的三维时频谱图(图 6c)也可发现同频噪声和随机噪声在时域和频域与信号完全重合，经典数据处理方法难以去除噪声提取精准的信号.

下面，采用本文提出的CMDT算法对上述原始数据进行处理，并以CMDT算法中的同相通道为例进行附图说明.首先将包含同频噪声、随机噪声和SNMR信号的原始数据做相关检测处理，即原始数据与同相参考信号相乘实现频谱迁移，结果如图 7所示.时域结果(图 7a)表明该步骤输出仍含有大量噪声，频域结果(图 7a)表明信号和同频噪声的频域由拉莫尔频率迁移到低频和2倍拉莫尔频率，三维时频谱图(图 7c)低频分量含有更大的幅度，但噪声和信号在时频域上仍是交叉重叠在一起的.

然后，再进行积分变换或低通滤波获取低频分量，结果如图 8所示.时域结果(图 8a)表明该步骤得到一条类指数曲线，且幅度远小于处理之前的幅度，频域结果(图 8a)表明高频成分都被滤除、仅剩下低频分量，三维时频谱图(图 8c)也进一步证明了相关检测滤除噪声的特点，随机噪声基本被滤除干净，只剩下同频噪声和信号.

最后，根据上一章节介绍的相关建模方法，即对两相通道的输出结果求导以消除同频噪声干扰项，再与指数因子相乘后求导变换，获取只含有SNMR信号参数的方程；同理，获取正交通道的输出方程，建立方程组，利用最小二乘指数拟合方法提取SNMR信号的所有关键参数.

提取的关键参数拟合结果为：=200.34 nV、= 108.49 ms、=27.58°，与原始信号相比，三个参数的提取误差分别是IAEE=0.17%、RTEE=8.49%、IPEE=8.06%，提取误差均小于10%，提取的信号三维时频图如图 9所示，验证了CMDT方法在含有同频噪声和随机噪声干扰时可以有效检测SNMR信号.

然而，与仅含有同频噪声的数据相比，引入随机噪声后，CMDT算法对信号参数的提取误差会增大，RTEE和IPEE尤其显著.因此，为了进一步论证在不同强度随机噪声干扰下CMDT算法性能，我们进行下面不同信噪比的四组对比实验.实验参数与上面相同，仅改变随机噪声的强度，以调节信噪比.噪声强度定义为：SNR≥-10 dB为弱噪声、-10 dB>SNR≥ -15 dB为普通噪声、SNR < -15 dB为强噪声.设定四组不同噪声强度数据进行实验，信噪比分别为SNR=-3.46 dB、SNR=-10.12 dB、SNR=-13.05 dB和SNR=-16.94 dB.CMDT算法对这四组原始数据处理后，提取的信号结果如表 2所示.

当SNR=-3.46 dB时，CMDT算法提取的关键参数拟合结果为：=197.59 nV、=102.37 ms、=27.96°，与原始信号相比，三个参数的提取误差分别是IAEE=1.21%、RTEE=2.37%、IPEE=6.80%，提取误差较低，能够获取有效信号.

当SNR=-10.12 dB时，与之前实验相同，CMDT算法提取的关键参数拟合结果为：= 200.34 nV、= 108.49 ms、=27.58°，与原始信号相比，三个参数的提取误差分别是IAEE=0.17%、RTEE=8.49%、IPEE=8.06%，提取误差均小于10%，能够获取有效信号.

当SNR=-13.05 dB时，CMDT算法提取的关键参数拟合结果为：₀=232.28 nV、=161.16 ms、=22.09°，与原始信号相比，三个参数的提取误差分别是IAEE=16.14%、RTEE=61.16%、IPEE=26.36%，各参数均有较大提取误差，尤其是平均弛豫时间与真实值偏离严重，在此信噪比情况下，不能获取有效信号.

当SNR=-16.94 dB时，CMDT算法提取的关键参数拟合结果为：= 505.15 nV、=512.98 ms、=50.83°，三个参数的提取误差分别是IAEE=152.5%、RTEE=412.9%、IPEE=69.43%，与原始信号相比，提取的信号完全失真，在此信噪比情况下，无法获取有效信号.

实验结果表明，随着随机噪声强度增大、信噪比降低，CMDT算法检测含有同频噪声干扰的SNMR信号的精度也逐渐降低，且在信噪比低于-13.05 dB时无法获取有效的信号.

采用CMDT算法对四组不同信噪比的原始数据处理后，提取的信号包络与原始SNMR信号包络如图 10所示，随着信噪比降低，信号参数误差越来越大，提取的信号与原始信号偏离程度越发严重，甚至完全失真.

通过本章节实验的定性定量分析表明，本文提出CMDT算法在仅含有同频噪声干扰下，能够压制同频噪声干扰的同时，提取出高精度的磁共振信号，而不受同频噪声强度的限制；且存在随机噪声时，CMDT算法能够处理弱噪声干扰下的同频噪声和SNMR信号数据，并提取出有效信号；但当随机噪声强度较大时，CMDT算法难以提取有效信号.因此，在强噪声环境中，要通过经验和其他消噪方法进行噪声压制以提高信噪比，最后再采用CMDT算法进行处理，可有效压制同频噪声的同时，提取出高精度信号参数.

4 结论与展望

在强噪声干扰下如何压制噪声，检测出有效SNMR信号是当前磁共振探水领域的重点和难点，尤其是受同频噪声干扰时，当前的方法难以准确提取关键参数.本文针对这一难题，结合磁共振信号与同频噪声特征，提出了相关建模检测算法，通过相关检测滤除与拉莫尔频率偏离较远的噪声，采用模型变换实现同频噪声参数的抑制，并结合最小二乘指数拟合方法提取高精度SNMR信号.通过数值仿真实验对新方法进行了定性定量分析，结果表明CMDT算法可有效处理含有同频噪声的SNMR全波信号，且效果不受同频噪声强度的影响，验证了该算法的有效性.另外，在随机噪声干扰的情况下，CMDT处理含有同频噪声的SNMR全波信号的性能受随机噪声强度影响较大，当信噪比大于-10.12 dB时，仍可以获取有效的信号，但随着信噪比降低，提取误差增大，甚至无法检测到有效信号.因此，下一步将研究含有更强随机噪声的磁共振数据处理问题，并将本文提出的方法应用实际野外数据处理中.