摘要:磁共振测深法(magnetic resonance sounding, MRS)具有无需钻探即可直接探测地下水含量的优势, 但是极低的信噪比(signal-to-noise ratio, SNR)限制了该方法的大范围应用, 目前的研究工作主要集中在消除MRS信号中的尖峰噪声和工频谐波噪声上, 而随机噪声由于其无规律性导致难以消除, 但是它的影响不容忽视. 目前MRS随机噪声的消减常采用叠加法, 需要重复采集数据叠加平均来达到消噪目的, 探测时间长且消噪效果有限. 针对这一问题, 本文提出了一种改进的短时傅里叶变换方法, 该方法通过处理单次采集的MRS包络信号来降低数据量, 提高数据处理效率. 改进的短时傅里叶变换方法采用解析信号代替常规短时傅里叶变换中的实值信号, 提高MRS信号时频域瞬时幅度的准确度, 得到MRS信号的高精度时频分布, 然后提取时频域峰值幅度和峰值相位重构信号来消除随机噪声. 仿真实验和实测数据处理结果表明, 该方法能够直接处理单次采集数据, 在信噪比高于–17.21 dB的情况下可有效提取MRS信号, 实现随机噪声的压制, 且与传统叠加法相比, 信噪比最多可提高27.88 dB, 均方根误差最多缩小36.44倍, 参数估计值更加准确. 本文的研究结果为利用MRS获取准确的地下水分布情况奠定了良好的基础.
关键词: 地面磁共振技术/
随机噪声/
改进短时傅里叶变换/
数据处理

English Abstract

全文HTML

--> --> -->

磁共振测深法(magnetic resonance sounding, MRS)是一种用于水文地质调查的非侵入式地球物理方法^[1-3], 由于其对质子具有独特敏感性, 因此可以直接定量地估计地下水含量和孔隙结构^[4-6]. 近年来, 这种方法在测井^[7]、实验室研究^[8]和野外测量^[9]中得到了快速发展. 但是MRS在地磁场(0.05—0.06 mT)中进行测量时, 微弱的MRS信号易被环境噪声淹没, 导致信噪比(signal-to-noise ratio, SNR)降低, 水文参数估计错误^[10,11].
为了提高测量结果的准确性, 通常情况下会采用信号增强和噪声抑制两种途径来提高信噪比. Davis等^[12,13]采用善于检测微弱信号的超导量子干涉器件(SQUID)代替传统接收线圈来提高磁共振仪器灵敏度; 2016年, Grunewald等^[9]通过在满足绝热条件时增加扳倒角来提高地下水中氢质子磁化强度的有效分量;2018年, Lin等^[14]在磁共振测深法中使用直流预极化场增强氢质子宏观磁矩, 测得了冰层下的流动水信号. 这些方法均在一定程度上提高了磁共振信号强度, 但是环境噪声的影响仍不容小觑. 干扰MRS信号的环境噪声主要包括尖峰噪声、工频谐波噪声和随机噪声, 在消噪过程中需要对这些噪声进行分类逐步消除, 常规步骤为: 1)去尖峰; 2)消除工频谐波噪声; 3)抑制随机噪声^[3,15]. 目前在MRS领域行之有效的降噪方法主要针对尖峰噪声和工频谐波噪声, 因为它们由特定噪声源产生, 具有一定的规律性, 可以总结出对应的数学表达式从而进行消除. Jiang等^[16]于2011年提出了一种基于罗曼诺夫斯基准则的统计叠加方法来丢弃尖峰噪声; Dalgaard等^[17]于2012年实现了基于非线性能量算子的尖峰噪声检测算法; Larsen^[18]于2016年的一项研究表明, 由电围网产生的尖峰噪声可建模为可调谐四阶带通滤波器的脉冲激励并进行消减. 陷波滤波法由Legchenko和Valla^[19]于2003年提出, 是消除单通道信号中工频谐波噪声的最常见并有效的方法. 之后, 基于Walsh^[20]于2008年发明的多通道仪器, 远程参考技术被应用于工频谐波噪声的消除^[17,20,21]. 当存在多个噪声源时, Larsen等^[22]于2014年提出了一种基于建模的降噪策略, 与多通道维纳滤波相结合来降低多源干扰. 然而随机噪声是一种来源广泛且不确定的前后独立平衡随机过程, 由于其无规律性、无法预测并且可能与有效信号发生频谱重叠而无法得到有效抑制. 目前最常见的随机噪声消除方法为叠加法^[2,17,23], 通过对同一地点多次采集信号进行叠加平均来提高信噪比, 但是该方法消耗时间长, 仪器探测效率低, 并且只能抑制部分随机噪声. 林婷婷等^[24]提出了一种分段时频峰值滤波法来处理磁共振全波信号中的随机噪声. 磁共振全波信号采样率高, 能够提供更加丰富的地下含水层信息, 但在实际探测过程中, 为了实现实时观测而对探测速度要求较高时, 仍然需要采集包络信号并对其中的随机噪声进行抑制.
为了简化MRS探测中的随机噪声消减方法, 实现随机噪声的即时处理, 本文提出了一种改进的短时傅里叶变换(modified short-time Fourier transform, MSTFT)方法, 应用于包含地下水探测关键参数的MRS包络信号中来消减随机噪声. 采集包络信号可以降低实测数据总量, 实现仪器的快速测量. 该方法首先对MRS包络信号作希尔伯特变换获取解析信号, 然后用解析信号代替MRS实值信号作对称短时傅里叶变换, 得到高精度的MRS信号时频分布, 最后提取峰值幅度和峰值相位重构信号, 消除随机噪声. 本文结构安排如下: 首先分析MRS信号中的随机噪声和传统叠加法消除随机噪声的局限性; 然后介绍本文提出的改进短时傅里叶变换算法的原理, 并阐述该算法提高MRS信号时频分布精度的原因; 最后通过仿真实验和处理实测数据, 验证该算法能够有效抑制随机噪声, 在不同信噪比下提取准确的信号参数.

