自1990年起，地面磁共振技术(Magnetic Resonance Sounding, MRS)已经成为非侵入式直接探测地下水资源最有效的技术之一, 在全球范围内广泛应用(Yaramanci, 2000; Legchenko et al., 2002, 2010).该方法利用拉莫尔频率进行激励，探测地下水的弛豫效应，能够进行地下水层的定性定量解释，系统评价含水量、孔隙度、渗透率等水文地质参数.然而，由于地面核磁共振探测获得的信号及其微弱，仅能到达纳伏级(10^-9 V)，低信噪比成为了限制该方法广泛应用的主要障碍.大量研究发现工频噪声干扰和尖峰噪声干扰是对磁共振信号影响最大的两类噪声.仪器系统和相关处理方法在不断的改进以提高系统探测信噪比.

第一代磁共振探测仪器仅有一个探测通道(Legchenko and Valla, 2003)，即激发和采集使用同一线圈.陷波器、8字形线圈被用来进行工频噪声和尖峰噪声的消除，取得了一定的效果，但也存在显著缺陷：陷波器会使信号产生失真，而8字型线圈则会降低探测深度并导致核函数正演计算难于求解.随着多通道仪器的诞生(Walsh, 2008;林婷婷等, 2013；Dalgaard et al., 2014), 可有效避免单通道仪器滤波方法的缺陷.在多通道仪器中，通过发射线圈激励，主通道仍作为信号采集通道，获取地下水磁共振弛豫信号.而其他的一个或多个通道远离激发源线圈(通常为线圈直径2倍以上距离)，因此仅能采集噪声.一般在实验前，可通过噪声分析仪器分析探测场地周边主要噪声干扰源，如电力线、电网等.再有针对性的实施参考通道布设.在后续数据处理中，当参考通道的噪声和主通道的噪声具有相关性时，噪声可以被有效消减，地面磁共振信号就可以被提取出来.多通道滤波技术对工频及固定频率干扰噪声消除效果显著，在复杂环境下受到限制，如尖峰噪声干扰严重, 会影响数据之间的相关性(田宝凤等, 2015).因此，去除尖峰噪声对获得高信噪比的地面磁共振信号显得极为重要.

针对尖峰噪声，国内外学者开展了一系列相关研究.蒋川东等(2011) 使用统计叠加的方法去除尖峰噪声时发现信号质量有所改善，但该方法必须要有大量的数据才能够应用，并且对幅值较小的尖峰噪声的消除效果并不理想，影响后续反演解释.Dalgaard(2012)和万玲等(2016) 采用非线性能量算子方法识别及去除尖峰噪声.该方法利用不含尖峰噪声其他组数据来替换尖峰噪声所污染的地面磁共振信号，也需要多组数据才能使用.为仅采用一组数据进行尖峰噪声消减，提高处理效率，Costabel和Müller-Petke(2014) 研究了小波多尺度分解硬阈值算法，证实了小波方法的先进性，取得了一定效果.但是由于小波多尺度分解硬阈值方法中硬阈值算法设定阈值的固有缺陷，会导致数据引起震荡，引发的伪吉布斯效应，在全波数据处理中会导致弛豫数据较大的反演拟合误差.因此，学者们预测，更为先进的小波处理手段将在消除尖峰噪声方具有更好的性能.

本文中，我们首次提出了结合压缩小波变换(Synchrosqueezing Wavelet Transform，SWT)和非线性阈值处理(Nonlinear Thresholding，NT)方法联合处理尖峰噪声.通过SWT获得的压缩小波系数来识别尖峰噪声，而NT则用来剔除尖峰噪声并恢复由于尖峰噪声干扰抑制而丢失的地面磁共振信号.这种新算法仅需要一组数据即可去除尖峰噪声并且重构信号.在仿真数据和实测数据中，本文提出的方法均展现了良好的消噪性能.

傅里叶变换是传统的信号分析处理工具，主要用于全局变换，不适用于分析非平稳信号(如MRS信号).因此，研究者提出了小波变换理论，在时域和频域的采样密度可以自动调节(Chavez-Roman and Ponomaryov, 2014).在高频时，调整时域的窗口变窄来提高分辨率，低频时，时域窗口变宽增加频域分辨率(Iqbal et al., 2013).因此，小波变换是理想的非平稳时变信号处理方法，可以用来处理地面磁共振信号.

