2018, Vol. 61
富水介质中氢质子在外磁场的作用下发生能级跃迁产生磁共振(Magnetic Resonance Sounding，简称MRS)信号(Weichman et al., 2000；Legchenko and Valla, 2002；Behroozmand et al., 2015)，该信号可表征氢质子丰度和水体赋存状态等信息(林君等，2011；潘玉玲和张昌达，2000；Legchenko et al., 2002), 磁共振地下水探测方法因具有对地下水非侵入性直接探测和定量评估的能力，近年来得到了快速的发展(林君，2010；林君等，2012；林君等，2016；陆其鹄等，2009；张荣等，2006).获取可靠的MRS信号是有效进行磁共振地下水探测的前提，但MRS信号极其微弱，仅为nV级.高灵敏度的探测系统在地磁场环境下开展工作，无法通过屏蔽等方式隔绝环境中的电磁噪声干扰，导致信噪比低，这一问题严重制约着磁共振地下水探测方法的应用和发展.
为提高MRS信号信噪比，国内外专家和学者开展了相关研究.Packard和Varian(1954)提出直流预极化场概念，通过将样品放置于螺线管或永磁体中以增强氢质子的宏观磁矩.De Pasquale和Mohnke(2014)通过模型仿真验证了施加预极化场可有效增强地下水的信号强度.Lin等(2018)于国际上首次使用预极化场磁共振方法在野外测得了冰层下流动的水信号.鉴于超导量子干涉器件(Superconducting Quantum Interference Device，简称SQUID)的灵敏度可达fT级，在检测微弱信号方面具有独特的优势，Davis等(2014, 2015)采用SQUID代替传统接收线圈提高对信号的接收灵敏度.另外，满足绝热条件时，通过增加扳倒角也能提高磁化强度的有效分量(Tannús and Garwood, 1997；Garwood and DelaBarre, 2001)，Grunewald等(2016)将绝热脉冲理念应用到了磁共振地下水探测中，得到了绝热激发的磁共振响应信号及反演解释结果.尽管这些方法在增加信号强度或提高仪器接收灵敏度方面取得了良好的效果，但在城市周边或村庄附近探测时仍不可避免的受到严重的电磁噪声干扰，需要根据信号和噪声的特点采取有效的数字滤波手段进一步提高信噪比，得到可靠的MRS信号.
干扰MRS信号的噪声主要分为工频及其谐波噪声、尖脉冲噪声和随机噪声.工频及其谐波噪声来源于电力线、发电机和变电器等(林君等，2011)，频率固定为50 Hz或60 Hz的整数倍，基波频率会有一定的偏差，但偏差通常小于0.1 Hz；尖脉冲噪声由太阳磁暴、雷暴或任何物体突然放电等产生，持续时间很短，幅度大于甚至远大于信号幅度，这些噪声都具有一定的统计规律，属于确定性噪声，众多研究成果表明工频及其谐波噪声和尖脉冲噪声可被有效的消除(Dalgaard et al., 2012；Müller-Petke and Costabel, 2014；Larsen et al., 2014；Costabel and Müller-Petke，2014；Larsen，2016；万玲等，2016).然而，随机噪声是由各种不可预知因素综合作用而成，没有统一的规律，在统计特征上具有随机性.由于随机噪声频带很宽，有效信号淹没其中难以辨识.针对磁共振探测随机噪声的抑制，通常采用叠加法(Legchenko and Valla, 2002；Dalgaard et al., 2012；Dalgaard et al., 2016)，将多次采集的信号叠加以抑制该类型的噪声干扰.但是，该方法只能抑制部分随机噪声，且多次采集信号会降低探测效率.除此之外，尚未发现有效的MRS信号随机噪声抑制方法.
时频峰值滤波具有良好的时频聚集性和较强的随机噪声消减能力，该方法最早由Boashash的博士生Arnold和Roessgen发明(Arnold et al., 1994；Arnold and Boashash, 1994；Roessgen and Boashash, 1994)，随后被用于新生儿脑电信号的增强(Boashash and Mesbah, 2004)，近年来，时频峰值滤波在地震信号处理领域也取得了良好的应用效果(李月等，2009；林红波等，2011；Wu et al., 2016).Lin等(2018)借鉴该方法的成功应用，将其用于抑制MRS包络信号的随机噪声.纵然该方法在MRS包络信号的随机噪声抑制方面表现出了良好的特性，但是，为更准确提取信号的有效信息，近年来多采用宽频带磁共振地下水探测仪收录MRS全波信号(Walsh，2008；林婷婷等，2014；曲永星，2015).由于MRS全波信号数据点数多，直接采用时频峰值滤波法对信号进行时频分析会导致计算机内存溢出.即便使用极高性能的计算机，但由于信号非线性强，通常不能得到信号的无偏估计，滤波结果与真实信号偏差很大.
针对上述问题，本文开展分段时频峰值滤波(Segment Time-Frequency Peak Filtering，简称STFPF)研究，分析随机噪声对MRS全波信号的影响，提出基于STFPF的全波磁共振信号的随机噪声抑制方法.首先，根据信号采样点数将其分成若干段，分别进行瞬时频率编码；其次，采用短窗长伪魏格纳维尔分布(Pseudo Wigner-Ville Distribution，简称PWVD)对解析信号进行时频分析，利用解析信号的时频分布沿其瞬时频率集中的特性，通过时频分布的峰值得到信号的无偏估计，该方法为提高特征参数提取精度和反演解释准确性奠定了基础.
1 随机噪声对全波磁共振信号影响分析磁共振地下水探测技术是一种能够实现对地下水定性定量探测的地球物理方法.利用“百米级”天线，通过发射拉莫尔频率脉冲场，测量得到随时间呈指数衰减的谐振信号为
|
(1) |
该信号与地层下含水量和水体的赋存状态有关.其中，信号初始振幅E0的大小与含水层的含水量成正比，信号的平均横向弛豫时间T2*与地下介质的孔隙度有关，信号的振荡频率ω0与探测地点的地磁场强度有关，信号的初始相位φ0与地下介质的电导率有关.早期为提高磁共振地下水探测仪的测量速度，主要获取MRS包络信号(IRIS Instruments，2004；王中兴等，2010).随着探测条件越来越复杂，噪声导致甚微弱磁共振信号难以被有效提取，迫切需要提供能够携带丰富信息的MRS全波信号，以便采用更先进的数据处理方法消除复杂噪声干扰.此后，伴随采集卡的更新换代，在满足测量速度的前提下，提高仪器的采样频率，增加采集数据量，磁共振地下水探测仪实现了MRS全波信号的收录(Walsh，2008；林婷婷等，2014；曲永星，2015).然而，较宽的通频带会引入更多的随机噪声，导致目前仪器采集的信号中随机噪声干扰严重，影响磁共振信号的检测及特征参数的提取.
随机噪声对MRS全波信号的影响如图 1所示.图 1a蓝色曲线表示纯净的MRS全波信号，它的初始振幅为200 nV，平均横向弛豫时间为150 ms，信号频率为2330 Hz，初始相位为1.05 rad；红色衰减曲线是该信号的拟合曲线，根据拟合结果获取信号的特征参数.由于MRS全波信号采样频率高，数据量大，为了进一步观察信号特征，选取100~105 ms区间段对信号进行放大(图 1c)，可见MRS全波信号是按固定频率衰减的正弦曲线.图 1b是被均值为0 nV、标准差为150 nV随机噪声干扰的MRS全波信号，从图中可以看出，随机噪声干扰后信号的拟合曲线与图 1a相比差别较大，信号的初始振幅和衰减趋势也有较大的偏差.图 1d是含噪信号在同样100~105 ms范围内的局部放大，可见正弦曲线在随机噪声干扰下波形发生畸变，这将导致信号的特征参数难以被准确获取，反演结果偏离实际.
|
图 1 纯净的MRS全波信号和随机噪声干扰下的MRS全波信号 (a)纯净的MRS全波信号及其拟合曲线；(b)随机噪声干扰下的MRS全波信号及其拟合曲线；(c)纯净的MRS全波信号在100~105 ms范围内的局部放大；(d)随机噪声干扰下的MRS全波信号在100~105 ms范围内的局部放大. Fig. 1 Pure MRS oscillating signal and random-noise-corrupted signal (a) Pure MRS oscillating signal and its fitting line; (b) Random-noise-corrupted signal and its fitting line; (c) Zoom in the pure signal at 100~105 ms; (d) Zoom in the random-noise-corrupted signal at 100~105 ms. |
由于全波磁共振信号采样频率高，数据点数多，直接运用时频峰值滤波法，时频分析矩阵维度大，计算机难以负荷.假设含随机噪声的全波磁共振信号可用式(2)的形式表示，其中m表示采样时刻，s(m)表示有效信号，即理想的MRS信号；n(m)表示随机噪声，它可能与有效信号具有重叠的频谱.本文提出基于分段时频峰值滤波的全波磁共振信号随机噪声抑制方法，从含噪声信号中恢复出有效信号s(m)，即
|
(2) |
基于分段时频峰值滤波的全波磁共振信号消噪流程见图 2，其过程可分为三个步骤：首先将全部数据分成若干段，分别进行时频峰值滤波处理；其次，将每一段信号编码为解析信号的瞬时频率，通过时频变换求取解析信号的时频分布，利用时频分布沿瞬时频率集中的特性，通过提取时频分布的峰值获得信号的无偏估计；最后将信号重组，实现对全波磁共振信号随机噪声的抑制处理.
|
图 2 全波MRS信号随机噪声抑制流程 Fig. 2 Flowchart of the random noise reduction of MRS oscillating signal |
假设MRS全波信号的采样点数为M，将信号分成I段，每段信号的长度为W，第i段信号表示为
|
(3) |
式中I为正整数，表示分段数.当I×W大于信号的采样点数M时，第I段信号由剩余的采样点数组成，即xI(m)=x(m), (I－1)W+1≤m≤M.
分别对每一段信号进行尺度变换，将信号的幅度值限定在归一化的频率范围之间，避免信号在进行频率调制时出现失真现象.xi(m), m=1, 2, …, W; i=1, 2, …, I，经过尺度变换后得到信号xic(m)，公式为
|
(4) |
式中xic(m)是尺度变换后的信号，系数a和b满足0.5≥a=max[xic(m)]，b=min[xic(m)]≥0，它们分别是变换后信号的最大值和最小值，可以通过选择a和b来确定编码信号的频率限制.
将尺度变换后的信号进行频率调制，编码为单位幅度解析信号的瞬时频率，得到的解析信号表示为
|
(5) |
使用魏格纳维尔分布(Wigner-Ville Distribution，简称WVD)计算由式(5)表示的解析信号的时频分布，可以得出
|
(6) |
再按频率变量求取时频分布Wzxic(m, k)的峰值作为瞬时频率的估计值，得到
|
(7) |
式中，argfmax[·]表示沿频率取最大值算子.
由式(7)得到的有效信号估计值的幅度仍在归一化频率范围内，需通过反尺度变换得到信号的实际幅度，可以得到
|
(8) |
最后，将分段滤波后的结果移位求和以得到完整的抑制随机噪声后的MRS全波信号
|
(9) |
将式(2)代入式(5)，可以得到
|
(10) |
式中，sic(m)和nic(m)分别是si(m)和ni(m)尺度变换后的有效信号和随机噪声；zsic(m)和znic(m)分别是以sic(m)和nic(m)为瞬时频率的解析信号.
同样，使用WVD得到的解析信号的时频分布也可以表示为
|
(11) |
再计算该解析信号的魏格纳维尔谱(Wigner-Ville-Spectrum，WVS)，可以得出
|
(12) |
其中，E[·]是求数学期望算子.由式(12)可知，瞬时频率的估计误差主要来源于信号zsic(m)的时频分布所引起的确定性误差和噪声znic(m)的WVS所引起的随机性误差.
由于
|
(13) |
可见，利用Wzx(t, f)的峰值估计瞬时频率时误差只来源于Wzs(t, f).此时，假设信号是时间的线性函数，即s(t)=αt+c，α和c是常数，得到
|
(14) |
因此，信号为线性是保证STFPF算法得到滤波结果无偏估计的前提条件.然而，MRS全波信号不是时间的线性函数，需要对信号进行局部线性化处理，采用加窗的WVD，即PWVD使信号在窗长内近似线性，计算解析信号的时频分布，表示为
|
(15) |
式中，w(l)是窗函数，宽度W为2L+1.为尽可能的保证信号在窗长内为线性，Boashash和Mesbah(2004)推导了相对误差小于0.2时，与信号频率及信号采样频率有关的窗长公式，得出
|
(16) |
式中，f和fs分别是信号频率和信号的采样频率.对于JLMRS-阵列式宽频带全波信号接收仪，其采样频率为5×104 Hz，信号频率与探测地点的拉莫尔频率相同(吉林省长春市附近为2330 Hz)，故根据式(16)，经计算，运用STFPF算法抑制MRS全波信号的随机噪声时，式(15)中窗长取7.
为验证理论推导，分别构造线性信号和非线性信号：将图 3a中蓝色实线所示的斜率-2、过点(0，1)的线性信号编码为单位幅度解析信号的瞬时频率，使用WVD计算解析信号的时频分布如图 3b所示.从图中可以看出，解析信号时频分布的能量沿着瞬时频率集中，而瞬时频率恰好是构造的线性信号，故沿频率取解析信号时频分布的峰值即可恢复信号，如图 3a中红色虚线所示.对比滤波前后的滤波结果可知，当信号是时间的线性函数时，使用WVD对信号做时频分布可获得信号的无偏估计.然而，对于频率为2 Hz、峰值为±1 V的非线性正弦信号(图 3c蓝色实线)，同样使用WVD对信号进行时频分析，如图 3d所示.从图中可以看出，解析信号时频分布能量仍沿着瞬时频率集中，但由于交叉干扰项的存在，按频率取解析信号时频分布的峰值恢复出的信号波形(图 3c红色虚线)与滤波前的信号波形在过零点处偏差较大，且波峰和波谷失真.可见，当信号为非线性时，通过WVD得到的瞬时频率估计值有偏.再采用PWVD对信号进行处理，图 3e和f给出了运用PWVD对信号进行时频分析得到的滤波结果.从图 3f中可以看出，解析信号时频分布的能量沿瞬时频率集中，此时没有交叉项，沿频率取解析信号时频分布的峰值即可恢复信号，如图 3e中红色虚线所示.从滤波前后信号的对比结果可以得出：对编码后的解析信号进行时频分析时，为信号加适当的时窗函数，使信号在时窗内尽可能线性，运用STFPF算法可获得非线性信号滤波结果的无偏估计.
|
图 3 滤波前后的信号及其解析信号的时频分布图 (a)采用基于WVD的时频峰值滤波法滤波前后的线性信号；(b)解析信号瞬时频率为线性信号时的时频分布；(c)采用基于WVD的时频峰值滤波法滤波前后的非线性信号；(d)解析信号瞬时频率为非线性信号时的时频分布；(e)采用基于PWVD的时频峰值滤波法滤波前后的非线性信号；(f)解析信号瞬时频率为非线性信号时的时频分布. Fig. 3 Ideal and estimated signals and the time-frequency distributions of the corresponding analytic signals (a) A linear signal before and after filtering; (b) WVD of the analytic signal whose IF is the linear signal; (c) a sine signal before and after filtering; (d) WVD of the analytic signal whose IF is the sine signal; (e) a sine signal before and after filtering; (f) PWVD of the analytic signal whose IF is the sine signal. |
当采集的MRS信号中存在随机噪声时，通常采用叠加法减小该噪声(Legchenko and Valla, 2002)，其基本思想是在某一激发脉冲矩下多次采集MRS信号，假设采集次数为N，对N次采集的信号相加求平均即为叠加后的MRS信号.采用叠加法抑制随机噪声后，信噪比可提高N1/2倍.但是，为了得到较好的消噪效果，需要增加叠加次数，导致探测效率随之降低.并且，数据叠加法对随机噪声的抑制效果有限，当叠加次数增加到一定程度后(通常为128次)，再继续增加叠加次数，信噪比几乎不再提高.
通过算法的数值仿真分别验证传统叠加方法和STFPF对MRS全波信号随机噪声的抑制效果，除了能从信号波形上观察随机噪声的消除情况，本文综合了信噪比、初始振幅拟合误差和平均横向弛豫时间拟合误差进行量化评估，信噪比定义为
|
(17) |
式中，信噪比SNR采用了dB的表示方式.
初始振幅拟合误差(Initial Amplitude Percentage Error，IAPE)定义为
|
(18) |
式中，ΔE0为初始振幅拟合值和真实值的偏差.
平均横向弛豫时间拟合误差(Relaxation Time Percentage Error，RTPE)定义为
|
(19) |
式中，ΔT2*为平均横向弛豫时间拟合值和真实值的偏差.
仿真信号模型参数如下：MRS全波信号的初始振幅E0=200 nV，平均横向弛豫时间T2*=150 ms，信号的中心频率fL=2128 Hz，相位φ0=-1.05 rad，加入平均值为0 nV、标准差为100 nV的随机噪声干扰，使得信噪比SNR=－5.17 dB，初始振幅拟合值为254.98 nV，初始振幅拟合误差IAPE=27.49%，平均横向弛豫时间拟合值为998 ms，平均横向弛豫时间拟合误差RTPE=565.3%，图 4a为具有25000个等精度采样点的时间域仿真含噪信号，图 4b为对应的频谱图.
|
图 4 仿真的含随机噪声的全波磁共振信号及其频谱 (a)时域信号图；(b)频谱图. Fig. 4 Random-noise-corrupted signal and the spectrum (a) Signal in time domain; (b) Signal in frequency domain. |
分别采用叠加法和STFPF方法对图 4所示的仿真含随机噪声信号进行消噪处理.图 5为叠加8次时，数据叠加法的消噪结果.图中蓝色线表示含随机噪声的信号，绿色线表示叠加后的信号.为了观察消噪前后信号波形，任意选取2~5 ms和112.6~115.2 ms对信号进行局部放大，见图 5b和c，可以明显看出，叠加后的信号仍含有部分随机噪声.运用叠加法后，信噪比为3.79 dB，提高了8.96 dB.另外，初始振幅E0=217.19 nV，初始振幅的拟合误差IAPE=8.6%，平均横向弛豫时间T2*为306 ms，平均横向弛豫时间的拟合误差为104%.
|
图 5 运用叠加法抑制随机噪声前后信号的时域图 (a)含随机噪声的MRS全波信号；(b) 2~5 ms范围内的局部放大信号；(c) 112.6~115.2 ms范围内的局部放大信号. Fig. 5 The signal in time domain before and after random noise attenuation by traditional stacking method (a) The noisy MRS oscillating signal; (b) MRS oscillating signal at 2~5 ms; (c) MRS oscillating signal at 112.6~115.2 ms. |
再运用STFPF方法抑制该信号中的随机噪声，将全波信号均匀分为100段，每一段包含250个等精度采样点，分别编码为单位幅度解析信号的瞬时频率，根据式(15)和式(16)使用时窗长度为7的PWVD做解析信号的时频分布，图 6为任意选取的第10段、第40段和第80段解析信号的时频分布情况，可见信号能量沿着瞬时频率集中，而噪声的能量分散在时频平面上，故取时频分布瞬时频率最大值恢复有效信号时可避免随机噪声干扰.
|
图 6 解析信号的时频分布 (a)第10段；(b)第40段；(c)第80段. Fig. 6 Time-frequency distribution of the analytic signal (a) The tenth segment; (b) The fortieth segment; (c) The eightieth segment. |
图 7为运用STFPF方法抑制随机噪声前后信号的时域图.图中蓝色线表示消除随机噪声前的信号，红色线表示消除随机噪声后的信号.图 7a为具有全部采样点的MRS全波信号消噪前后时域图，同样选取2~5 ms和112.6~115.2 ms对信号进行局部放大，见图 7b和c.经计算，消除随机噪声后信号的信噪比为18.02 dB，信噪比提高了23.19 dB，初始振幅E0=206.06 nV，初始振幅的拟合误差从滤波前的27.49 %提高到滤波后的3.03 %，平均横向弛豫时间T2*为154 ms，平均横向弛豫时间的拟合误差从滤波前的565.3%提高到滤波后的2.7%.可见，运用STFPF方法，随机噪声得到了有效的抑制，消噪效果优于传统叠加法.
|
图 7 运用STFPF方法抑制随机噪声前后信号的时域图 (a)具有全部采样点的含噪MRS全波信号；(b) 2~5 ms范围内的局部放大信号；(c) 112.6~115.2 ms范围内的局部放大信号. Fig. 7 The signal in time domain before and after random noise attenuation with STFPF algorithm (a) The complete MRS oscillating signal; (b) MRS oscillating signal at 2~5 ms; (c) MRS oscillating signal at 112.6~115.2 ms. |
为进一步比较两种方法消噪前后信号的细节特征，采用短时傅里叶变换对全波信号进行分析，得到信号的三维时频分布(图 8).
|
图 8 运用叠加法和STFPF方法抑制随机噪声前后信号的三维时频域图 (a)消噪前的信号；(b)运用叠加法消噪后的信号；(c)运用STFPF方法消噪后的信号. Fig. 8 Time-frequency domain of the signal before and after de-noising by traditional stacking method and STFPF method (a) Signal before de-noising; (b) Signal after de-noising by traditional stacking method; (c) Signal after de-noising by STFPF method. |
从图 8可以看出，信号中心频率为2128 Hz，随时间衰减.随机噪声能量均匀分布在整个时频平面上，对含随机噪声的信号进行8次叠加后，部分随机噪声被抑制，而运用STFPF算法对单次信号进行处理，随机噪声即得到了最佳抑制，信号成分被有效保留.可见，相比于传统的数据叠加方法，运用STFPF方法抑制MRS信号的随机噪声，不仅能取得较好的消噪效果，还能提高探测效率.
4 模型数据的反演解释为进一步讨论随机噪声对反演解释的影响，通过对比STFPF算法抑制随机噪声前后的地下水成像结果，分析有效抑制随机噪声对最终反演解释结果的意义.
首先，建立地下含水层仿真模型，见表 1.分别在地面以下5~15 m和25~35 m的位置处设置含水量为40%、平均横向弛豫时间为370 ms和520 ms的含水层.使用50 m单匝方形收发分离线圈，发射系统共设置15个脉冲矩，范围为0.01~5 As，按对数分布；假设测量时的地磁场为54, 721 nT，即Larmor频率为2330 Hz；地磁倾角为60°，地磁偏角为0°；进行反演解释时假设地下电阻率为500 Ωm，采用均匀半空间模型.首先计算地下水灵敏度核函数，再根据地下含水模型和灵敏度核函数，通过正演分析可获取实际地下含水层模型信号；随后，向信号中加入均值为0 nV，标准差分别为10~300 nV之间强度不等的随机噪声，运用SFTPF算法消除随机噪声，采用QT反演方法(Müller-Petke and Yaramanci, 2010)，通过共轭梯度算法(Scales, 1987)及边界限制条件(地下深度为0~70 m、弛豫时间为0.01~1 s、含水量为0~100%)对消噪前后的含水层深度、厚度、含水量大小及弛豫时间进行反演解释.
|
表 1 地下含水层模型 Table 1 Theoretical model consists of two-layer of aquifer |
图 9所示为受随机噪声干扰的MRS全波信号及其反演结果.可见，加入随机噪声后，信号淹没在随机噪声中(图 9a)，含水层的位置与实际存在偏差，浅部含水层反演得到了两个不同的弛豫时间值，深部含水层的弛豫时间也偏高，见图 9c.同时，数据拟合精度低.此时，运用STFPF算法抑制随机噪声，消除随机噪声后的MRS全波信号及反演结果如图 10所示.
|
图 9 随机噪声干扰下的MRS全波信号及反演结果 (a)含随机噪声的信号与脉冲矩关系；(b)拟合信号与脉冲矩关系；(c) QT反演结果. Fig. 9 Noisy MRS oscillating signal and best-fitting data estimation after inversion (a) The original noisy data; (b) The best-fitting data estimation after inversion; (c) The inversion result for the relaxation time and water content. |
|
图 10 抑制随机噪声后的MRS全波信号及反演结果 (a)抑制随机噪声后的信号与脉冲矩关系；(b)拟合信号与脉冲矩关系；(c) QT反演结果. Fig. 10 De-noised data by STFPF algorithm and best-fitting data estimation after inversion (a) The data after random noise reduction; (b) The best-fitting data estimation after inversion; (c) The inversion result for the relaxation time and water content. |
从图 10a可以看出，抑制随机噪声后各脉冲矩下的时域信号幅度有一个明显的逐渐变强再逐渐变弱的趋势，由于对应的脉冲矩较小，表明地下浅部有一个含水层.随后，再次出现变强再变弱的趋势，其对应的脉冲矩较大，说明深部也存在一个含水层，综合了弛豫时间和含水量大小的QT反演解释结果如图 10c所示.可见，地下5~15 m和25~35 m位置处各存在一个含水层，浅部含水层的弛豫时间以370 ms为主，在270~520 ms之间变化，含水量最高可达40%；深部含水层的弛豫时间以520 ms为主，在370~720 ms之间变化，含水量最高可达40%，反演得到的含水量与平均横向弛豫时间和仿真模型相符.根据反演结果得到的拟合信号如图 10b，拟合信号与消噪后的信号近乎相同.综合以上分析可知，运用STFPF方法消除信号中的随机噪声干扰后，解释结果可信度较高，与仿真的地下含水层模型基本一致.
5 结论复杂探测条件下如何有效地抑制随机噪声，实现微弱MRS信号可靠提取是当前磁共振数据处理的重点和难点，它将直接影响特征参数的准确提取和后续反演解释的准确性.本文基于时频峰值滤波法，针对全波磁共振信号数据量大、信号非线性强的特点，研究了分段时频峰值滤波算法，通过将信号分段处理，采用短窗长PWVD进行时频分析，实现了全波磁共振信号随机噪声的有效抑制.通过数值仿真实验和模型数据的反演对本文方法进行初步研究，得出如下结论：
(1) 针对随机噪声干扰问题，在时频分析的基础上，将时域信号变换为解析信号，根据信号和噪声的时频特征消除随机噪声，为磁共振地下水探测随机噪声的压制提供了理论依据；同时，针对传统叠加方法抑制随机噪声导致测量效率低、去噪效果有限问题，靶向单次测量信号，分段时频峰值滤波方法处理单次数据即可实现纯净磁共振信号的可靠提取，有效压制随机噪声，可提高探测效率.
(2) 对于高采样率全波磁共振信号数据点数多、直接做时频分析容易导致计算机内存溢出问题，将信号分段处理切实可行地实现了时频峰值滤波技术在磁共振全波信号中的应用.并针对磁共振全波信号非线性强的特点，采用短窗长PWVD对信号进行时频分析，仿真实验表明，当信噪比低至-5 dB时可有效提取信号，消除随机噪声后信噪比提高23.19 dB，信号的初始振幅拟合误差为3.03%，平均横向弛豫时间拟合误差为2.7%，消噪效果优于叠加法.
(3) 模型数据的反演表明，受随机噪声干扰，反演得到的含水层信息、平均横向弛豫时间等均与实际存在偏差；运用本文STFPF方法抑制随机噪声后含水量与平均横向弛豫时间和仿真模型相符，STFPF方法可较好的增强反演结果的可信度，为复杂噪声条件下应用磁共振方法提供技术支撑.
对于非线性的全波磁共振地下水探测信号，在运用本文STFPF方法时，通过适当的加窗处理，保证了信号在时窗内近似线性，窗长越小，线性度越高.但是，当信号的信噪比很低时，需增大窗长才能提高噪声抑制效果.二者矛盾难以调解.大量的仿真实验证明目前STFPF方法可较好的抑制信噪比-5 dB以上磁共振全波信号中的随机噪声.下一步将研究噪声干扰更严重情况下的随机噪声抑制问题.
致谢 感谢审稿专家和编辑部老师的大力支持.
Arnold M J, Roessgen M, Boashash B. 1994. Filtering real signals through frequency modulation and peak detection in the time-frequency plane. //Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing. Adelaide, SA, Australia: IEEE, 3: 345-348. https://ieeexplore.ieee.org/document/390019/?reload=true&arnumber=390019
|
Arnold M J, Boashash B. 1994. Time-frequency peak filtering: Analysis and implementation details. //Proceedings of the IEEE-SP International Symposium on Time-Frequency and Time-Scale Analysis. Philadelphia, PA, USA, USA: IEEE, 256-259. https://www.researchgate.net/publication/3613374_Time-frequency_peak_filtering_analysis_and_implementation_details
|
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. DOI:10.1007/s10712-014-9304-0 |
Boashash B, Mesbah M. 2004. Signal enhancement by time-frequency peak filtering. IEEE Transactions on Signal Processing, 52(4): 929-937. DOI:10.1109/TSP.2004.823510 |
Costabel S, Müller-Petke M. 2014. Despiking of magnetic resonance signals in time and wavelet domains. Near Surface Geophysics, 12(4): 185-197. |
Dalgaard E, Auken E, Larsen J J. 2012. Adaptive noise cancelling of multichannel magnetic resonance sounding signals. Geophysical Journal International, 191(1): 88-100. DOI:10.1111/gji.2012.191.issue-1 |
Dalgaard E, Müller-Petke M, Auken E. 2016. Enhancing SNMR model resolution by selecting an optimum combination of pulse moments, stacking, and gating. Near Surface Geophysics, 14(3): 243-253. |
Davis A C, Dlugosch R, Queitsch M, et al. 2014. First evidence of detecting surface nuclear magnetic resonance signals using a compact B-field senso. Geophysical Research Letters, 41(12): 4222-4229. DOI:10.1002/2014GL060150 |
Davis A C, Müller-Petke M, Dlugosch R, et al. 2015. First evidence of T2* in SNMR measurements with SQUID sensors. //85th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
DePasquale G, Mohnke O. 2014. Numerical study of prepolarized surface nuclear magnetic resonance in the vadose zone. Vadose Zone Journal, 13(11): 1-9. |
Garwood M, DelaBarre L. 2001. The return of the frequency sweep:Designing adiabatic pulses for contemporary NMR. Journal of Magnetic Resonance, 153(2): 155-177. DOI:10.1006/jmre.2001.2340 |
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 |
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. 2016. Model-based subtraction of spikes from surface nuclear magnetic resonance data. Geophysics, 81(4): WB1-WB8. DOI:10.1190/geo2015-0442.1 |
Legchenko A V, Valla P. 2002. A review of the basic principles for proton magnetic resonance sounding measurements. Journal of Applied Geophysics, 50(1-2): 3-19. DOI:10.1016/S0926-9851(02)00127-1 |
Legchenko A V, Baltassat J M, Beauce A, et al. 2002. Nuclear magnetic resonance as a geophysical tool for hydrogeologists. Journal of Applied Geophysics, 50(1-2): 21-46. DOI:10.1016/S0926-9851(02)00128-3 |
Li Y, Lin H B, Yang B J, et al. 2009. The influence of limited linearization of time window on TFPT under the strong noise background. Chinese Journal of Geophysics (in Chinese), 52(7): 1899-1906. |
Lin H B, Li Y, Xu X C, et al. 2011. Time-frequency peak filtering algorithm for reducing the discrete error. Journal of Jilin University (Earth Science Edition) (in Chinese), 41(2): 572-578. |
Lin J. 2010. Situation and progress of nuclear magnetic resonance technique for groundwater investigations. Progress in Geophysics (in Chinese), 25(2): 681-691. |
Lin J, Duan Q M, Wang Y J, et al. 2011. Theory and Design of Magnetic Resonance Sounding Instrument for Groundwater Detection and Its Applications. Beijing: Science Press.
|
Lin J, Jiang C D, Duan Q M, et al. 2012. The situation and progress of magnetic resonance sounding for groundwater investigations and underground applications. Journal of Jilin University (Earth Science Edition) (in Chinese), 42(5): 1560-1570. |
Lin J, Zhang Y, Zhang S Y, et al. 2016. Progress of magnetic resonance sounding for groundwater investigation under high-level electromagnetic interference. Journal of Jilin University (Earth Science Edition) (in Chinese), 46(4): 1221-1230. |
Lin T T, Shi W L, Qi X, et al. 2014. Full-wave recording of MRS groundwater exploration system. Journal of Jilin University (Engineering and Science Edition) (in Chinese), 44(4): 95-102. |
Lin T T, Yang Y J, Teng F, et al. 2018a. Enabling surface nuclear magnetic resonance at high-noise environments using a pre-polarization pulse. Geophysical Journal International, 212(2): 1463-1467. DOI:10.1093/gji/ggx490 |
Lin T T, Zhang Y, Yi X F, et al. 2018b. Time-frequency peak filtering for random noise attenuation of magnetic resonance sounding signal. Geophysical Journal International, 213(2): 723-738. |
Lu Q H, Wu T B, Lin J. 2009. A research report on development of instrument science for geophysics. Progress in Geophysics (in Chinese), 24(2): 750-758. |
Müller-Petke M, Yaramanci U. 2010. QT inversion-Comprehensive use of the complete surface NMR data set. Geophysics, 75(4): WA199-WA209. DOI:10.1190/1.3471523 |
Müller-Petke M, Costabel S. 2014. Comparison and optimal parameter settings of referencebased harmonic noise cancellation in time and frequency domains for surface-NMR. Near Surface Geophysics, 12(2): 199-210. |
Packard M E. 1954. Free nuclear induction in the Earth's magnetic field. Phys Rev, 93(4): 941. |
Pan Y L, Zhang C D. 2000. Theories and Methods of Surface Nuclear Magnetic Resonance. Beijing: Publishing House of Geology University of China.
|
Qu Y X. 2015. Development of broadband array full wave MRS acquisition system[Master's thesis](in Chinese). Changchun: Jilin University.
|
Roessgen M, Boashash B. 1994. Time-frequency peak filtering applied to FSK signals. //Proceedings of the IEEE-SP International Symposium on Time-Frequency and Time-Scale Analysis. Philadelphia, PA, USA, USA: IEEE, 516-519.
|
Scales J A. 1987. Tomographic inversion via the conjugate gradient method. Geophysics, 52(2): 179-185. DOI:10.1190/1.1442293 |
Tannús A, Garwood M. 1997. Adiabatic pulses. NMR in Biomedicine, 10(8): 423-434. DOI:10.1002/(ISSN)1099-1492 |
Walsh D O. 2008. Multi-channel surface NMR instrumentation and software for 1D/2D groundwater investigation. 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 (in Chinese), 59(6): 2290-2301. DOI:10.6038/cjg20160631 |
Wang Z X, Rong L L, Lin J, et al. 2010. FID signal detection based on DLIA sampled by quadruple larmor frequency. Journal of Data Acquisition and Processing, 25(5): 626-630. |
Weichman P B, Lavely E M, Ritzwoller M H. 2000. Theory of surface nuclear magnetic resonance with applications to geophysical imaging problems. Physical Review E, 62(1): 1290-1312. DOI:10.1103/PhysRevE.62.1290 |
Wu N, Li Y, Tian Y N, et al. 2016. Trace-transform-based time-frequency filtering for seismic signal enhancement in Northeast China. Comptes Rendus Geoscience, 348(5): 360-367. DOI:10.1016/j.crte.2016.02.001 |
Zhang R, Hu X Y, Yang D K, et al. 2006. Review of development of surface nuclear magnetic resonance. Progress in Geophysics (in Chinese), 2006, 21(1): 284-289. |
李月, 林红波, 杨宝俊, 等. 2009. 强随机噪声条件下时窗类型局部线性化对TFPF技术的影响. 地球物理学报, 52(7): 1899-1906. DOI:10.3969/j.issn.0001-5733.2009.07.025 |
林红波, 李月, 徐学纯, 等. 2011. 减小离散误差的时频峰值滤波算法. 吉林大学学报(地球科学版), 41(2): 572-578. |
林君. 2010. 核磁共振找水技术的研究现状与发展趋势. 地球物理学进展, 25(2): 681-691. DOI:10.3969/j.issn.1004-2903.2010.02.043 |
林君, 段清明, 王应吉, 等. 2011. 核磁共振找水仪原理与应用. 北京: 科学出版社.
|
林君, 蒋川东, 段清明, 等. 2012. 复杂条件下地下水磁共振探测与灾害水源探查研究进展. 吉林大学学报(地球科学版), 42(5): 1560-1570. |
林君, 张扬, 张思远, 等. 2016. 强电磁干扰下磁共振地下水探测噪声压制方法研究进展. 吉林大学学报(地球科学版), 46(4): 1221-1230. |
林婷婷, 史文龙, 齐鑫, 等. 2014. 核磁共振地下水探测全波收录系统. 吉林大学学报(工学版), 44(4): 95-102. |
陆其鹄, 吴天彪, 林君. 2009. 地球物理仪器学科发展研究报告. 地球物理学进展, 24(2): 750-758. DOI:10.3969/j.issn.1004-2903.2009.02.053 |
潘玉玲, 张昌达. 2000. 地面核磁共振找水理论和方法. 北京: 中国地质大学出版社.
|
曲永星. 2015. 阵列式宽带核磁共振全波采集系统的研制[硕士论文]. 长春: 吉林大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-1015591426.htm
|
万玲, 张扬, 林君, 等. 2016. 基于能量运算的磁共振信号尖峰噪声抑制方法. 地球物理学报, 59(6): 2290-2301. DOI:10.6038/cjg20160631 |
王中兴, 荣亮亮, 林君, 等. 2010. 基于4倍频采样的数字正交FID信号检测技术. 数据采集与处理, 25(5): 626-630. DOI:10.3969/j.issn.1004-9037.2010.05.016 |
张荣, 胡祥云, 杨迪琨, 等. 2006. 地面核磁共振技术发展述评. 地球物理学进展, 21(1): 284-289. DOI:10.3969/j.issn.1004-2903.2006.01.043 |