2018, Vol. 61
磁共振(Magnetic Resonance Sounding，MRS)作为唯一一种直接定量表征含水量和孔隙大小的方法，已经成为探测地下水的重要手段之一(Behroozmand et al., 2015；Hertrich，2008).近年来，随着MRS理论方法和数据处理技术的迅速发展(Jiang et al., 2016；Grunewald et al., 2014, 2016；Dalgaard et al., 2012；Müller-Petke and Costabel 2014)，MRS的应用扩展到区域水资源调查、水文环境评价和灾害预警等领域(Jiang et al., 2017；Lin J et al., 2013).然而，仪器采集数据信噪比低和信号检测难度大仍然是制约MRS方法在城镇、村庄等人居环境和矿井、隧道等地下工程开展有效探测的主要因素(张荣等，2006；潘玉玲和张昌达，2000).
MRS方法中的背景磁场是地磁场，强度在(5~6)×10-5 T，与大于0.5 T的中高场核磁共振(Nuclear Magnetic Resonance，NMR)相比，信号非常微弱，通常只有几到几百纳伏(nV，10-9 V)，要求仪器具有较高的灵敏度.同时，MRS方法都是在室外进行测量，无法屏蔽外界环境噪声干扰，导致MRS仪器测量数据信噪比低，检测信号的难度大.最初，MRS仪器采用收发共线探测方式(Roy and Lubczynski, 2003)，由于仅有一个数据采集通道，且采样率和分辨率较低，很多滤波和叠加等消噪方法的效果受限(Legchenko and Valla, 2003；Strehl，2006). “8”字形线圈包含两个尺寸相等的线圈，反向相接可以抵消环境噪声，但由于噪声分布不匀，在不同测量地点消噪效果差异大，而且探测深度比常规线圈减小一半(Trushkin et al., 1994).多通道MRS仪器的出现(Walsh，2008；Lin T T et al., 2013)改善了单通道的缺陷，利用一个主通道和其他参考通道同步采集信号和环境噪声，利用自适应算法进行噪声抵消(Larsen et al., 2014；Dalgaard et al., 2012).但当多通道仪器用于二维和三维MRS探测时，所有的通道均被用于接收信号，无法进行参考消噪.因此，多通道MRS仪器开展多维探测时仍然需要有效的噪声消除方法.
通常，MRS数据中主要包含三类噪声：尖峰噪声、工频谐波噪声和随机噪声(Larsen and Behroozmand, 2016).尖峰噪声是由云层放电、汽车打火等引起的瞬时大幅度干扰，严重影响后续的滤波和叠加处理，因此应首先将其去除.目前，比较成熟的方法有能量运算方法(万玲等，2016)、尖峰建模方法(Larsen，2016)、统计叠加方法(Jiang et al., 2011)和小波模极大值方法(孙淑琴等，2015).工频谐波噪声是由高压输电线、变压器和其他电器设备产生的，在测量数据中最常出现，影响最大.由于工频谐波噪声在多通道采集时具有较强的相关性，因此可以采用维纳滤波(Dalgaard et al., 2012)或自适应噪声抵消(Müller-Petke and Costabel, 2014；田宝凤等，2012)等方法进行消除.但是，这些方法的消噪效果取决于多通道采集噪声的相关程度.对于随机噪声，通常采用叠加的方法进行减小，叠加次数越多，消噪效果越明显，但是消耗大量的测量时间.
除了降低噪声，从测量数据中直接提取有效MRS信号也是目前研究的热点.测量数据经过Hilbert变换后得到MRS信号的复包络曲线(Müller-Petke et al., 2016)，然后运用非线性拟合方法进行参数提取(Legchenko and Valla, 1998).此外，最新的MRS信号检测方法还包括独立成分分析(ICA)、经验模态分解(EMD)和奇异谱分析(SSA)(田宝凤等，2015；Ghanati et al., 2014, 2016)，但这些方法都是基于盲源分析理论，因此提取的MRS信号均存在部分失真的问题.
本文针对MRS仪器采集数据中的工频谐波噪声和随机噪声，鉴于MRS信号的表达式和基本特征是已知的，可以利用其自身在不同时间点的相关性(即自相关)进行有效检测；同时，工频谐波的特征也是已知的，利用这些特征进行建模，因而提出工频谐波噪声建模和自相关相结合的MRS信号提取方法和数据处理流程，摆脱对参考通道采集相关噪声的依赖性，同时解决以往只能通过多次叠加降低随机噪声的测量效率问题.首先，建立工频谐波噪声模型，改进基频搜索算法，自适应达到有效建模精度.其次，推导MRS信号自相关表达式，讨论其与测量数据自相关序列的关系.然后，建立MRS信号自相关参数提取的非线性拟合方程.通过仿真MRS信号测量数据，验证谐波建模和自相关结合方法对数据信噪比的提升效果，并分析MRS信号自相关参数提取的稳定性.最后，在长春市太平池水库，对已知地质信息的测量地点进行MRS实验，经过新的数据处理流程验证本文提出的方法对后续反演算法的准确性和稳定性具有重要的意义.
1 磁共振信号消噪方法 1.1 消噪处理流程根据MRS探测原理(Behroozmand et al., 2015)，地下水中的氢质子受到激发，其自旋发生偏转，发生核磁共振现象.激发停止后，在地面上的线圈接收到自由感应衰减信号，即MRS信号
|
(1) |
其中，V0称为MRS信号的初始振幅，与地下水中的氢质子含量呈正比；T2*称为平均弛豫时间(磁场均匀时与横向弛豫时间T2相等)，与赋水介质的孔隙大小相关；fL称为拉摩尔(Larmor)频率；φ0称为初始相位，与地下介质的电阻率、发射和接收频率偏差以及线圈的形状和位置相关.
MRS信号非常微弱，在1~1000 nV范围内，而外界环境的噪声和干扰严重，最大可达1×105 nV，因此接收数据的信噪比低，MRS信号提取的难度大.按照噪声的特征可以将其分成三类：尖峰噪声Vspike、谐波噪声Vharmonic和随机噪声Vrandom.此时，接收数据表示为
|
(2) |
为了从接收数据中准确提取MRS信号，常规的数据处理流程见图 1所示：1)去除尖峰噪声；2)去除工频谐波噪声；3)Hilbert变换和低通滤波；4)叠加和参数提取.参数提取的结果为MRS信号的四个关键参数(V0、T2*、φ0和df=fL-fT，其中fT是发射频率)，用于计算地下含水层的水文信息.
|
图 1 常规磁共振信号数据处理流程 Fig. 1 Conventional flow of data processing for MRS signal |
常规的消噪处理流程需要多通道仪器同步采集参考线圈的噪声数据，才能进行工频谐波噪声抵消，而且在信噪比较低时，需要多次重复测量数据进行叠加减小随机噪声.本文将用谐波建模消噪方法代替消噪流程的第2步，而在第3步和第4步之间加入自相关处理，降低随机噪声.
1.2 谐波建模原理谐波噪声Vharmonic是以工频(50±0.1 Hz)为基频f0的n个谐波的累加结果：
|
(3) |
其中，An和φn分别是第n个谐波的幅度和相位；N是谐波个数，根据接收信号的带宽(1~3 kHz)，取N=100.如图 2所示为一组接收信号数据，图 2a是其时域信号，MRS信号淹没在噪声之中；图 2c是其功率谱，在谐波的频点上均存在不同幅度的谐波噪声.
|
图 2 谐波建模原理和消噪效果 (a)信号时域图；(b)消噪后信号时域图；(c)信号频域图；(d)消噪后信号频域图. Fig. 2 Principle of powerline harmonic modeling and the effectcomparison between before and after de-noising (a) Signal in time domain; (b) Signal after de-noising in time domain; (c) Spectrum of signal; (d) Spectrum of signal after de-noising. |
谐波建模的原理是根据测量数据VR估计基频f0和每个谐波的An和φn，然后建立谐波噪声Vharmonic的模型数据.求解谐波参数(An，φn和f0)是一个非线性的优化问题，首先假设f0已知，将表达式(3)中的第n项线性化为
|
(4) |
其中，
|
(5) |
其中，Vp(p=1, 2, …, P)是tp时刻的接收数据.表达式(5)整理成矩阵形式为
|
(6) |
其中 x =[α1, …, αN, β1, …, βN]T，b=[V1, V2, …, Vp]T.表达式(6)是标准的线性方程组，可以采用最小二乘方法，无约束优化方法和SVD(奇异值分解)方法等求解(Lay，2012；Aster et al., 2005).
谐波建模后，在测量数据VR中减去Vharmonic即完成了消噪过程，如图 2b和2d所示.谐波建模消噪后，MRS信号的轮廓呈现衰减趋势，在功率谱中只存在MRS信号的频点.但是在实际中，工频的基频f0不能保证是精确的50 Hz，而是在50±0.1 Hz范围内变化，比如f0=50.01 Hz，其50次谐波为2500.5 Hz.如果只采用f0=50 Hz建模，谐波误差随着阶次的增加逐渐增大，其建模结果也会存在较大误差，因此谐波f0的数值必须进行精确估计.
1.3 自适应搜索谐波基频方法以上的计算过程都是基于f0已知的情况.当f0未知时，求解表达式(4)仍然是一个非线性问题.为了避免求解非线性优化问题，本文采用扫描搜索f0的方案，分为均匀扫描和自适应扫描两种方法.均匀扫描方法是在49.9~50.1 Hz的频率范围内，f0以0.001 Hz步进量通过表达式(6)分别进行建模计算，获得谐波的估计值Vharmonic(fl)(fl为第l个扫描值)，然后寻找与测量数据VR差值的2阶范数的最小值，即
|
(7) |
均匀扫描方式测量数据与建模数据的2阶范数随扫描频率fl的变化曲线见图 3a所示.扫描曲线在50.02 Hz附近出现凹陷，当fl=50.018 Hz时达到最小值，说明谐波的基频为50.018 Hz.采用均匀扫描方式，0.2 Hz的带宽需要重复计算200次，且只能精确到1 mHz，如果精度要求更高需要的计算次数更多，将消耗大量的时间.
|
图 3 工频基频搜索方法 (a)均匀扫描方法；(b)自适应扫描方法. Fig. 3 Two search methods for the base frequency of powerline harmonic noise (a) Uniform scanning method; (b) Adaptive scanning method. |
自适应扫描方式首先在49.9~50.1 Hz的范围内以0.03 Hz为步长进行粗扫，利用表达式(7)获得2阶范数最小的3个频率点作为下一次扫描的范围，见图 3b中的方块所示.在第2组扫描中，扫描范围为49.99~50.05 Hz，扫描步长为0.0075 Hz，再次获得范数最小的3个频点.当第3组扫描时，扫面范围缩小为50.0125~50.0275 Hz，扫描步长为0.001875 Hz.此时2阶范数的最小值对应的频点为50.018125 Hz，即经过3组扫描即可达到精度要求，所需要的计算次数和时间均小于均匀扫描方式.
2 自相关信号提取方法 2.1 自相关表达式自相关函数描述信号一个时刻与另一个时刻的依赖关系，即研究不同时刻两个信号变量的相关性.表达式(2)的仪器测量数据VR中，MRS信号VMRS和谐波噪声Vharmonic在不同时刻具有相关性，而尖峰噪声Vspike和随机噪声Vrandom不具备相关性.经过去尖峰和谐波建模消噪以后，测量数据中只剩下MRS信号VMRS和随机噪声Vrandom，因此可以利用自相关函数来压制随机噪声，从而进一步提高MRS信号的信噪比.
MRS信号VMRS的自相关函数表达式定义如下：
|
(8) |
式中τ是两个时刻的时间间隔.由于自相关函数是关于自变量τ的实偶函数，即R(τ)=R(-τ)，τ只取其正半轴.表达式(8)中的积分上下限分别取τ和tmax，tmax是测量时间长度，将MRS信号表达式(1)代入表达式(8)得到
|
(9) |
根据三角公式推导以及积分公式，将表达式(9)化简成
|
(10) |
其中，
|
(11) |
由表达式(10)和(11)可知，MRS信号自相关表达式包含四个分量，其中
|
(12) |
为了进一步化简MRS信号的自相关表达式(10)，将其四个分量表达式的计算结果分别与MRS离散数据VMRS(n)自相关序列的数值计算结果进行对比分析，见图 4所示.图 4中数据自相关序列的数值计算结果根据自相关定义计算：
|
(13) |
|
图 4 MRS信号自相关各分量表达式(RC)与数据自相关序列(AC)对比图 (a)、(d)、(g)和(j)RC各分量与AC对比；(b)、(e)、(h)和(k)局部放大图；(c) RC与AC的相对误差；(e)、(h)和(k) RC各分量占AC的比例. Fig. 4 Comparison of components of MRS signal autocorrelation expression (RC) and data autocorrelation sequence (AC) (a), (d), (g) and (j) components of RC and AC. (b), (e), (h) and (k) partial enlarged view; (c) The relative error of RC and AC. (e), (h) and (k) the proportion of componentsof RC to AC. |
根据表达式(10)计算得到的MRS信号的自相关结果RMRS与AC基本一致，见图 4a所示.图 4b是MRS信号的自相关结果的尾部(1~2 s)细节放大图，图 4c是RMRS对AC的相对误差曲线，可见两者的相对误差在0~2 s内均小于2%.
假设只利用自相关表达式的第一个分量RC1作为数据自相关序列的近似，由于不考虑tmax的影响，RC1与AC存在一定的偏差，主要出现在1 s以后的尾部，见图 4d和尾部细节放大图 4e所示.图 4f是RC1与AC的比值，即在0~1 s范围内，两者的比值接近1；当时间大于1 s后，两者的比值逐渐增大，由于AC(tmax=2s)=0，此时两者的比值为正无穷.因此，RC1只能作为数据自相关序列前半段结果的近似.同理，图 4g、4h和4i分别描绘了RC2与AC的关系，可见RC2在0~1 s内数值很小，可以忽略不计；在1~2 s内数值逐渐增大，与AC的比值也逐渐增大，最终趋近于负无穷.最后，图 4j、4k和4i分别描绘了RC3+RC4与AC的关系，可见RC3+RC4在0~2 s范围内的数值均较小，只占AC中0.017%左右，可以忽略.综上所述，当采集时间tmax足够长，数据自相关序列的前半段可以用表达式(12)中RC1近似.如果包含数据自相关序列的后半段，则可以用表达式(12)中RC1+RC2近似.
MRS信号经过自相关处理后仍然呈现衰减趋势，而随机噪声Vrandom由于不具备相关性，经过自相关处理后主要集中在0时刻，见图 5所示.图 5a中黑色线是理想的MRS信号，参数为V0=100 nV，T2*=400 ms，f=2326 Hz和φ0=π/4 rad；灰色线是加入了100 nV的高斯白噪声后的测量信号，此时数据的信噪比为0 dB，MRS信号基本淹没在噪声中.图 5b是测量信号和MRS信号的频谱，可见MRS信号的频谱峰值是噪声频谱均值的10倍.经过自相关函数后，测量信号(灰色线)和MRS信号(黑色线)的结果见图 5b所示.此时，灰色线与黑色线类似，也呈现明显的衰减趋势，只出现了少量的随机噪声，信噪比为23.75 dB，与自相关之前相比大幅提高.从图 5d中的频谱分析也可以得出，MRS信号自相关后的频谱峰值约为噪声频谱均值的100倍，与自相关之前图 5c相比提高了一个数量级.因此，自相关函数能够有效压制随机噪声，提高信噪比.
|
图 5 自相关处理减小随机噪声的效果 (a) MRS信号和加噪信号；(b)自相关后的MRS信号和加噪信号；(c) MRS信号和加噪信号的频谱；(d)自相关后MRS信号和加噪信号的频谱. Fig. 5 Effect of random noise reduction based on the autocorrelation method (a) MRS signal and noisy data; (b) MRS signal and noisy data after autocorrelation; (c) Spectrum of MRS signal and noisy data; (d) Spectrum of MRS signal and noisy data after autocorrelation. |
由2.1节可知，当采集时间tmax足够长时，MRS数据自相关序列的前半段可以用表达式(12)中RC1近似.为了提取MRS信号的关键参数V0，T2*和df(忽略相位φ0)，根据图 1中的数据处理流程，对RC1进行Hilbert变换后，再通过低通滤波，转化为两个分量VRC=Vreal+i·Vimag.最后，利用非线性拟合方法求解：
|
(14) |
式中
|
(15) |
其中
图 5中的仿真数据经过Hilbert变换后的两个分量见图 6所示，其中图 6a是未经过自相关处理的曲线，图 6b是经过自相关处理后的曲线.利用非线性拟合方法求解得到的两个分量的拟合曲线见图 6中黑色实线和黑色虚线所示.未经过自相关处理的数据信噪比低，参数提取结果分别为V0=108.3 nV，T2*=364.2 ms，f=2326.1 Hz；而经过自相关处理的数据信噪比提高明显，参数提取结果为V0=100.5 nV，T2*=404.1 ms，f=2326.0 Hz，准确度分别提高了16.6倍、8.7和10倍.因此，自相关处理能够提高数据的信噪比，从而提高参数提取的准确度.
|
图 6 自相关处理前后非线性拟合结果 (a)加噪信号(实部虚部)和拟合结果；(b)自相关处理后的信号(实部虚部)和拟合结果. Fig. 6 Results of nonlinear fitting before and after autocorrelation (a) Real and imagery party of noisy signal and the corresponding fitted results; (b) Real and imagery party of noisy signal after autocorrelation and the corresponding fitted results. |
为了验证谐波建模和自相关的消噪和信号提取效果，根据野外实验中仪器实际的测量过程，本文仿真了一组包含16次独立采集的测量数据集，见图 7a所示.仿真MRS信号的参数为V0= 100 nV，T2*= 500 ms，f=2326 Hz和φ0=π/4 rad；每次采集的谐波噪声基频在49.9~50.1 Hz随机产生，谐波个数为100，谐波幅度在200 nV内随机分布；每次采集的随机噪声为200 nV的高斯白噪声；采集数据中未加入尖峰噪声，假设所有的尖峰噪声已通过常规的方法消除.图 7a中颜色由深到浅的曲线分别代表 1~16次采集数据，采集时间为2 s，可见测量数据在±2000 nV之间变化，噪声较大，无法观测到明显的MRS信号衰减趋势，信噪比在-17.91~-13.34 dB之间.
|
图 7 仿真数据和数据处理结果 (a) 16组仿真加噪信号；(b)谐波建模消噪后；(c)自相关处理后；(d)和(e)自相关处理前和后的Hilbert变换结果；(f)和(g)自相关处理前和后的非线性拟合结果. Fig. 7 Synthetic data and the results after data processing (a) 16 sets of synthetic noisy data; (b) Denoising results of harmonic modelling; (c) Autocorrelation results; (d) and (e) Results of Hilbert transformation before and after autocorrelation; (f) and (g) Results of nonlinear fitting before and after autocorrelation. |
根据图 1给出的MRS信号处理的一般流程，对16次采集的数据分别进行处理.首先进行谐波建模消噪方法处理，利用第2.3小节给出的自适应基频频点搜索算法，在49.9~50.1 Hz范围内搜索谐波噪声的基频频率.然后利用表达式(4)进行谐波建模，在采集数据中将谐波建模结果减掉，消噪结果见图 7b所示.谐波建模消噪后，数据的信噪比大幅提升，平均提高了17.03 dB.数据幅度减小到±500 nV范围内，但是由于存在大量的随机噪声，仍然不能观察到衰减的MRS信号，数据信噪比在0.71~0.99 dB之间，需要进一步提高信噪比.
为了压制建模消噪后剩余的随机噪声，利用第2.1小节中表达式(13)对数据进行自相关处理，结果见图 7c所示.自相关处理后，数据呈现明显的衰减趋势，除了0时刻存在脉冲响应(随机噪声的自相关特性)，其他时刻的数据幅度比自相关之前的数据明显提高，由于随机噪声的影响明显减小，因此信噪比大幅度增加到15.64~16.43 dB，平均提高了16.10 dB.
按照数据处理流程，分别对建模消噪后和自相关后的数据进行Hilbert变换(参考频率2325 Hz)和低通滤波(截止频率200 Hz)，得到MRS信号的实部和虚部见图 7d和7e所示.建模消噪处理后，由于还存在大量的随机噪声，数据的实部和虚部曲线区分不明显，只能观察到曲线的整体趋势稍有起伏(1 Hz的振荡).经过自相关进一步处理后，数据的实部和虚部能够清晰地区分开，两者都呈现出1 Hz的振荡衰减过程.最后，对16次采集和处理的数据进行叠加，得到MRS信号实部和虚部的叠加曲线见图 7f和7g所示.叠加处理能够进一步降低随机噪声，16次叠加信噪比提高4倍(12.04 dB)，此时建模消噪后和自相关后的数据均呈现出振荡衰减趋势，但是自相关后的结果明显优于建模消噪后的结果.利用第2.2小节的非线性拟合方法对两组数据进行MRS信号参数提取，得到建模消噪后的参数提取结果：V0=98.27 nV，T2*= 470.2 ms，f=2325.9 Hz；而自相关处理后的参数提取结果为：V0=100.45 nV，T2*= 489.3 ms，f=2326.0 Hz，参数估计的准确度分别提高了3.8倍，2.8倍和10倍.因此，仪器采集数据经过谐波建模消噪后能够有效消除谐波噪声，提升信噪比，而自相关处理后，数据信噪比得以进一步改善，最终提高了MRS信号参数的估计精度.
3.3 信号提取的稳定性分析为了分析噪声水平和MRS信号参数变化时，自相关参数提取结果的稳定性，本文分别建立了4组噪声水平(50 nV、100 nV、200 nV和500 nV)和4组弛豫时间T2*(100 ms、200 ms、300 ms和400 ms)的数据，其他信号参数为e0=100 nV，f=2326 Hz和φ0=π/4 rad.每组数据重复进行100次计算，每次计算随机产生不同的白噪声，假设数据中已经消除了尖峰噪声和谐波噪声，经过自相关处理和Hilbert变换后，再进行非线性拟合，得到V0和T2*提取结果的统计见图 8a和8b所示.图中红线代表统计结果的中位数，外框代表参数提取结果出现的几率在25%~75%范围之内，黑色虚线代表结果的最大值和最小值的范围，“+”号代表明显存在错误的提取结果.
|
图 8 100组仿真数据的MRS信号参数提取结果统计 (a)不同噪声水平的e0统计；(b)不同噪声水平的T2*统计；(c)不同弛豫时间的e0统计；(d)不同弛豫时间的T2*统计. Fig. 8 Statistics of parameter extraction results from 100 groups of synthetic (a) Statistics of e0 with different noise levels; (b) Statistics of T2* with different noise levels; (c) Statistics of e0 with different relaxation times; (d) Statistics of T2* with different relaxation times. |
从图 8a和8b中的参数提取统计分析得到，虽然噪声变化，100次参数提取结果中e0和T2*的中位数均与设定参数(100 nV和500 ms)相近，但是两者的25%~75%范围框和最大最小范围随着噪声的增加而增大.当噪声为50 nV时，e0和T2*提取结果50%的几率相对误差小于2.36%和2.55%，最大相对误差仅为7.20%和9.67%，只有1次出现提取错误.当噪声为200 nV时，e0和T2*提取结果50%的几率相对误差小于8.68%和11.60%，最大相对误差分别为26.77%和36.01%，出现了8次提取错误.而当噪声为500 nV时，50%几率相对误差小于27.03%和43.25%，最大相对误差均超过了90%.因此，随噪声水平增加，MRS信号自相关参数提取结果的稳定性逐渐变差，当噪声大于200 nV以后，参数提取存在较大几率出现偏差较大的结果.
当弛豫时间T2*不同，而噪声水平相同时(200 nV)，MRS信号自相关处理后的参数提取结果见图 8c和8d所示.从图中的参数提取统计分析得到，100次参数提取结果中e0的25%~75%范围框和最大最小范围随着T2*的增加而减小，T2*的提取结果范围基本保持不变.由于环境噪声保持不变(200 nV)，T2*增大，MRS信号自相关幅度变大且衰减变慢，使得数据信噪比增加，因此对于较大T2*，自相关参数提取的效果更好.当T2*为200 ms时，e0和T2*提取结果50%的几率相对误差小于25.06%和39.47%，最大相对误差分别为60.03%和114.23%，结果明显比T2*=500 ms时(图 8a和8b)的准确度低.
4 实测数据验证为了进一步验证谐波建模消噪和自相关提取信号参数方法的应用效果，在长春市农安县烧锅镇的太平池水库旁，利用自主研制的MRS地下水探测仪器进行了野外数据采集与数据处理实验.当地的地磁场为54720 nT，对应的Larmor频率在2330 Hz，变化范围小于2 Hz.为了从浅到深地对含水层进行扫描，在0.2~8.5 As的范围内由小到大按对数分布设置了20组发射脉冲矩，每组脉冲矩重复叠加16次，采样率为50 kHz，采集时间为1 s.根据MRS信号处理流程，对每次采集数据分别进行了剔除尖峰噪声、谐波建模消噪、Hilbert变换和低通滤波、自相关、叠加和非线性拟合等处理，其中第10、16和20个脉冲矩的部分处理结果见图 9所示.
|
图 9 实测数据第10、16和20个脉冲矩的数据处理结果 (a)实测信号频谱和谐波建模消噪后频谱；(b)自相关处理前的信号(实部和虚部)和拟合结果；(c)自相关处理后的信号(实部和虚部)和拟合结果. Fig. 9 Measured data and the results after data processing at the 10th、16th and 20th pulse moments (a) Spectrum of measured data and processed after harmonic modelling; (b) Real and imagery part of signal and fitted results before autocorrelation; (c) Real and imagery part of signal and fitted results after autocorrelation. |
实测数据经过谐波建模消噪处理前后的频谱见图 9a所示，其中灰色曲线为实测数据的频谱，黑色曲线为谐波建模消噪后的频谱.由图可见，实测数据中存在大量的谐波噪声，噪声水平在2500~3000 nV范围内，而MRS信号幅度在80~375 nV，因此信噪比较低.经过谐波建模消噪处理后，大多数谐波噪声显著降低，谐波基频在49.975~50.041 Hz之间，MRS信号(2330 Hz)的峰值比较明显，其中q=1.4 As的信号较大，而q=4.2 As和8.6 As的信号较小.虽然在部分谐波频点处仍然能观察到峰值，但是对后续的MRS信号参数提取结果影响不大，总体噪声水平降至860~950 nV，即信噪比提高了近20 dB.较高频率的谐波噪声滤波效果一般，是因为谐波基频随时间存在微小变化，引起较高阶次的谐波频率估计偏差变大，因此消噪效果变差.
谐波消噪处理后的数据直接经过Hilbert变换和非线性拟合后的结果见图 9b所示，而加入自相关处理后再经过Hilbert变换和非线性拟合后的结果见图 9c所示，其中浅颜色的两条线分别是Hilbert变换后的实部和虚部曲线，深颜色的两条线分别是非线性拟合后的实部和虚部曲线.由图 9b可以看出，经过多种消噪处理后的数据信噪比仍然较低，参数提取结果的不确定度较大，见表 1所示.实测信号经过自相关处理后，信噪比明显提高，如图 9c和表 1所示，3个脉冲矩下信噪比增量分别为12.05 dB、2.98 dB和7.71 dB，此时经过非线性拟合后，e0和T2*结果的不确定度均得到改善.其中q=4.2 As时，MRS信号幅度最小，拟合结果的e0不确定度由36.1 nV减小到15.9 nV，而T2*拟合结果的不确定度由117.9 ms减小到38.2 ms，说明增加SNR，能够提高非线性拟合结果的准确度.
|
表 1 实测数据处理结果和参数提取结果 Table 1 Results of data processing and parameter extraction from the measured data |
谐波建模消噪和自相关参数提取方法能够显著提高实测数据的信噪比以及非线性拟合的准确度，从而可以实现对地下水文信息精确解释，见图 10所示.采集的实测数据自相关处理前后的参数提取结果见图 10a所示.20个脉冲矩下未经过自相关处理前的不确定度(误差棒)均较大，反演曲线(黑色)与测量曲线(浅色)偏差较大；而经过自相关处理后不确定度(误差棒)较小，反演曲线(黑色)与测量曲线(浅色)的偏差较小.图 10b是利用不同叠加次数测量曲线经过反演得到的地下含水量信息，较少的叠加次数对应较大的噪声水平.图 10c是当地的钻探资料，用于验证反演结果.对于16次叠加的测量曲线，未经过自相关处理的反演结果与钻探资料基本相符，但是在5 m处出现异常含水层，14~20 m处的含水层含水量偏低，而40~75 m的含水层与细砂层和砾质砂层相对应，但是分辨率较低，含水层边界模糊，含水量偏低.对于8次叠加和4次叠加的测量曲线，未经自相关处理的反演结果噪声较大，在10 m内出现多个异常含水层，而且主要含水层和深部含水层均存在偏差.相比之下，自相关处理后的反演结果与钻探资料相符，不同叠加次数的测量曲线得到的反演结果基本一致，其中10 m以内不存在异常的含水层，14~22 m的含水层与砾质砂层相对应，含水量为16.6%，而45~65 m的含水层与细砂层和砾质砂层相对应，含水量最大的位置出现在58~60 m的砾质砂层中.但是，随着叠加次数减小，噪声水平增加，对深部含水层的分辨率也逐渐下降。因此，通过对比可以证明，自相关处理能够明显提高测量曲线的质量，从而改善反演结果的准确度和稳定性.
|
图 10 长春太平池水库实测数据反演结果和钻探结果对比 (a)自相关处理前和处理后的MRS信号初始振幅和反演结果计算数据；(b)不同叠加次数自相关处理前和处理后的含水量反演结果；(c)钻探结果. Fig. 10 Inversion results and the borehole logging in Taipingchi reservoir, Changchun (a) Initial amplitudes of MRS signal and calculated data based on the inversion result before and after autocorrelation; (b) Inversion results before and after autocorrelation with different stacks; (c) Results of borehole logging. |
针对MRS方法存在采集数据信噪比低和信号检测难度大的问题，本文提出了基于谐波建模和自相关参数提取相结合的MRS数据处理方法和流程，以消除或减小影响MRS信号质量最大的两类噪声：工频谐波噪声和随机噪声.常规数据处理流程中采用参考对消的方式消除工频谐波噪声，需要多通道仪器同步采集噪声数据；采用多次测量数据后叠加的方式减小随机噪声，消耗大量的时间且效果不明显.本文提出的方法避免了多通道仪器和多次测量的限制，解决了在城镇和农村附近强噪声干扰下获取有效MRS信号的难题.
本文首先建立了工频谐波的模型表达式，推导了谐波系数的求解方程，针对建模算法严重依赖工频基频精确度的问题，提出了自适应搜索工频基频的方案，相比均匀扫描方法大幅提高搜索的准确度和速度.通过仿真数据消噪实验证明，建模消噪后多数的工频谐波都被有效消除，数据信噪比平均提升了17.03 dB.与常规参考对消方法相比，谐波建模方法不需要额外的采集通道，也不需要单独采集噪声数据，只需采集的信号数据即可建立工频谐波模型，从而消除工频谐波噪声.
其次，本文提出了利用自相关处理压制随机噪声并进行参数提取的方法，推导了MRS信号的自相关表达式，通过仿真证明了MRS数据自相关结果序列的前半段可以用RC1代替，而后半段必须用RC1+RC2代替.同时，实现了针对MRS自相关参数提取的非线性拟合方法.经过仿真数据消噪实验证明，自相关处理后，数据信噪比再次提升16.10 dB，初始振幅和弛豫时间的参数提取结果的准确度比自相关处理前分别提高了3.8倍和2.8倍.通过大量不同信噪比，不同弛豫时间的仿真数据参数提取实验对比，得到在噪声水平小于200 nV时，自相关参数提取方法具有较高的稳定性，随着噪声增加和弛豫时间变短后，稳定性逐渐变差.
最后，通过在长春市太平池水库旁的野外数据采集与数据处理实验，进一步验证了谐波建模消噪和自相关参数提取方法的应用效果：经过谐波建模消噪后，采集数据的信噪比提高了20 dB，再经过自相关处理后信噪比提高了2.98 dB以上，且参数提取的不确定度均大幅降低.因此，两种方法相结合可以显著提高实测信号的信噪比和参数提取的准确度.利用实测数据进行地下含水量反演，并与钻探资料对比，结果进一步表明本文所提方法能够提高测量曲线的质量，使得反演结果的准确性和稳定性获得明显改善.
Aster R C, Brochers B, Thurber C H. 2005.
Parameter Estimation and Inverse Problems. Burlington: Elsevier Academic Press.
|
|
Behroozmand A A, Keating K, Auken E.
2015. A review of the principles and applications of the NMR technique for near-surface characterization. Surveys in Geophysics, 36(1): 27-85.
|
|
Dalgaard E, Auken E, Larsen J J.
2012. Adaptive noise cancelling of multichannel magnetic resonance sounding signals. Geophys. J. Int., 191(1): 88-100.
DOI:10.1111/gji.2012.191.issue-1 |
|
Ghanati R, Fallahsafari M, Hafizi M K.
2014. Joint application of a statistical optimization process and empirical mode decomposition to magnetic resonance sounding noise cancelation. Journal of Applied Geophysics, 111: 110-120.
DOI:10.1016/j.jappgeo.2014.09.023 |
|
Ghanati R, Hafizi M K, Mahmoudvand R, et al.
2016. Filtering and parameter estimation of surface-NMR data using singular spectrum analysis. Journal of Applied Geophysics, 130: 118-130.
DOI:10.1016/j.jappgeo.2016.04.005 |
|
Grunewald E, Knight R, Walsh D.
2014. Advancement and validation of surface NMR spin echo measurements of T2. Geophysics, 79(2): EN15-EN23.
DOI:10.1190/geo2013-0105.1 |
|
Grunewald E, Grombacher D, Walsh D.
2016. Adiabatic pulses enhance surface nuclear magnetic resonance measurement and survey speed for groundwater investigations. Geophysics, 81(4): WB85-WB96.
DOI:10.1190/geo2015-0527.1 |
|
Hertrich M.
2008. Imaging of groundwater with nuclear magnetic resonance. Progress in Nuclear Magnetic Resonance Spectroscopy, 53(4): 227-248.
DOI:10.1016/j.pnmrs.2008.01.002 |
|
Jiang C D, Lin J, Duan Q M, et al.
2011. Statistical stacking and adaptive notch filter to remove high-level electromagnetic noise from MRS measurements. Near Surface Geophysics, 9(5): 459-468.
|
|
Jiang C D, Liu J Y, Tian B F, et al.
2016. 3D block QT inversion of surface nuclear magnetic resonance data. Geophysics, 81(5): E363-E376.
DOI:10.1190/geo2015-0582.1 |
|
Jiang C D, Shang X L, Lin T T, et al.
2017. Quasi-2D block inversion of large-scale surface nuclear magnetic resonance profile data using a laterally constrained model. Geophysics, 82(2): JM1-JM14.
DOI:10.1190/geo2015-0455.1 |
|
Larsen J J, Dalgaard E, Auken E.
2014. Noise cancelling of MRS signals combining model-based removal of powerline harmonics and multichannel Wiener filtering. Geophysical Journal International, 196(2): 828-836.
DOI:10.1093/gji/ggt422 |
|
Larsen J J, Behroozmand A A.
2016. Processing of surface-nuclear magnetic resonance data from sites with high noise levels. Geophysics, 81(4): WB75-WB83.
DOI:10.1190/geo2015-0441.1 |
|
Larsen J J.
2016. Model-based subtraction of spikes from surface nuclear magnetic resonance data. Geophysics, 81(4): WB1-WB8.
DOI:10.1190/geo2015-0442.1 |
|
Lay D C. 2012.
Linear Algebra and its Applications. (4th ed). New York: Pearson Education.
|
|
Legchenko A, Valla P.
1998. Processing of surface proton magnetic resonance signals using non-linear fitting. Journal of Applied Geophysics, 39(2): 77-83.
DOI:10.1016/S0926-9851(98)00011-1 |
|
Legchenko A, Valla P.
2003. Removal of power-line harmonics from proton magnetic resonance mea-surements. Journal of Applied Geophysics, 53(2-3): 103-120.
DOI:10.1016/S0926-9851(03)00041-7 |
|
Lin J, Jiang C D, Lin T T.
2013. Underground magnetic resonance Sounding (UMRS) for detection of disastrous water in mining and tunneling. Chinese Journal of Geophysics, 56(11): 3619-3628.
DOI:10.6038/cjg20131103 |
|
Lin T T, Jiang C D, Qi X, et al.
2013. Theories and key technologies of distributed surface magnetic resonance sounding. Chinese Journal of Geophysics, 56(11): 3651-3662.
DOI:10.6038/cjg20131106 |
|
Müller-Petke M, Costabel S.
2014. Comparison and optimal parameter setting of reference-based harmonic noise cancellation in time and frequency domain for surface-NMR. Near Surface Geophysics, 12(2): 199-210.
|
|
Müller-Petke M, Braun M, Hertrich M, et al.
2016. MRSmatlab-A software tool for processing, modeling, and inversion of magnetic resonance sounding data. Geophysics, 81(4): WB9-WB21.
DOI:10.1190/geo2015-0461.1 |
|
Pan Y L, Zhang C D. 2000.
Theories and Methods of Surface Nuclear Magnetic Resonance. Beijing: China University of Geosciences Press.
|
|
Roy J, Lubczynski M.
2003. The magnetic resonance sounding technique and its use for groundwater investigations. Hydrogeology Journal, 11(4): 455-465.
DOI:10.1007/s10040-003-0254-8 |
|
Strehl S. 2006. Development of strategies for improved filtering and fitting of SNMR-signals[Ph. D. thesis]. Berlin:Technical University of Berlin.
|
|
Sun S Q, Meng Q Y, Fang X C, et al.
2015. Noise cancellation algorithm for surface magnetic resonance signals based on self-adaption and reconstruction of wavelet modulus maximum value. Journal of Jilin University (Engineering and Technology Edition), 45(5): 1642-1651.
|
|
Tian B F, Lin J, Duan Q M, et al.
2012. Variable step adaptive noise cancellation algorithm for magnetic resonance sounding signal with a reference coil. Chinese Journal of Geophysics, 55(7): 2462-2472.
DOI:10.6038/j.issn.0001-5733.2012.07.030 |
|
Trushkin D V, Shushakov O A, Legchenko A V.
1994. The potential of a noise-reducing antenna for surface NMR groundwater surveys in the earth's magnetic field. Geophysical Prospecting, 42(8): 855-862.
DOI:10.1111/gpr.1994.42.issue-8 |
|
Tian B F, Zhou Y Y, Wang Y, et al.
2015. Noise cancellation method for full-wave magnetic resonance sounding signal based on independent component analysis. Acta Physica Sinica, 64(22): 229301.
DOI:10.7498/aps.64.229301 |
|
Walsh D O.
2008. Multi-channel surface NMR instrumentation and software for 1D/2D groundwater investigations. Journal of Applied Geophysics, 66(3-4): 140-150.
DOI:10.1016/j.jappgeo.2008.03.006 |
|
Wan L, Zhang Y, Lin J, et al.
2016. Spikes removal of magnetic resonance sounding data based on energy calculation. Chinese Journal of Geophysics, 59(6): 2290-2301.
DOI:10.6038/cjg20160631 |
|
Zhang R, Hu X Y, Yang D K, et al.
2006. Review of development of surface nuclear magnetic resonance. Progress in Geophysics, 21(1): 284-289.
|
|
林君, 蒋川东, 林婷婷, 等.
2013. 地下工程灾害水源的磁共振探测研究. 地球物理学报, 56(11): 3619–3628.
|
|
林婷婷, 蒋川东, 齐鑫, 等.
2013. 地面磁共振测深分布式探测方法与关键技术. 地球物理学报, 56(11): 3651–3662.
|
|
潘玉玲, 张昌达. 2000.
地面核磁共振找水理论和方法. 北京: 中国地质大学出版社.
|
|
孙淑琴, 孟庆云, 方秀成, 等.
2015. 基于自适应和小波模极大值重构的地面核磁共振信号噪声压制方法. 吉林大学学报(工学版), 45(5): 1642–1651.
|
|
田宝凤, 林君, 段清明, 等.
2012. 基于参考线圈和变步长自适应的磁共振信号噪声压制方法. 地球物理学报, 55(7): 2462–2472.
|
|
田宝凤, 周媛媛, 王悦, 等.
2015. 基于独立成分分析的全波核磁共振信号噪声滤除方法研究. 物理学报, 64(22): 229301.
DOI:10.7498/aps.64.229301 |
|
万玲, 张扬, 林君, 等.
2016. 基于能量运算的磁共振信号尖峰噪声抑制方法. 地球物理学报, 59(6): 2290–2301.
DOI:10.6038/cjg20160631 |
|
张荣, 胡祥云, 杨迪琨, 等.
2006. 地面核磁共振技术发展述评. 地球物理学进展, 21(1): 284–289.
|
|