在L²(R)空间中，设ψ(t)为一平方可积函数，即ψ(t)∈L²(R)，若其傅里叶变换满足条件：

(1)

则称ψ(t)为一个基本小波或小波母函数，其中为ψ(t)的傅里叶变换.式(1) 为对ψ(t)提出的容许条件，该条件保证了小波变换的反变换存在，因为任何变换都必须存在反变换才有实际意义.

将小波母函数ψ(t)进行伸缩和平移，就可以得到小波基函数ψ_{a, τ}(t)，公式为

(2)

其中a为伸缩因子，τ为平移因子，ψ_{a, τ}(t)为依赖于参数a和τ的小波基函数.伸缩因子a和平移因子τ是连续变化的值，ψ_{a, τ}(t)是由同一母函数ψ(t)经伸缩和平移后得到的一组函数序列.小波基函数正是依靠伸缩和平移来实现不同时刻信号的局部分析的.

SWT是Daubechies等(2011) 提出的基于连续小波变换的时频再分配方法，首先用于地震数据处理(Herrera et al., 2014).SWT可以使时频曲线更加精细和清晰，从而改善时频曲线的精度(汪祥莉等, 2016).常见的基小波类型有Haar小波、Daubechies小波、Meyer小波和Morlet小波等.本文采用Morlet小波作为基小波，原因如下：(1) Morlet小波是由ψ(t)=e^i_μ0te^－t2/2给出，其中μ₀是Morlet小波基频.Morlet小波是由高斯函数调制而来的，在时域有较强的衰减性.(2) Morlet小波的时域窗口中心t³=0，窗口半径.频域窗口中心ω³=ω₀，窗口半径由Heisenberg测不准原理，Morlet小波的窗口面积达到测不准原理中的下界1/2.故Morlet小波更适合于刻画信号时域局部特性，适用于磁共振信号的处理，会在非线性处理阶段取得良好的效果.

由接收线圈接收到的地面磁共振信号被描述为

(3)

其中t为时间，E₀、T₂^*、ω₀和φ分别为地面磁共振信号的初始振幅，横向弛豫时间，拉莫尔频率和初始相位.

将L²(R)空间中的函数f(t)在小波基下展开，称为函数f(t)的连续小波变换, 其表达式为

(4)

式(4) 的WT_f(a, τ)可理解为用一组分析宽度不断变化的基函数对f(t)作分析，这一变化正好适应了我们对信号分析时在不同频率范围所需要不同的分辨率这一要求.

通过证明可以得到，若采用的小波函数满足容许条件，则连续小波变换存在逆变换，逆变换公式为

(5)

式(5) 中，为对ψ(t)提出的容许条件.

而小波基函数ψ_{a, τ}(t)式中的因子是为了保证在不同的尺度a时，ψ_{a, τ}(t)始终能和母函数ψ(t)有着相同的能量，即

(6)

令，则dt=adt′，这样上式的积分即等于

令f(t)的傅里叶变换为f(Ω)，ψ(t)的傅里叶变换为ψ(Ω)，由傅里叶变换的性质，ψ_{a, τ}(t)的傅里叶变换为

(7)

由Parsevals定理, 式(4) 可重新表示为

(8)

将地面磁共振信号f(t)在小波基下展开，得到f(t)的压缩小波系数为

(9)

其中WT_f(a_k, τ)为地面磁共振信号f(t)经过连续小波变换后获得的小波系数，a_k是伸缩因子a的离散值，(Δa)_k=a_k-a_k-1，ω_l是f(t)的中心频率，SWT在(ω_l, τ)坐标下重组小波系数W_f(a, τ)≠0，ω_f(a, τ)是在范围内，其中Δω=ω_l－ω_l－1.

传统的小波变换阈值方法中包括软阈值和硬阈值.其中软阈值函数为

(10)

其中f(t)为输入信号，T为预先设定好的阈值，sign是符号函数，

与软阈值相比，硬阈值收缩性能更强，在去除尖峰噪声时，信号重建过程中会产生伪吉布斯现象.硬阈值函数为

(11)

其中f(t)为输入信号，T为预先设定好的阈值.