原始MRS信号是由受交变磁场激发的氢质子在拉莫尔频率处震荡产生的指数衰减信号, 其计算公式为

$E\left( {t,q} \right) = {E_0}\left( q \right) \cdot {{\rm{e}}^{ - {t / {T_2^*\left( q \right)}}}} \cdot \cos \left( {2{\rm{\pi }}{f_{\rm{L}}}t + {\varphi _0}\left( q \right)} \right),$

其中t是采样时间, ${E_0}$是MRS信号初始振幅, 与含水量成正比, $T_2^*$为横向弛豫时间, 表征含水层介质孔隙度, 拉莫尔频率${f_{\rm{L}}}$由地磁场强度决定, 会随地理位置的变换而稍有改变, ${\varphi _0}$表示由导电介质和仪器相移产生的初始相位信息, q为激发脉冲矩.
通常情况下, 测量(1)式中的MRS全波信号需要高达几十kHz的采样频率, 记录的数据量非常大. 为了减少对数据采集的需求, 提高探测效率, 现有的仪器通常采用双通道正交探测系统辅以低通滤波器来探测MRS包络信号(JLMRS-I (Jiang et al. 2011), Numis^Plus和Numis^Lite (Bernard 2006))^[16,25].

$ \begin{split}\;& {E}_{\rm{envelope}}(t)\\=\;&{E}_{0}\cdot {\rm{e}}^{-t/{T}_{2}^{*}}\cdot {\rm{cos}}\left(2{\rm{\pi}}\Delta ft+{\varphi }_{0}\right)+{\sigma }_{x}\left(t\right)\\&+ {\rm{i}} \cdot \Big[{E}_{0}\cdot {\rm{e}}^{-t/{T}_{2}^{*}}\cdot {\rm{sin}}\left(2{\rm{\pi}}\Delta ft+{\varphi }_{0}\right) +{\sigma }_{y}\left(t\right)\Big],\\[-15pt] \end{split}$

其中, ${E_{{\rm{envelope}}}}\left( t \right)$表示MRS包络信号, $\Delta f$为拉莫尔频率和激发电场之间的频率偏移, 理想条件下$\Delta f$的值趋近于零. ${\sigma _x}\left( t \right)$和${\sigma _y}\left( t \right)$分别为MRS信号的同相噪声和异相噪声.
MRS信号中的随机噪声来源于宽带背景噪声和仪器电子元件噪声, 图1所示为对一组野外实测记录分别进行1, 4, 8和16次叠加法处理的结果. 常用的磁共振测深法工作方案为先测量空采噪声, 然后再发射电流采集信号, 因此该组实测记录包含噪声数据和含噪信号两部分, 分别如图1中灰色曲线和黑色曲线所示. 测量地点的拉莫尔频率为2330 Hz, 尖峰噪声和工频谐波噪声在使用叠加法之前已经消除. 如图1(a)和图1(b)所示, 1次叠加时信号的时间序列被噪声破坏, 频谱中几乎看不到信号的频率响应. 从图1(c)和图1(d)可以看出, 当叠加次数为4时, 随机噪声被抑制, 信号明显增强. 随着叠加次数$Ns$增大, 噪声水平越来越小. 8次叠加(图1(e)和图1(f))和16次叠加(图1(g)和图1(h))的结果对比表明, 虽然SNR随着叠加次数的增大而增大了, 但是随机噪声仍然存在.
图 1 叠加法分别叠加1, 4, 8和16次时消除随机噪声效果时域和频域图
Figure1. Time-domain and frequency-domain diagrams after stacking 1, 4, 8 and 16 times to suppress random noise, respectively.

假设S_i和σ_i分别为单次记录的信号和噪声, 那么经Ns次叠加后的信噪比为

${\rm{SNR}} = \frac{{\displaystyle\sum\nolimits_{i = 1}^{Ns} {{S_i}} }}{{\sqrt {\displaystyle\sum\nolimits_{i = 1}^{Ns} {\sigma _i^2} } }} = \frac{{Ns \cdot S}}{{\sqrt {Ns} \cdot \sigma }} = \sqrt {Ns} \cdot \frac{S}{\sigma }, $