为克服硬阈值和软阈值的固有缺陷，本文引入非线性阈值处理技术(Li et al., 2000).该方法相对于硬阈值算法及软阈值算法的优势在于，它能聚焦信号局部特性，利用调整阈值衰减相应的系数对地面磁共振信号进行处理，方法更为平滑.不会有软、硬阈值算法设置阈值导致的部分信号失真.但目前，利用本文算法在处理尖峰噪声时(特别是采集尾部信号，磁共振衰减平缓)要求噪声持续时间在40 ms以内.针对持续时间过长的尖峰噪声，精准的阈值系数获取较为困难，尚需要进一步研究和算法改进.

NT算法的函数可以表示为

(12)

其中ω_{j, k}是压缩小波系数，是去除尖峰和进行处理后的新小波系数，k是地面磁共振信号的长度，λ是所设定的阈值，公式为

(13)

其中σ是地面磁共振信号的标准偏差(Donoho et al., 1995).

用NT方法获得新的小波系数后对地面磁共振信号进行小波重构，就会使被尖峰噪声干扰的数据得到相应恢复.

为验证本文提出算法的有效性，首先，作者用仅含有一个尖峰的数据说明算法的执行过程.MRS信号用MATLAB 2015软件通过式(3) 产生，产生的仿真信号的特征参数为初始振幅E₀=300 nV, 拉莫尔频率为2326 Hz, 初始相位φ=2 rad和横向弛豫时间T₂^*=300 ms.而尖峰噪声产生公式为

(14)

其中A是尖峰噪声f_s的幅值，T_w是尖峰噪声f_s的持续时间，ω是尖峰噪声f_s的频率.

为了观察SWT-NT算法处理尖峰噪声的效果，引入信噪比计算函数为

(15)

其中，s(i)代表信号，n(i)代表噪声，M代表采样点数.

仿真数据含有10000个等精度采样点、采样频率为66667 Hz的单次采集的原始信号.如图 1a所示，在图中的200 ms处加一个幅值2000 nV，持续时间5 ms的尖峰噪声，使得信噪比SNR=1.1602 dB.使用SWT-NT法联合去除尖峰噪声及处理该信号，主要包括三个步骤：

图 1

Fig. 1

图 1 仿真了本算法对单个尖峰进行处理的情况 (a)加入一个尖峰噪声的地面磁共振信号；(b)经过本文提出的算法处理后的时域图；(c)利用连续小波变换对(a)进行时频分析的结果; (d)利用连续小波变换对(b)进行时频分析结果. Fig. 1 Simulations of a single spike rejection based on the SWTNT algorithm (a) Time domain of the SNMR signal with one spike added; (b) Time domain of the SNMR signal after processing by the SWTNT algorithm; (c) Time-frequency analysis by CWT to the signals shown in (a); (d) Time-frequency analysis by CWT to the signals shown in (b).

(1) 针对带有尖峰噪声的仿真信号，进行小波分解，能够将信号分为不同尺度，采用Morlet小波来得到压缩小波系数，信号的压缩小波系数与尖峰噪声的压缩小波系数特征不同，即尖峰噪声对应的小波系数幅值小.

(2) 将(1) 中得到的不同小波尺度下的压缩小波系数通过式(12) 的方程获得相应的非线性阈值处理函数，利用NT算法在不同的小波域中去除尖峰噪声同时通过获得的新阈值衰减相应的系数对地面磁共振信号进行复原.

(3) 将各个尺度上的信号进行小波重构，如图 1b所示，得到本算法处理后的信号.应用SWT-NT算法处理后，信噪比SNR= 24.2300 dB.

图 1c是经过时频域变换后的结果，能够很好的观测到尖峰噪声的集中情况, 而图 1d是消除尖峰后的信号数据，可以观测到信号能量随弛豫衰减过程均匀分布，表明尖峰噪声已被有效滤除.

在实际探测中，磁共振信号易被多个尖峰噪声干扰.为验证该算法在多个尖峰噪声干扰情况下的消噪能力：图 2a、c和e中加入三个幅值以及持续时间不同的尖峰噪声.其中，产生的仿真信号的特征参数为初始振幅E₀=300 nV, 拉莫尔频率为2326 Hz, 初始相位ϕ=2 rad和横向弛豫时间T₂^*=300 ms.在图 2a的50 ms、100 ms和200 ms处分别加入幅值为2000 nV、1000 nV和800 nV以及持续时间为5 ms、5 ms和10 ms的尖峰噪声，信噪比SNR=0.3593 dB.在图 2c中的噪声分布与图 2a相同，但在50 ms处的尖峰噪声持续时间为10 ms，信噪比SNR=0.2640 dB，考察一个尖峰脉冲持续时间延长对消噪效果的影响.而图 2e在图 2c的基础上，又放大了第一个尖脉冲的幅值，达到4000 nV，信噪比SNR=0.1192 dB.

图 2

Fig. 2

图 2 对本算法处理三个尖峰噪声的情况进行了仿真 左边栏是加入了三个尖峰噪声的地面磁共振信号，而右边栏是经过本算法处理后的信号. (a)加入三个尖峰噪声的地面磁共振信号；(c)将(a)中第一个尖峰噪声的持续时间由5 ms改为10 ms；(e)将(c)中第一个尖峰噪声的幅度由2000 nV改为4000 nV；(b)，(d)和(f)是本算法分别对(a)，(c)和(e)进行处理后所得. Fig. 2 Simulations of three spike rejection based on the SWTNT algorithms The left-hand column is the SNMR signal with three spikes added. The right-hand column is the SNMR signal after processing with the SWTNT algorithm; (a) The SNMR signal with three spikes added; (c) The duration of first spike duration in (a) is changed from from 5 ms to 10 ms; (e) Amplitude of the first spike amplitude in (c) is changed from 2000 nV to 4000nV; (b), (d) and (f) are the SNMR signals after being processed by the SWTNT algorithm.

因此，图 2不仅用来验证算法在多个尖峰噪声干扰情况的消噪能力，而且比较了干扰幅度不同以及持续时间不同的尖峰噪声的消噪能力.采用SWT-NT算法对含尖峰噪声信号进行处理，图 2b、d和f为对应的消噪结果.应用SWT-NT算法处理后，信噪比分别提升至17.0150 dB、11.4423 dB和8.0507 dB.从消噪结果图中可以明显看出，应用SWT-NT算法处理后，不论尖峰噪声幅值的大小、持续时间的长短均被去除.同时，通过阈值衰减相应的系数对地面磁共振缺失的信号进行重构，恢复了原始信号.然而从信噪比的情况来看，信号中尖峰噪声的持续时间和幅度会在一定程度影响本文提出算法的效果，即尖峰持续时间越长，对算法拟合要求精度越高.

为了更接近实际探测中的噪声环境，我们进一步进行了多尖峰噪声干扰的仿真.仿真了被幅值不同、持续时间不同、多个复杂的尖峰噪声环境所破坏的信号.信号仍为含有10000个等精度采样点、采样频率为66667 Hz的单次采集信号.如图 3a所示.其中，在仿真信号的10、50、100、150、200、230 ms和280 ms处分别加入幅值为140、180、40、40、90、45 uV和35 uV以及持续时间为3、3、5、3、5、40 ms和5 ms的尖峰噪声，信噪比SNR=-20.6405 dB.

图 3

Fig. 3

图 3 利用本算法处理含有七个尖峰噪声的仿真 (a)是加入尖峰噪声的地面磁共振信号；(b), (d)和(f)是分别应用SWT-NT算法，小波多尺度分解硬阈值方法和小波多尺度分解软阈值方法对(a)进行处理结果；(c), (e)和(g)是分别对(b), (d)和(f)进行时频分析后所得的时频分析结果. Fig. 3 Simulations of seven spike rejection based on the SWTNT algorithm (a) Time domain of the SNMR signal with seven spikes added; (b), (d) and (f) represent the time domains of the SNMR signal in (a) after processing by the SWTNT algorithm, the multi-scale wavelet hard threshold algorithm, the multi-scale wavelet soft threshold algorithm, respectively; (c), (e) and (g) represent the time-frequency analysis to the signals shown in (b), (d) and (f), respectively.