其中$\displaystyle\sum\nolimits_{i = 1}^{Ns} {{S_i}}$和$\sqrt {\displaystyle\sum\nolimits_{i = 1}^{Ns} {\sigma _i^2} }$分别为$Ns$次采集的信号加和和噪声强度, S和σ分别为$Ns$个信号和噪声的平均幅度, ${S}/{\sigma }$为叠加前的信噪比, 因此在叠加$Ns$次后信噪比可以提高$\sqrt {Ns} $倍. 由指数函数的变化规律可知, 随着叠加次数的增多, SNR的增速逐渐减小并趋于平缓, 也就是说一直增加叠加次数是没有意义的. 另一方面, 多次采集数据需要耗费很多时间, 为了获得不同深度的水文地质信息, 通常需要发射16个不同的脉冲矩, 如果对每个脉冲矩时刻的数据进行100次测量, 再加上移动设备和铺设测量线圈的时间, 则总测量时间将超过100 h. 因此, 使用叠加法提高信噪比存在重复测量浪费时间且消噪效果有限的局限. 针对这一问题, 本文提出了一种MSTFT方法来代替传统的叠加法, 该方法可以处理单次采集数据压制随机噪声, 提高MRS信号的采集效率.

2

-->

3.1.短时傅里叶变换

实值信号$x\left( t \right)$的短时傅里叶变换形式为

$F\left( {\tau,f} \right) = \int_{ - \infty }^\infty {x\left( t \right)} w\left( {t - \tau } \right) \cdot {{\rm{e}}^{ - {\rm{i}}2{\rm{\pi }}ft}}{\rm{d}}t,$

其中, ${\rm{i}} = \sqrt { - 1} $, $\tau $为时移, f表示频率, 窗函数$w\left( t \right)$满足

$w\left( { - t} \right) = w\left( t \right),$

$F\left( {\tau, f} \right)$为信号$x\left( t \right)w\left( {t - \tau } \right)$的频谱, 包含信号$x\left( t \right)w\left( {t - \tau } \right)$的初始幅度和初始相位信息. 窗函数$w\left( t \right)$滑动截取信号, 分别在频率域形成以f为中心的窄带, 时间窗函数越短, 短时傅里叶变换对信号的时间分辨率越高, $x\left( t \right)w\left( {t - \tau } \right)$的初始幅度和初始相位可近似看作信号$x\left( t \right)$在固定窗内的瞬时幅度和瞬时相位, 则$F\left( {\tau, f} \right)$可表示信号$x\left( t \right)$随时间变化的频谱. 但是相应的频域分辨率会降低, 因此由Heisenberg不确定准则的限制, 窗函数$w\left( t \right)$的面积要不小于2^[26,27], 若定义

$F\left( {\tau,f} \right) = {A_{\rm{s}}}\left( {\tau,f} \right) \cdot {{\rm{e}}^{{\rm{i}}{\theta _{\rm{s}}}\left( {\tau,f} \right)}},$

则${A_{\rm{s}}}\left( {\tau, f} \right)$和${\theta _{\rm{s}}}\left( {\tau, f} \right)$可分别表示信号的时频幅度谱和时频相位谱.
2 -->

3.2.改进的短时傅里叶变换

为了准确地描述非平稳信号, 获得高精度的时频振幅谱和时频相位谱, 本文在STFT的基础上, 以解析信号代替实信号, 提出了一种MSTFT方法来提取MRS信号, 消除随机噪声, 其变换形式为

$G\left( {\tau,f} \right) = \int_{ - \infty }^\infty {h\left( t \right)} w\left( {t - \tau } \right) \cdot {{\rm{e}}^{ - {\rm{i}}2{\rm{\pi }}ft}}{\rm{d}}t,$

其中, $h\left( t \right)$为$x\left( t \right)$的解析信号, 由(8)式给出.