在这一仿真中，我们对比了Costabel等提出的小波多尺度分解硬阈值方法, 小波多尺度分解软阈值方法和本文提出的SWT-NT算法.Costabel等提出的小波多尺度分解硬阈值方法是将含有尖峰噪声的核磁数据分解到多个小波尺度内，在所有的小波尺度中利用设定好的阈值对尖峰噪声进行去除.由于在不同的小波尺度中采用的阈值都相同，所以最后进行小波重构后难以获得好的噪声压制效果，部分尖峰噪声残余仍然突出，或过分压制导信号失真.而本文提出的SWT-NT算法能够在不同的压缩小波尺度内进行阈值变换，即利用通过公式(12) 得到不同阈值，最后再经过小波重构后得到较好的处理结果.

应用SWT-NT算法对图 3a所示的仿真含尖峰噪声信号进行处理，为对比本文方法和软、硬阈值小波变换的差异，在这里详细说明采用非线性阈值处理方法针对尖峰噪声处理的阈值取值情况.依据公式(12)，在第一个尖峰噪声阈值取值自267 nV减少至262 nV，保证磁共振信号的弛豫效应.同理，第二个尖峰阈值自202 nV减少至195 nV，第三个尖峰噪声阈值自121 nV减少至116 nV，第四个尖峰阈值自65 nV减少至64 nV，第五个尖峰噪声阈值自43 nV减少至40 nV，第六个尖峰噪声阈值自29 nV减少至20 nV，第七个尖峰噪声阈值为16 nV.改善了传统阈值选取方法固定阈值的弊端，同时兼顾了磁共振信号衰减特性.图 3b为对应的消噪结果，信噪比从之前的SNR=-20.6405 dB，提高至4.2426 dB，从弛豫衰减的拟合情况看，信号良好拟合.图 3c时频分析的结果进一步解释信号的能量随弛豫衰减过程均匀分布，表明尖峰噪声被有效滤除.但是在尾部的一个持续时间为40 ms (图 3b红色线框中)，幅值为45 μV的尖峰噪声处重构信号的效果并不十分理想.证明后续研究中需要进一步改进阈值函数，比如结合尖峰噪声出现前及结束后的信号幅值特性改进补偿函数的参数进而提高重构数据的效果.相比之下，应用小波多尺度分解硬阈值方法，阈值设置为363 nV(该阈值为信噪比最优化结果)，根据公式(11) 对含尖峰噪声信号进行处理时(图 3d)，伪吉布斯效应明显，引起信号的震荡，信噪比提高至-8.3050 dB.从图 3e时频分析结果来看，信号的能量随弛豫衰减过程分布紊乱，表明尖峰噪声未被有效滤除，同时核磁共振信号丢失较多细节.利用小波多尺度分解软阈值方法，阈值设置为246 nV，根据公式(10) 对仿真含尖峰噪声信号进行处理(图 3f)，信噪比仅提高至-11.4862 dB，从图 3g时频分析的结果来看，信号的能量随弛豫衰减过程分布失衡，磁共振信号畸变失真.综上所述，本文提出的SWT-NT算法可以针对复杂噪声环境的多个尖峰噪声进行处理, 并且最大程度还原信号.

在地面MRS方法应用中，一般首先利用消噪方法处理磁共振数据，之后再进行后续反演解释，所以利用消噪技术提取的特征参数对反演解释十分重要.在表 1中针对图 3的仿真信号以及经过新算法处理后的信号的特征参数(初始振幅E₀，拉莫尔频率f_L，初始相位Ф和横向弛豫时间T₂^*)进行对比.表格中的数据表明，经过SWT-NT处理后信号的特征参数恢复具有很高的精度，为后续高精度的反演解释奠定基础.

表 1(Table 1)

表 1 比较了仿真的地面磁共振信号和经过本算法处理尖峰噪声后的信号的特征参数 Table 1 Comparison of characteristic parameters of simulated ground MRS signal and the signal after processing by SWTNT

表 1 比较了仿真的地面磁共振信号和经过本算法处理尖峰噪声后的信号的特征参数 Table 1 Comparison of characteristic parameters of simulated ground MRS signal and the signal after processing by SWTNT

为了便于读者明确本文算法的处理流程，图 4详细描述了SWT-NT算法原理.首先，针对带有尖峰噪声的MRS信号, 利用SWT算法获得不同小波尺度下的压缩小波系数，同时得到NT算法所需要的参数.利用尖峰噪声的压缩小波系数与地面磁共振信号的压缩小波系数具有不同特点(噪声对应的小波系数幅值小)，可以在小波域内对尖峰噪声进行识别，获得不同信号处理的局部阈值.在NT算法处理后，SWT算法会对地面磁共振信号进行重构.通过以上步骤就会对混入地面磁共振信号中的尖峰噪声去除的同时通过阈值衰减相应的系数对地面磁共振缺失的信号进行重构，恢复原始信号.

图 4

Fig. 4

图 4 应用SWT-NT算法处理尖峰噪声的原理图 Fig. 4 Schematic showing principle of the SWTNT algorithm applied to despiking

为了验证本文提出的算法在实际应用中的有效性，课题组进行了多次野外实验的数据采集与处理.测试数据均由吉林大学自主研制的多通道地面磁共振探测仪器完成(Lin et al., 2015).选择的实验场地位于长春市郊区太平水库附近，测区地表有防护林植被覆盖.在测试地点附近有砖厂施工会引入工频噪声，此外，有农用拖拉机在探测线圈附近经过引入尖峰噪声.线圈铺设方式为方型同一线圈，边长100 m×100 m，采用16个脉冲矩探测，发射范围从0.2 As到6 As，单次实验发射时间为40 ms，采集时间为256 ms.当地拉莫尔频率为2330 Hz.

图 5a和图 5c为第八组脉冲矩0.792 As的测量数据，蓝色线为仪器采集到的原始信号时域波形图.图中可见，图 5a为单个尖峰噪声干扰数据，振幅达2000 nV，持续时间为30 ms.经过式(15) 对图 5a实测信号进行信噪比的计算，得知信噪比SNR= -37.6274 dB.图 5c为多个尖峰噪声干扰数据，此外，还混入部分工频噪声.信噪比SNR=-29.4156 dB.频率为2330的信号均被噪声干扰及淹没.

图 5

Fig. 5

图 5 利用本算法来处理实测数据中的尖峰噪声 (a)单个尖峰噪声情况下，实测数据(蓝色)和经过本算法处理尖峰噪声以及自适应算法去除工频噪声后的数据(红色)的对比；(b)是与(a)对应的频谱对比；(c)多个尖峰噪声情况下，实测数据(蓝色)和经过本算法处理尖峰噪声以及自适应算法去除工频噪声后的数据(红色)的对比；(d)是与(c)对应的频谱对比. Fig. 5 SNMR signal despiking to field data (a) In the condition of single spike, comparison of measured data (blue) and the data after despiking by the SWTNT and self-adaptive algorithms (red); (b) Spectrum comparison corresponding to (a); (c) In the condition of multiple spikes, comparison of measured data (blue) and the data after despiking by the SWTNT and self-adaptive algorithms (red); (d) Spectrum comparison corresponding to (c).

依据本文所述算法，针对实测数据进行尖峰噪声剔除.在去除尖峰噪声和工频噪声后，信噪比分别提升至1.1357 dB和8.5482 dB.消噪后的时域信号如图 5a和c中的红色线所示，实测信号中的尖峰噪声已经明显被剔除，通过非线性阈值处理方法后的MRS信号呈现指数形式衰减.再从频域上进行对比分析，如图 5b和d所示，信号频率2330 Hz清晰可见.此外，在频谱中部分工频点还存在噪声残余，可配合成熟的陷波器、自适应陷波器进行进一步处理.本次实验证实了本文提出的算法能够有效用于实际数据处理中.在应用SWT-NT算法处理实测数据的尖峰噪声时，应注意它的适用条件，即尖峰噪声持续时间在40 ms以内.

本次实验是在长春市郊区烧锅地区进行磁共振地下水探测实验.测试地点附近有工地施工会引入工频噪声，此外，探测过程中，摩托车在线圈附近经过，导致尖峰噪声引入.发射线圈采用单匝25 m×75 m矩形线圈，接收线圈采用3个25 m×25 m双匝方形线圈在发射线圈内边对边摆放.此种线圈布设方法能够实现二维地下水探测.在试验中发射脉冲矩为20个, 发射范围从0.2 As到10 As，单次实验发射时间为40 ms，采集时间为256 ms.当地拉莫尔频率2332 Hz.