$ \left\{ \begin{aligned}&\hat{x}\left(t\right)=x\left(t\right)\ast \frac{1}{{\rm{\pi}}t}=\frac{1}{{\rm{\pi}}}{\displaystyle {\int }_{-\infty }^{\infty }\frac{x\left(\rho \right)}{t-\rho }}{\rm{d}}\rho, \\ &h\left(t\right)=x\left(t\right)+{\rm{i}}\cdot \widehat{x}\left(t\right).\end{aligned} \right.$

同样地, $G\left( {\tau, f} \right)$表示信号$ h\left( t \right)w\left( {t - \tau } \right) $的频谱, 如果定义

$G\left( {\tau,f} \right) = A\left( {\tau,f} \right) \cdot {{\rm{e}}^{{\rm{i}}\theta \left( {\tau,f} \right)}},$

$ \left\{ \begin{aligned}&A\left(\tau,{f}_{\tau }\right)={\rm{max}}\left[A\left(\tau,f\right)\right], \\ &\theta \left(\tau,{f}_{\tau }\right)={\rm{max}}\left[\theta \left(\tau,f\right)\right], \end{aligned} \right.$

则$A\left( {\tau, f} \right)$和$\theta \left( {\tau, f} \right)$分别表示信号的高精度时频幅度谱和高精度时频相位谱, $A\left( {\tau, {f_\tau }} \right)$和$\theta \left( {\tau, {f_\tau }} \right)$为给定时间$\tau $处信号频谱的峰值幅度和峰值相位, ${f_\tau }$为取得峰值幅度和峰值相位处的频率. 由于峰值幅度和峰值相位受随机噪声影响小, 因此可以用来重构信号:

$r\left( t \right) = {\rm{Re}} \left[ {G\left( {\tau,{f_\tau }} \right)} \right] = A\left( {\tau,{f_\tau }} \right)\cos \left[ {\theta \left( {\tau,{f_\tau }} \right)} \right],$

$r\left( t \right)$为抑制随机噪声影响后的MRS信号, 其中${\rm{Re}}\left[ {\,} \right]$表示选取复值信号的实数部分.
通常情况下傅里叶变换是在频域内计算的, 为了更好地理解MSTFT与常规短时傅里叶变换之间的区别, 将变换形式改写为频域形式:

$ \left\{ \begin{aligned}&G\left(\tau,f\right)={\displaystyle {\int }_{-\infty }^{\infty }H\left(\omega \right)}W\left(\omega -f\right)\cdot {\rm{e}}^{{\rm{i}}2{\rm{\pi}}\omega \tau }{\rm{d}}\omega, \\ &F\left(\tau,f\right)={\displaystyle {\int }_{-\infty }^{\infty }X\left(\omega \right)}W\left(\omega -f\right)\cdot {\rm{e}}^{{\rm{i}}2{\rm{\pi}}\omega \tau }{\rm{d}}\omega, \end{aligned} \right.$

其中, $H\left( \omega \right)$, $X\left( \omega \right)$和$W\left( \omega \right)$分别为$h\left( t \right)$, $x\left( t \right)$和$w\left( t \right)$的傅里叶变换形式, 且窗函数$W\left( \omega \right)$满足

$W\left( {\omega - f} \right) = \left\{ \begin{aligned} &{W_1}\left( \omega \right),~~~~\omega \geqslant 0 \\ &{W_2}\left( \omega \right),~~~~\omega < 0 .\end{aligned} \right.$

若以窗函数${W_1}\left( \omega \right)$代替$W\left( \omega \right)$, 则实值信号

$R\left( {\tau,f} \right) = \int_{ - \infty }^\infty {X\left( \omega \right)} {W_1}\left( \omega \right) \cdot {{\rm{e}}^{{\rm{i}}2{\rm{\pi }}\omega \tau }}{\rm{d}}\omega, $

表示信号$x\left( t \right)$的零相带通滤波响应, $G\left( {\tau, f} \right)$为$R\left( {\tau, f} \right)$的解析信号, 因此采用MSTFT计算得到的$G\left( {\tau, f} \right)$包含信号$R\left( {\tau, f} \right)$的瞬时幅度和瞬时相位信息. 但对于常规的短时傅里叶变换$F\left( {\tau, f} \right)$, 由Hilbert变换和Fourier变换的性质可知, 只有满足

$\left\{ \begin{aligned} &X\left( 0 \right){W_1}\left( 0 \right) = 0, \\ &{W_2}\left( 0 \right) = 0, \end{aligned} \right.$

时, 能得到

$\left\{ \begin{aligned} &G\left( {\tau,f} \right) = 2F\left( {\tau,f} \right), \\ &A\left( {\tau,f} \right) = 2{A_{\rm{s}}}\left( {\tau,f} \right), \\ &\theta \left( {\tau,f} \right) = {\theta _{\rm{s}}}\left( {\tau,f} \right). \end{aligned} \right.$

这时常规短时傅里叶变换能够准确获取信号$R\left( {\tau, f} \right)$的瞬时相位, 但是瞬时幅度仅为真实值的一半, 且在$X\left( 0 \right){W_1}\left( 0 \right) \ne 0$或${W_2}\left( 0 \right) \ne 0$时, (4)式所描述的短时傅里叶变换获取信号$R\left( {\tau, f} \right)$的瞬时幅度和瞬时相位均为近似值. 因此, 在使用相同时间窗函数的情况下, MSTFT比常规短时傅里叶变换获取的信号参数变化情况更为精确.
2 -->

3.3.信号实例

为了说明MSTFT与常规傅里叶变换之间的区别, 以及MSTFT消除随机噪声提取MRS信号的过程, 仿真一组受随机噪声干扰的MRS信号, 信号参数E₀ = 200 nV, $T_2^* $ = 150 ms, f_L = 2330 Hz, 信号持续时间500 ms. 加入的随机噪声均值为0, 标准差为50 nV, 使信噪比为0.51 dB. 以拉莫尔频率对信号重采样, 信号长度为1165. 选取对称窗函数为定义在$\left[ { - {L / {2, {L / 2}}}} \right]$区间内的海明窗, 其中L为窗长, 根据Heisenberg不确定准则和MRS信号特征, 选取窗长$L = 6$. 对包络信号和其解析信号分别进行对称窗短时傅里叶变换, 结果如图2所示. 图2(a)和图2(c)分别为时频幅度谱${A_{\rm{s}}}\left( {\tau, f} \right)$和高精度时频幅度谱$A\left( {\tau, f} \right)$, 从图中可知${A_{\rm{s}}}\left( {\tau, f} \right)$和$A\left( {\tau, f} \right)$的变化趋势基本一致, 但是$A\left( {\tau, f} \right)$的变化范围在0—300 nV之间, 与仿真信号的瞬时幅度一致, 而${A_{\rm{s}}}\left( {\tau, f} \right)$在0—150 nV之间变化. 图2(e)描述了${A_{\rm{s}}}\left( {\tau, f} \right)$和$A\left( {\tau, f} \right)$之间的差异, 与图2(a)的变化趋势相近, 说明$A\left( {\tau, f} \right)$约为${A_{\rm{s}}}\left( {\tau, f} \right)$的2倍. 图2(b)和图2(d)为时频相位谱${\theta _{\rm{s}}}\left( {\tau, f} \right)$和高精度时频相位谱$\theta \left( {\tau, f} \right)$, 二者之间的差异如图2(f)所示, 可以看出${\theta _{\rm{s}}}\left( {\tau, f} \right)$和$\theta \left( {\tau, f} \right)$的差异几乎为0, 与(16)式得出的结论相同.
图 2 MSTFT与常规短时傅里叶变换之间的区别　(a) 时频幅度谱${A_{\rm{s}}}\left( {\tau, f} \right)$; (b) 时频相位谱${\theta _{\rm{s}}}\left( {\tau, f} \right)$; (c) 高精度时频幅度谱$A\left( {\tau, f} \right)$; (d) 高精度时频相位谱$\theta \left( {\tau, f} \right)$; (e) ${A_{\rm{s}}}\left( {\tau, f} \right)$和$A\left( {\tau, f} \right)$的差异; (f) ${\theta _{\rm{s}}}\left( {\tau, f} \right)$和$\theta \left( {\tau, f} \right)$的差异
Figure2. Difference between the MSTFT and the conventional short-time Fourier transform: (a) time–frequency amplitude spectrum${A_{\rm{s}}}\left( {\tau, f} \right)$; (b) time–frequency phase spectrum${\theta _{\rm{s}}}\left( {\tau, f} \right)$; (c) high-precision time–frequency amplitude spectrum$A\left( {\tau, f} \right)$; (d) high-precision time–frequency phase spectrum$\theta \left( {\tau, f} \right)$; (e) difference between${A_{\rm{s}}}\left( {\tau, f} \right)$ and$A\left( {\tau, f} \right)$; (f) difference between${\theta _{\rm{s}}}\left( {\tau, f} \right)$ and$\theta \left( {\tau, f} \right)$.

图3(a)为含噪信号的时域波形, 其瞬时幅度(图3(b))和瞬时相位(图3(c))受随机噪声影响严重, 无法准确提取原信号信息. 采用前文所述的MSTFT方法获取高精度时频幅度谱$A\left( {\tau, f} \right)$和高精度时频相位谱$\theta \left( {\tau, f} \right)$后, 提取频谱峰值幅度$A\left( {\tau, {f_\tau }} \right) $和峰值相位$ \theta \left( {\tau, {f_\tau }} \right) $分别如图3(d)和图3(e)所示, 随机噪声的影响被消除, 信号的瞬时幅度展现出明显的指数衰减趋势, 瞬时相位维持在0°水平上. 重构MRS包络信号由图3(f)给出, 提取信号参数E₀ = 199.14 nV, $T_2^* $ = 149.2 ms, 信噪比提高为30.69 dB.
图 3 MSTFT方法消除随机噪声过程示例　(a) 仿真含噪数据; (b)瞬时幅度; (c) 瞬时相位; (d) 峰值幅度; (e) 峰值相位; (f)重构信号
Figure3. Example of the basic procedure for the random noise suppression using MSTFT: (a) Synthetic noisy signal; (b) instantaneous amplitude; (c) instantaneous phase; (d) peak amplitude; (e) peak phase; (f) reconstructed signal.

为了验证改进的短时傅里叶方法消除随机噪声, 提取MRS信号的有效性, 在仿真随机噪声和实测随机噪声中分别加入模拟的单指数衰减信号, 除了在波形上观察随机噪声的消噪效果外, 采用Müller-Petke等^[15]于2016提出的MRSmatlab软件中的同步检测方法提取信号参数, 并基于信噪比和均方根误差(root mean square error, RMSE)来量化评估本文所提方法的消噪效果. 均方根误差定义为

${\rm{RMSE}} = \sqrt {\frac{1}{k}{{\sum\limits_{k = 1}^K {\left( {r\left( k \right) - x\left( k \right)} \right)}^2 }}} ,$

其中K为信号长度, $x\left( k \right)$为理想信号, $r\left( k \right)$为重构信号.
2

-->

4.1.仿真噪声与仿真信号的消噪示例

单指数信号的仿真参数E₀ = 200 nV, $T_2^\ast$ = 150 ms, $\Delta f$ = 0 Hz, ${\varphi _0}$ = 60°, 信号持续时间为500 ms, 信号长度为1165. 为了使仿真的随机噪声更接近野外实际环境中的随机噪声, 根据(18)式进行仿真^[28]:

$\begin{split} {E_{{\rm{noise}}}}\left( k \right) =\;& {N_{\left( {0,1} \right)}}\left( k \right) \cdot x\left( k \right) \\ &\times \sqrt {{{\left( {\frac{{{V_{{\rm{bac}}}}\left( k \right)}}{{x\left( k \right)}}} \right)}^2} + V_{{\rm{uni}}}^2\left( k \right)} ,\end{split}$

其中, k为离散时间样本, ${N_{\left( {0, 1} \right)}}\left( k \right)$表示均值为0标准差为1 nV的高斯噪声, ${V_{{\rm{bac}}}}\left( m \right)$表示背景噪声, 考虑到结构噪声等非特定噪声的存在, ${V_{{\rm{uni}}}}\left( m \right)$为加入的均匀噪声.
图4给出了不同噪声条件下MSTFT方法对随机噪声的消除效果. 信号1中背景噪声${V_{{\rm{bac}}}}\left( m \right)$是均值为0标准差为50 nV的高斯噪声, 均匀噪声${V_{{\rm{uni}}}}\left( m \right)$是1%的理想信号值, 加入噪声后信噪比为0.68 dB. 然后应用MSTFT算法处理噪声, 图4(a)为含噪信号的高精度时频幅度谱, 噪声在时频域内随机分布, 但是对信号峰值影响较弱, 峰值出现在0 Hz附近, 幅度在300 nV以内. 沿频率轴提取各个时刻幅度和相位的最大值后重构信号, 信号时域和频域的处理结果分别如图4(b)和图4(c)所示, 其中灰色曲线表示含噪信号, 黑色曲线为采用MSTFT方法处理数据后重构的MRS信号, 红色曲线为仿真信号, 可以看出, 随机噪声成分被消除, 信号得以保留. 消噪之后提取的信号参数和SNR, RMSE如表1首行所列, 参数提取结果E₀ = (200.19 ± 5.01) nV, $T_2^* $ = (149.2 ± 2.8)ms, $\Delta f$ = (–0.03 ± 0.01) Hz, 信号SNR提升为32.67 dB, RMSE为(1.03 ± 0.55) nV. 为了测试该方法在较高噪声水平下的有效性, ${V_{{\rm{bac}}}}\left( m \right)$的标准差增大为100 nV, ${V_{{\rm{uni}}}}\left( m \right)$增大为理想信号的3%, 使信号2的信噪比为–5.20 dB. 表1的中间行为MSTFT方法对信号2的消噪结果, 提取参数E₀, $T_2^* $ 和$\Delta f$分别为(202.61 ± 7.90) nV, (152.0 ± 11.8) ms和(0.04 ± 0.03) Hz, SNR提高到24.01 dB, RMSE为(3.79 ± 1.89) nV. 最后, 加入干扰程度严重的随机噪声, ${V_{{\rm{bac}}}}\left( m \right)$标准差为200 nV, ${V_{{\rm{uni}}}}\left( m \right)$为5%理想信号, 信噪比降低至–11.22 dB. 从图4(h)的时间序列上看, 在随机噪声干扰较大的情况下, MSTFT方法仍然具有良好的消噪性能, 由图4(g)可以看出, 虽然随机噪声水平增大了, 但是由于其随机分布而不能在时频域产生聚集性, 因此对信号峰值处幅度和相位影响较小, 可以直接提取而不损耗信号信息. 消噪后MRS信号的参数提取结果E₀为(204.12 ± 14.07) nV, $T_2^* $ 为(154.4 ± 14.5) ms, $\Delta f$为(0.04 ± 0.06) Hz, 信噪比SNR和均方根误差RMSE分别为20.81 dB和(5.81 ± 2.42) nV.
图 4 3组仿真随机噪声消噪结果图　(a)低噪声强度下仿真数据高精度时频域振幅; (b) 低噪声强度下消噪结果时域图; (c) 低噪声强度下消噪结果频域图; (d) 中噪声强度下仿真数据高精度时频域振幅; (e) 中噪声强度下消噪结果时域图; (f) 中噪声强度下消噪结果频域图; (g) 高噪声强度下仿真数据高精度时频域振幅; (h) 高噪声强度下消噪结果时域图; (i) 高噪声强度下消噪结果频域图
Figure4. The de-noising results of 3 sets of random noise simulation: (a) High-precision time-frequency domain amplitude of simulated data under low noise intensity; (b) time domain results under low noise intensity; (c) frequency domain results under low noise intensity; (d) high-precision time-frequency domain amplitude of simulated data under moderate noise intensity; (e) time domain results under moderate noise intensity; (f) frequency domain results under moderate noise intensity; (g) high-precision time-frequency domain amplitude of simulated data under high noise intensity; (h) time domain results under high noise intensity; (i) frequency domain results under high noise intensity.

	E0/nV	$T_2^* $ /ms	$\Delta f$ / Hz	SNR/dB	RMSE/nV
信号1 (SNR = 0.68 dB)	200.19 ± 3.01	149.2 ± 2.8	–0.03 ± 0.01	32.67	1.03 ± 0.55
信号2 (SNR = –5.20 dB)	202.61 ± 4.90	152.0 ± 6.7	0.04 ± 0.03	24.01	3.79 ± 1.89
信号3 (SNR = –11.22 dB)	204.12 ± 5.96	154.4 ± 12.8	0.04 ± 0.06	20.81	5.81 ± 2.42

表13组仿真随机噪声消噪后参数估计情况表
Table1.The parameter estimation after de-noising of 3 sets of simulated random noise.

2 -->

4.2.实测噪声与仿真信号的消噪示例

为了进一步验证MSTFT方法的有效性, 本节进行了另一组仿真实验, 由(2)式给出的仿真信号中E₀ = 200 nV, $T_2^* $ = 100 ms, $\Delta f$ = 0 Hz, ${\varphi _0}$ = 57°, 信号持续时间为256 ms, 信号长度为596. 加入在两个不同地点进行MRS实验采集的两组实测噪声数据, 由于噪声记录中既包含尖峰噪声和工频谐波噪声, 又包含随机噪声, 因此先采用统计叠加法和自适应陷波法消除尖峰噪声和工频谐波噪声^[16], 只保留随机噪声, 并对随机噪声进行重采样, 使其与仿真信号长度相等. 然后采用MSTFT方法分别对两个测点数据中的一次噪声记录进行随机噪声压制, 并将结果与叠加法处理同一脉冲矩数据的结果进行对比. 由于测点2的噪声水平略高于测点1, 因此叠加法处理随机噪声时, 对测点1的每组脉冲矩数据叠加16次, 测点2的每组脉冲矩数据叠加32次.
图5给出了传统叠加法和改进短时傅里叶变换方法压制随机噪声的对比结果, 从第一列图中可以看出, 叠加法处理数据后的结果中仍残留随机噪声成分, 而图5中间列中, MSTFT方法处理单次噪声记录, 随机噪声抑制效果明显, 信号衰减趋势与仿真信号接近, 并且与传统叠加法相比, 该方法在较低信噪比下仍能提取所需信号.
图 5 叠加法和MSTFT方法消除随机噪声效果对比图　(a) 叠加法消除随机噪声时域图; (b) MSTFT方法消除随机噪声时域图; (c) 叠加法和MSTFT方法消除随机噪声频域对比图
Figure5. Comparison of the de-noising results by using stacking and MSTFT methods: (a) Results of random noise elimination by stacking in time domain; (b) results of random noise elimination by MSTFT in time domain; (c) comparison of the de-noising results by using stacking and MSTFT methods in frequency domain.

表2列出了实测数据中选定脉冲矩的单指数信号拟合参数结果, 并且计算了每种情况下信号的SNR和RMSE. 可以看出, 叠加法处理随机噪声后, 提取参数仍具有较大的拟合误差, 信噪比水平较低, 且每个信号T₂^* 参数的估计误差较大. 采用MSTFT方法处理单次记录数据提取的信号参数与仿真信号参数更为接近. 实验结果表明在信噪比高于–17.21 dB时均可有效提取信号, 与叠加法消噪结果相比, 信号的信噪比SNR最多可提高27.88 dB, 均方根误差RMSE最多缩小36.44倍, 信号拟合结果更精确.

		E0/nV	$T_2^* $ /ms	$\Delta f$/Hz	SNR/dB	RMSE/nV
测点1: q1	16次叠加	167.49	270.1	–0.53	–16.81	484.95
	MSTFT	197.42 ± 4.41	102.6 ± 5.3	–0.31 ± 0.09	11.07	13.31 ± 8.09
测点1: q2	16次叠加	163.78	285.5	–0.65	–16.85	487.52
	MSTFT	195.77 ± 5.02	103.9 ± 5.9	–0.34 ± 0.10	9.83	19.10 ± 10.83
测点2: q1	32次叠加	217.43	336.9	–0.74	–17.05	498.41
	MSTFT	205.20 ± 6.25	105.7 ± 8.7	–0.36 ± 0.11	8.08	24.31 ± 12.36
测点2: q2	32次叠加	221.51	359.4	–0.92	–17.21	507.79
	MSTFT	206.21 ± 7.47	107.4 ± 10.9	–0.37 ± 0.11	6.91	27.70 ± 18.64

表2叠加法和MSTFT消除随机噪声提取参数结果对比
Table2.Comparison of the parameter estimation after random noise elimination by stacking and MSTFT methods.

本节分别采用MSTFT方法和叠加法对同一组实测MRS数据进行处理, 对比两种方法下该组实测数据的消噪结果、参数估计和MRS反演情况. 所有数据均使用吉林大学自主开发的JLMRS-I系统在吉林省长春市郊区进行采集, 该地的地磁场强度为54720 nT, 拉莫尔频率为2330 Hz. 接收线圈为边长100 m的正方形, 16个发射脉冲矩分布在240—8350 A·ms范围内, 每组脉冲矩数据采集16次.
图6所示为不同脉冲矩下原始数据及经叠加法和MSTFT方法处理后的数据随时间的变化情况. 从图6(a)中可以看出, 原始数据受噪声干扰严重, 无法辨别信号的存在. 因此采用预处理方法预先消除尖峰噪声和工频谐波噪声, 只保留随机噪声, 并对含噪数据进行重采样. 然后采用传统叠加法消除随机噪声, 结果如图6(b)所示, 可以看出不同脉冲矩下信号幅度的变化情况和衰减趋势, 但未被完全消除的随机噪声仍会对信号产生影响. 采用MSTFT方法消除单次记录的预处理数据中的随机噪声, 结果如图6(c)所示, 随机噪声抑制效果更加明显, 信号更为突出.
图 6 实测数据处理结果图　(a) 野外探测原始数据; (b) 叠加法消除随机噪声后数据; (c) MSTFT方法消除随机噪声后数据
Figure6. De-noising results of field data: (a) Original field data; (b) random noise elimination by stacking; (c) random noise elimination by MSTFT.

图7为16个发射脉冲矩中任意4个脉冲矩信号的消噪结果. 其中黑色曲线为叠加法消噪后的MRS包络信号, 红色曲线为MSTFT方法处理单次数据记录后得到的MRS包络曲线. 图7(a)和图7(c)中可以观察到指数衰减趋势, 但是图7(e)和图7(g)中SNR仍然很低. 由时域和频域消噪结果可以看出, 与叠加法相比, 经MSTFT方法压制后信号中随机噪声水平更低. 并且如图7(f)和图7(h)的内置图所示, 当信噪比越低时, MSTFT方法的优势越明显. 另外, 虽然实际测量环境中的随机噪声不是纯白噪声数据, 但MSTFT方法仍然可行.
图 7 叠加法和MSTFT方法处理4个脉冲矩信号结果对比图　(a) 时域对比结果; (b) 频域对比结果
Figure7. Comparison of the de-noising results of 4 pulse moments by using stacking and MSTFT methods: (a) Time domain comparison; (b) frequency domain comparison.

16个脉冲矩信号的参数估计情况如图8所示. 可以看出, 采用传统叠加法和MSTFT方法处理数据后提取的信号初始振幅差距很小, 但是MSTFT方法消噪后提取的信号弛豫时间比叠加法消噪后提取的弛豫时间更为准确, 两种方法消噪后信号的频率偏移在零附近波动并且变化较为平滑, 因此可以认为采用MSTFT方法抑制随机噪声后MRS信号的提取结果更为准确.
图 8 叠加法和MSTFT处理实测数据提取参数结果对比　(a) 初始振幅; (b) 弛豫时间; (c) 频率偏移
Figure8. Parameter estimation after random noise elimination of field data by using stacking and MSTFT methods: (a) Initial amplitude; (b) relaxation time; (c) frequency offset.

对两种方法消噪后提取的MRS信号分别进行初始振幅反演, 反演结果和实验地点附近的钻孔结果对比如图9所示. 可以看出, 两组数据的反演结果得到的含水层信息均与钻孔结果显示的地层结构一致. 该地点存在两个含水层, 其中主含水层位于地下12.5—24 m的细-中-粗砂孔隙介质中, 最大含水量在11%—12%之间, 次含水层位于地下38—40 m的中砂孔隙介质中, 平均含水量约为4%. 但是相比于叠加法, 由MSTFT方法消噪后的数据进行反演, 反演结果得到的含水层信息与实际更为接近, 且实验数据采集中只需进行一次探测, 提高了野外工作效率.
图 9 叠加法和MSTFT处理实测数据提取MRS信号反演结果与钻孔结果对比　(a) 叠加法; (b) MSTFT方法
Figure9. Comparison of the inversion of the field data processd by stacking and MSTFT methods respectively with the borehole log: (a) Stacking result; (b) MSTFT result.

随机噪声由于其不确定性、无法预测而成为MRS信号中难以抑制的噪声种类之一. 在强噪声干扰下有效抑制随机噪声, 准确提取MRS信号是磁共振测深法能够广泛应用的基础. 针对现有随机噪声消除方法存在的探测效率低且只能抑制部分随机噪声的问题, 本文研究了一种改进的短时傅里叶变换方法, 通过对单次采集的MRS包络信号进行Hilbert变换求取解析信号, 然后代替实值信号进行短时傅里叶变换后提取时频峰值幅度和峰值相位重构信号, 实现MRS信号中随机噪声的有效抑制. 通过理论推导和分析数值仿真和实测数据处理结果, 可以得出以下结论:
1)在使用相同时间窗函数的情况下, 与常规短时傅里叶变换方法相比, MSTFT方法获取信号的瞬时相位不变, 但是瞬时幅度更为精确, 且由于时频域中峰值幅度和峰值相位受随机噪声影响小, 可以用来重构信号抑制随机噪声;
2)与传统叠加法相比, MSTFT方法可以在较低信噪比下提取所需信号, 且信噪比越低时, MSTFT方法的优势越明显, 仿真实验表明, MSTFT方法在信噪比高于–17.21 dB时可有效提取信号, 与叠加法相比, 消噪后信噪比最多可提高27.88 dB, 均方根误差RMSE最多缩小到原来的1/36.44. 另外, MSTFT方法处理单次记录数据提取的信号参数更为准确, 尤其能提高$T_2^* $参数估计的稳定性. 采用MSTFT方法消除随机噪声可以优化消噪效果, 提高探测效率, 为MRS方法快速获取地下水信息提供可行性.