图 6a处理第六组脉冲矩576Ams的测量数据.蓝色线为仪器采集到的原始信号时域波形图，从原始信号的时域和频域的波形图上看出MRS信号由于尖峰噪声，工频噪声和随机噪声等噪声干扰而产生了较严重的失真.在图 6a中尖峰噪声数量多达5个，第一个尖峰噪声在10 ms处出现，振幅最大达600 nV，持续时间为20 ms，最后一个尖峰噪声在245 ms处出现，持续时间最长为35 ms，振幅为120 nV.本组测试结果的尖峰噪声极为复杂，且连续出现，均持续大于20 ms以上，严重导致信号畸变，信噪比SNR=-34.5872 dB.

图 6

Fig. 6

图 6 利用本算法来处理实测数据中的尖峰噪声以及对处理后的数据进行反演解释 (a)实测数据(蓝色)和经过SWT-NT算法处理尖峰噪声(红色)的对比；(b)是与(a)对应的频谱对比；(c)为初始振幅与脉冲矩的对应曲线；(d)为地下含水量二维全弛豫数据反演. Fig. 6 SNMR signal despiking for field data and inversion results (a) Comparison of measured data of and the data after despiking by the SWTNT algorithm; (b) Spectrum comparison to (a); (c) Curves of E0-Q; (d) 2D inversion for water content of subsurface.

依据本文所述的SWT-NT算法，进行尖峰噪声的剔除.消噪后的时域信号如图 6a中的红色线所示，时域上磁共振信号弛豫衰减特性已有恢复、但部分的信号依然受到工频噪声影响.再从频域上进行观察，如图 6b红色线所示，2330 Hz信号频率点经过尖峰噪声处理后，信号凸出.2350 Hz、2330 Hz处工频谐波噪声仍然为剩余噪声的主要频率成分，我们采用陷波器进行工频噪声的进一步处理.图 6c和图 6d是对图 6a中去除尖峰噪声后的信号进行反演解释.其中，图 6c为E₀-Q拟合图，coil1~3是做二维反演所需的边对边摆放的三个线圈，用来获得地下的平面以及深度组合成的数据来实现对地下水信息的二维全弛豫数据反演.图 6d为地下含水量二维全弛豫反演，其中的纵轴Z代表地下深度，而横轴X是以线圈的中心点为原点的水平距离.图中15~25 m、以及40~60 m可见两层层状含水层.第一含水层与当地水文地质资料情况完全相符(Lin et al., 2015).第二层水中coil1下未能完全反映层状信息，这与边对边铺设线圈的分辨率相关.通过图 6说明了本文所述的SWT-NT算法保障了后续反演解释的有效进行，对地面磁共振探测地下水来了解地下水的分布以及含量等水文地质信息具有重要的意义.

本文基于压缩小波变换和非线性阈值处理方法，对地面磁共振信号中的尖峰噪声的去除进行了研究，提出了一种可以在去除尖峰噪声的同时通过阈值衰减相应的系数对地面磁共振缺失的信号进行重构，恢复原始信号的方法.与国内外近年来所研究并推广的能量算子以及统计叠加相比，本文所提出的方法仅通过一次采集数据就可完成尖峰噪声剔除，是一种高效而准确的数据处理方法.与国际上采用硬阈值进行尖峰噪声消除等前沿技术相比较，能够最大程度还原磁共振信号，为后续工频噪声消除、地下含水层反演解释等，奠定了坚实基础.

尽管小波变换在地面磁共振信号领域已有报道，但基于小波多尺度分解硬阈值方法有明显缺陷，损失部分有效信号，特别在尖峰噪声弛豫时间较长时，信号损失尤为严重，而本文提出的新的联合算法克服了小波多尺度分解硬阈值方法阈值固定的缺点，是小波技术在磁共振领域应用的提升，具有着重要的研究意义.

对本文进行仿真测试或实际测试，测试结果均说明了新算法的有效性，既可以去除尖峰噪声又可以通过阈值衰减相应的系数对地面磁共振缺失的信号进行重构，恢复原始信号.处理后得到的信号的特征参数误差很小可以满足实际应用的要求，通过一次处理完成噪声消除可以显著提高本文方法的计算效率.然而，在尖峰噪声持续40 ms或以上的数据(实际中较少发生)处理时，本文所提出的方法仍有改进空间，未来可以进一步完善阈值函数，提升本文算法的性能.